Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4FTFAnnihilation.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27//
28
29// ------------------------------------------------------------
30// GEANT 4 class implemetation file
31//
32// ---------------- G4FTFAnnihilation --------------
33// by V. Uzhinsky, Spring 2011.
34// Take a projectile and a target
35// make annihilation or re-orangement of quarks and anti-quarks.
36// Ideas of Quark-Gluon-String model my A. Capella and A.B. Kaidalov
37// are implemented.
38// ---------------------------------------------------------------------
39
40#include "globals.hh"
41#include "Randomize.hh"
42#include "G4PhysicalConstants.hh"
43#include "G4SystemOfUnits.hh"
44
45#include "G4DiffractiveSplitableHadron.hh"
46#include "G4DiffractiveExcitation.hh"
47#include "G4FTFParameters.hh"
48#include "G4ElasticHNScattering.hh"
49#include "G4FTFAnnihilation.hh"
50
51#include "G4LorentzRotation.hh"
52#include "G4RotationMatrix.hh"
53#include "G4ThreeVector.hh"
54#include "G4ParticleDefinition.hh"
55#include "G4VSplitableHadron.hh"
56#include "G4ExcitedString.hh"
57#include "G4ParticleTable.hh"
58#include "G4Neutron.hh"
59#include "G4ParticleDefinition.hh"
60
61#include "G4Exp.hh"
62#include "G4Log.hh"
63#include "G4Pow.hh"
64
65//#include "UZHI_diffraction.hh"
66
67#include "G4ParticleTable.hh"
68
69//============================================================================
70
71//#define debugFTFannih
72
73
74//============================================================================
75
76G4FTFAnnihilation::G4FTFAnnihilation() {}
77
78
79//============================================================================
80
81G4FTFAnnihilation::~G4FTFAnnihilation() {}
82
83
84//============================================================================
85
86G4bool G4FTFAnnihilation::Annihilate( G4VSplitableHadron* projectile,
87 G4VSplitableHadron* target,
88 G4VSplitableHadron*& AdditionalString,
89 G4FTFParameters* theParameters ) const {
90
91 #ifdef debugFTFannih
92 G4cout << "---------------------------- Annihilation----------------" << G4endl;
93 #endif
94
95 CommonVariables common;
96
97 // Projectile parameters
98 common.Pprojectile = projectile->Get4Momentum();
99 G4int ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
100 if ( ProjectilePDGcode > 0 ) {
101 target->SetStatus( 3 );
102 return false;
103 }
104 G4double M0projectile2 = common.Pprojectile.mag2();
105
106 // Target parameters
107 G4int TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
108 common.Ptarget = target->Get4Momentum();
109 G4double M0target2 = common.Ptarget.mag2();
110
111 #ifdef debugFTFannih
112 G4cout << "PDG codes " << ProjectilePDGcode << " " << TargetPDGcode << G4endl
113 << "Pprojec " << common.Pprojectile << " " << common.Pprojectile.mag() << G4endl
114 << "Ptarget " << common.Ptarget << " " << common.Ptarget.mag() << G4endl
115 << "M0 proj target " << std::sqrt( M0projectile2 )
116 << " " << std::sqrt( M0target2 ) << G4endl;
117 #endif
118
119 // Kinematical properties of the interactions
120 G4LorentzVector Psum = common.Pprojectile + common.Ptarget; // 4-momentum in CMS
121 common.S = Psum.mag2();
122 common.SqrtS = std::sqrt( common.S );
123 #ifdef debugFTFannih
124 G4cout << "Psum SqrtS S " << Psum << " " << common.SqrtS << " " << common.S << G4endl;
125 #endif
126
127 // Transform momenta to cms and then rotate parallel to z axis
128 G4LorentzRotation toCms( -1*Psum.boostVector() );
129 G4LorentzVector Ptmp( toCms*common.Pprojectile );
130 toCms.rotateZ( -1*Ptmp.phi() );
131 toCms.rotateY( -1*Ptmp.theta() );
132 common.toLab = toCms.inverse();
133
134 if ( G4UniformRand() <= G4Pow::GetInstance()->powA( 1880.0/common.SqrtS, 4.0 ) ) {
135 common.RotateStrings = true;
136 common.RandomRotation.rotateZ( 2.0*pi*G4UniformRand() );
137 common.RandomRotation.rotateY( std::acos( 2.0*G4UniformRand() - 1.0 ) );
138 common.RandomRotation.rotateZ( 2.0*pi*G4UniformRand() ); //AR-Jun2018
139 }
140
141 G4double MesonProdThreshold = projectile->GetDefinition()->GetPDGMass() +
142 target->GetDefinition()->GetPDGMass() +
143 ( 2.0*140.0 + 16.0 )*MeV; // 2 Mpi + DeltaE
144 G4double Prel2 = sqr(common.S) + sqr(M0projectile2) + sqr(M0target2)
145 - 2.0*( common.S*(M0projectile2 + M0target2) + M0projectile2*M0target2 );
146 Prel2 /= common.S;
147 G4double X_a = 0.0, X_b = 0.0, X_c = 0.0, X_d = 0.0;
148 if ( Prel2 <= 0.0 ) {
149 // Annihilation at rest! Values are copied from Parameters
150 X_a = 625.1; // mb // 3-shirt diagram
151 X_b = 0.0; // mb // anti-quark-quark annihilation
152 X_c = 49.989; // mb // 2 Q-Qbar string creation
153 X_d = 6.614; // mb // One Q-Qbar string
154 #ifdef debugFTFannih
155 G4cout << "Annih at Rest X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d
156 << G4endl;
157 #endif
158 } else { // Annihilation in flight!
159 G4double FlowF = 1.0 / std::sqrt( Prel2 )*GeV;
160 // Process cross sections
161 X_a = 25.0*FlowF; // mb 3-shirt diagram
162 if ( common.SqrtS < MesonProdThreshold ) {
163 X_b = 3.13 + 140.0*G4Pow::GetInstance()->powA( ( MesonProdThreshold - common.SqrtS )/GeV, 2.5 );
164 } else {
165 X_b = 6.8*GeV / common.SqrtS; // mb anti-quark-quark annihilation
166 }
167 if ( projectile->GetDefinition()->GetPDGMass() + target->GetDefinition()->GetPDGMass()
168 > common.SqrtS ) {
169 X_b = 0.0;
170 }
171 // This can be in an interaction of low energy anti-baryon with off-shell nuclear nucleon
172 X_c = 2.0 * FlowF * sqr( projectile->GetDefinition()->GetPDGMass() +
173 target->GetDefinition()->GetPDGMass() ) / common.S; // mb re-arrangement of
174 // 2 quarks and 2 anti-quarks
175 X_d = 23.3*GeV*GeV / common.S; // mb anti-quark-quark string creation
176 #ifdef debugFTFannih
177 G4cout << "Annih in Flight X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d
178 << G4endl << "SqrtS MesonProdThreshold " << common.SqrtS << " " << MesonProdThreshold
179 << G4endl;
180 #endif
181 }
182
183 G4bool isUnknown = false;
184 if ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) { // Target proton or Delta+
185 if ( ProjectilePDGcode == -2212 || ProjectilePDGcode == -2214 ) { // anti_proton or anti_Delta+
186 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // Pbar P
187 } else if ( ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114 ) { // anti_neutron or anti_Delta0
188 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // NeutrBar P
189 } else if ( ProjectilePDGcode == -3122 ) { // anti_Lambda (no anti_Lambda* in PDG)
190 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar P
191 } else if ( ProjectilePDGcode == -3112 ) { // anti_Sigma- (no anti_Sigma*- in G4)
192 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma-Bar P
193 } else if ( ProjectilePDGcode == -3212 ) { // anti_Sigma0 (no anti_Sigma*0 in G4)
194 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar P
195 } else if ( ProjectilePDGcode == -3222 ) { // anti_Sigma+ (no anti_Sigma*+ in G4)
196 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma+Bar P
197 } else if ( ProjectilePDGcode == -3312 ) { // anti_Xi- (no anti_Xi*- in G4)
198 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi-Bar P
199 } else if ( ProjectilePDGcode == -3322 ) { // anti_Xi0 (no anti_Xi*0 in G4)
200 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi0Bar P
201 } else if ( ProjectilePDGcode == -3334 ) { // anti_Omega- (no anti_Omega*- in PDG)
202 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // Omega-Bar P
203 } else {
204 isUnknown = true;
205 }
206 } else if ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) { // Target neutron or Delta0
207 if ( ProjectilePDGcode == -2212 || ProjectilePDGcode == -2214 ) { // anti_proton or anti_Delta+
208 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // Pbar N
209 } else if ( ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114 ) { // anti_neutron or anti_Delta0
210 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // NeutrBar N
211 } else if ( ProjectilePDGcode == -3122 ) { // anti_Lambda (no anti_Lambda* in PDG)
212 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar N
213 } else if ( ProjectilePDGcode == -3112 ) { // anti_Sigma- (no anti_Sigma*- in G4)
214 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma-Bar N
215 } else if ( ProjectilePDGcode == -3212 ) { // anti_Sigma0 (no anti_Sigma*0 in G4)
216 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar N
217 } else if ( ProjectilePDGcode == -3222 ) { // anti_Sigma+ (no anti_Sigma*+ in G4)
218 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma+Bar N
219 } else if ( ProjectilePDGcode == -3312 ) { // anti_Xi- (no anti_Xi*- in G4)
220 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi-Bar N
221 } else if ( ProjectilePDGcode == -3322 ) { // anti_Xi0 (no anti_Xi*0 in G4)
222 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi0Bar N
223 } else if ( ProjectilePDGcode == -3334 ) { // anti_Omega- (no anti_Omega*- in PDG)
224 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // Omega-Bar N
225 } else {
226 isUnknown = true;
227 }
228 } else {
229 isUnknown = true;
230 }
231 if ( isUnknown ) {
232 G4cout << "Unknown anti-baryon for FTF annihilation: PDGcodes - "
233 << ProjectilePDGcode << " " << TargetPDGcode << G4endl;
234 }
235 #ifdef debugFTFannih
236 G4cout << "Annih Actual X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d << G4endl;
237 #endif
238
239 G4double Xannihilation = X_a + X_b + X_c + X_d;
240
241 // Projectile unpacking
242 UnpackBaryon( ProjectilePDGcode, common.AQ[0], common.AQ[1], common.AQ[2] );
243
244 // Target unpacking
245 UnpackBaryon( TargetPDGcode, common.Q[0], common.Q[1], common.Q[2] );
246
247 G4double Ksi = G4UniformRand();
248
249 if ( Ksi < X_a / Xannihilation ) {
250 return Create3QuarkAntiQuarkStrings( projectile, target, AdditionalString, theParameters, common );
251 }
252
253 G4int resultCode = 99;
254 if ( Ksi < (X_a + X_b) / Xannihilation ) {
255 resultCode = Create1DiquarkAntiDiquarkString( projectile, target, common );
256 if ( resultCode == 0 ) {
257 return true;
258 } else if ( resultCode == 99 ) {
259 return false;
260 }
261 }
262
263 if ( Ksi < ( X_a + X_b + X_c ) / Xannihilation ) {
264 resultCode = Create2QuarkAntiQuarkStrings( projectile, target, theParameters, common );
265 if ( resultCode == 0 ) {
266 return true;
267 } else if ( resultCode == 99 ) {
268 return false;
269 }
270 }
271
272 if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xannihilation ) {
273 return Create1QuarkAntiQuarkString( projectile, target, theParameters, common );
274 }
275
276 return true;
277}
278
279//-----------------------------------------------------------------------
280
281G4bool G4FTFAnnihilation::
282Create3QuarkAntiQuarkStrings( G4VSplitableHadron* projectile,
283 G4VSplitableHadron* target,
284 G4VSplitableHadron*& AdditionalString,
285 G4FTFParameters* theParameters,
286 G4FTFAnnihilation::CommonVariables& common ) const {
287 // Simulation of 3 anti-quark - quark strings creation
288
289 #ifdef debugFTFannih
290 G4cout << "Process a, 3 shirt diagram" << G4endl;
291 #endif
292
293 // Sampling of anti-quark order in projectile
294 G4int SampledCase = G4RandFlat::shootInt( G4long( 6 ) );
295 G4int Tmp1 = 0, Tmp2 = 0;
296 switch ( SampledCase ) {
297 case 1 : Tmp1 = common.AQ[1]; common.AQ[1] = common.AQ[2]; common.AQ[2] = Tmp1; break;
298 case 2 : Tmp1 = common.AQ[0]; common.AQ[0] = common.AQ[1]; common.AQ[1] = Tmp1; break;
299 case 3 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = common.AQ[2];
300 common.AQ[1] = Tmp1; common.AQ[2] = Tmp2; break;
301 case 4 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = Tmp2;
302 common.AQ[1] = common.AQ[2]; common.AQ[2] = Tmp1; break;
303 case 5 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = common.AQ[2];
304 common.AQ[1] = Tmp2; common.AQ[2] = Tmp1; break;
305 }
306
307 // Set the string properties
308 // An anti quark - quark pair can have the quantum number of either a scalar meson
309 // or a vector meson: the last digit of the PDG code is, respectively, 1 and 3.
310 // For simplicity only scalar is considered here.
311 G4int NewCode = 0, antiQuark = 0, quark = 0;
312 G4ParticleDefinition* TestParticle = nullptr;
313 for ( G4int iString = 0; iString < 3; ++iString ) { // Loop over the 3 string cases
314 if ( iString == 0 ) {
315 antiQuark = common.AQ[0]; quark = common.Q[0];
316 projectile->SplitUp();
317 projectile->SetFirstParton( antiQuark );
318 projectile->SetSecondParton( quark );
319 projectile->SetStatus( 0 );
320 } else if ( iString == 1 ) {
321 quark = common.Q[1]; antiQuark = common.AQ[1];
322 target->SplitUp();
323 target->SetFirstParton( quark );
324 target->SetSecondParton( antiQuark );
325 target->SetStatus( 0 );
326 } else { // iString == 2
327 antiQuark = common.AQ[2]; quark = common.Q[2];
328 }
329 G4int absAntiQuark = std::abs( antiQuark ), absQuark = std::abs( quark );
330 G4double aKsi = G4UniformRand();
331 if ( absAntiQuark == absQuark ) {
332 if ( absAntiQuark != 3 ) {
333 NewCode = 111; // Pi0-meson
334 if ( aKsi < 0.5 ) {
335 NewCode = 221; // Eta -meson
336 if ( aKsi < 0.25 ) {
337 NewCode = 331; // Eta'-meson
338 }
339 }
340 } else {
341 NewCode = 221; // Eta -meson
342 if ( aKsi < 0.5 ) {
343 NewCode = 331; // Eta'-meson
344 }
345 }
346 } else {
347 if ( absAntiQuark > absQuark ) {
348 NewCode = absAntiQuark*100 + absQuark*10 + 1; NewCode *= absAntiQuark/antiQuark;
349 } else {
350 NewCode = absQuark*100 + absAntiQuark*10 + 1; NewCode *= absQuark/quark;
351 }
352 }
353 if ( iString == 2 ) AdditionalString = new G4DiffractiveSplitableHadron();
354 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
355 if ( ! TestParticle ) return false;
356 if ( iString == 0 ) {
357 projectile->SetDefinition( TestParticle );
358 theParameters->SetProjMinDiffMass( 0.5 );
359 theParameters->SetProjMinNonDiffMass( 0.5 );
360 } else if ( iString == 1 ) {
361 target->SetDefinition( TestParticle );
362 theParameters->SetTarMinDiffMass( 0.5 );
363 theParameters->SetTarMinNonDiffMass( 0.5 );
364 } else { // iString == 2
365 AdditionalString->SetDefinition( TestParticle );
366 AdditionalString->SplitUp();
367 AdditionalString->SetFirstParton( common.AQ[2] );
368 AdditionalString->SetSecondParton( common.Q[2] );
369 AdditionalString->SetStatus( 0 );
370 }
371 } // End of the for loop over the 3 string cases
372
373 // Sampling kinematical properties:
374 // 1st string AQ[0]-Q[0], 2nd string AQ[1]-Q[1], 3rd string AQ[2]-Q[2]
375
376 G4ThreeVector Quark_Mom[6];
377 G4double ModMom2[6];
378 G4double AveragePt2 = 200.0*200.0, maxPtSquare = common.S, SumMt = 0.0, MassQ2 = 0.0,
379 ScaleFactor = 1.0;
380 G4int NumberOfTries = 0, loopCounter = 0;
381 const G4int maxNumberOfLoops = 1000;
382 do {
383 ++NumberOfTries;
384 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
385 // At large number of tries it would be better to reduce the values of <Pt^2>
386 ScaleFactor /= 2.0;
387 AveragePt2 *= ScaleFactor;
388 }
389 G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
390 for ( G4int i = 0; i < 6; ++i ) {
391 Quark_Mom [i] = GaussianPt( AveragePt2, maxPtSquare );
392 PtSum += Quark_Mom[i];
393 }
394 PtSum /= 6.0;
395 SumMt = 0.0;
396 for ( G4int i = 0; i < 6; ++i ) {
397 Quark_Mom[i] -= PtSum;
398 ModMom2[i] = Quark_Mom[i].mag2();
399 SumMt += std::sqrt( ModMom2[i] + MassQ2 );
400 }
401 } while ( ( SumMt > common.SqrtS ) &&
402 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
403 if ( loopCounter >= maxNumberOfLoops ) {
404 return false;
405 }
406
407 // Sampling X's of anti-baryon and baryon
408 G4double WminusTarget = 0.0, WplusProjectile = 0.0;
409 G4double Alfa_R = 0.5; ScaleFactor = 1.0;
410 G4bool Success = true;
411 NumberOfTries = 0; loopCounter = 0;
412 do {
413 Success = true;
414 ++NumberOfTries;
415 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
416 // At large number of tries it would be better to reduce the values of Pt's
417 ScaleFactor /= 2.0;
418 }
419 G4double Alfa = 0.0, Beta = 0.0;
420 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // anti-baryon (1st case), baryon (2nd case)
421 G4double x1 = 0.0, x2 = 0.0;
422 G4double r1 = G4UniformRand(), r2 = G4UniformRand();
423 if ( Alfa_R == 1.0 ) {
424 x1 = 1.0 - std::sqrt( r1 );
425 x2 = (1.0 - x1) * r2;
426 } else {
427 x1 = sqr( r1 );
428 x2 = (1.0 - x1) * sqr( std::sin( pi/2.0*r2 ) );
429 }
430 G4double x3 = 1.0 - x1 - x2;
431 G4int index = iCase*3; // 0 for anti-baryon, 3 for baryon
432 Quark_Mom[index].setZ( x1 ); Quark_Mom[index+1].setZ( x2 ); Quark_Mom[index+2].setZ( x3 );
433 for ( G4int i = 0; i < 3; ++i ) { // Loop over the 3 (anti-)quarks
434 if ( Quark_Mom[index+i].getZ() != 0.0 ) {
435 G4double val = ( ScaleFactor * ModMom2[index+i] + MassQ2 ) / Quark_Mom[index+i].getZ();
436 if ( iCase == 0 ) { // anti-baryon
437 Alfa += val;
438 } else { // baryon (iCase == 1)
439 Beta += val;
440 }
441 } else {
442 Success = false;
443 }
444 }
445 }
446 if ( ! Success ) continue;
447 if ( std::sqrt( Alfa ) + std::sqrt( Beta ) > common.SqrtS ) {
448 Success = false;
449 continue;
450 }
451 G4double DecayMomentum2 = sqr(common.S) + sqr(Alfa) + sqr(Beta)
452 - 2.0*( common.S*(Alfa + Beta) + Alfa*Beta );
453 WminusTarget = ( common.S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.SqrtS;
454 WplusProjectile = common.SqrtS - Beta/WminusTarget;
455 } while ( ( ! Success ) &&
456 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
457 if ( loopCounter >= maxNumberOfLoops ) {
458 return false;
459 }
460
461 G4double SqrtScaleF = std::sqrt( ScaleFactor );
462 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // anti-baryon (1st case), baryon (2nd case)
463 G4int index = iCase*3; // 0 for anti-baryon, 3 for baryon
464 G4double w = WplusProjectile; // for anti-baryon
465 if ( iCase == 1 ) w = - WminusTarget; // for baryon
466 for ( G4int i = 0; i < 3; ++i ) {
467 G4double Pz = w * Quark_Mom[index+i].getZ() / 2.0 -
468 ( ScaleFactor * ModMom2[index+i] + MassQ2 ) /
469 ( 2.0 * w * Quark_Mom[index+i].getZ() );
470 Quark_Mom[index+i].setZ( Pz );
471 if ( ScaleFactor != 1.0 ) {
472 Quark_Mom[index+i].setX( SqrtScaleF * Quark_Mom[index+i].getX() );
473 Quark_Mom[index+i].setY( SqrtScaleF * Quark_Mom[index+i].getY() );
474 }
475 }
476 }
477 G4LorentzVector Pstring1, Pstring2, Pstring3;
478 G4double YstringMax = 0.0, YstringMin = 0.0;
479 for ( G4int i = 0; i < 3; ++i ) {
480 G4ThreeVector tmp = Quark_Mom[i] + Quark_Mom[i+3];
481 G4LorentzVector Pstring( tmp, std::sqrt( Quark_Mom[i].mag2() + MassQ2 ) +
482 std::sqrt( Quark_Mom[i+3].mag2() + MassQ2 ) );
483 //AR-Jun2018 if ( common.RotateStrings ) Pstring *= common.RandomRotation;
484 // Add protection for rapidity = 0.5*ln( (E+Pz)/(E-Pz) )
485 G4double Ystring = 0.0;
486 if ( Pstring.e() > 1.0e-30 ) {
487 if ( Pstring.e() + Pstring.pz() < 1.0e-30 ) { // Very small numerator in the logarithm
488 Ystring = -1.0e30; // A very large negative value (E ~= -Pz)
489 if ( Pstring.e() - Pstring.pz() < 1.0e-30 ) { // Very small denominator in the logarithm
490 Ystring = 1.0e30; // A very large positive value (E ~= Pz)
491 } else { // Normal case
492 Ystring = Pstring.rapidity();
493 }
494 }
495 }
496 // Keep ordering in rapidity: "1" highest, "2" middle, "3" smallest
497 if ( i == 0 ) {
498 Pstring1 = Pstring; YstringMax = Ystring;
499 } else if ( i == 1 ) {
500 if ( Ystring > YstringMax ) {
501 Pstring2 = Pstring1; YstringMin = YstringMax;
502 Pstring1 = Pstring; YstringMax = Ystring;
503 } else {
504 Pstring2 = Pstring; YstringMin = Ystring;
505 }
506 } else { // i == 2
507 if ( Ystring > YstringMax ) {
508 Pstring3 = Pstring2;
509 Pstring2 = Pstring1;
510 Pstring1 = Pstring;
511 } else if ( Ystring > YstringMin ) {
512 Pstring3 = Pstring2;
513 Pstring2 = Pstring;
514 } else {
515 Pstring3 = Pstring;
516 }
517 }
518 }
519 common.Pprojectile = Pstring1; // Highest rapidity
520 common.Ptarget = Pstring3; // Lowest rapidity
521 G4LorentzVector LeftString( Pstring2 ); // Middle rapidity
522
523 if ( common.RotateStrings ) { //AR-Jun2018
524 common.Pprojectile *= common.RandomRotation;
525 LeftString *= common.RandomRotation;
526 common.Ptarget *= common.RandomRotation;
527 }
528
529 common.Pprojectile.transform( common.toLab );
530 LeftString.transform( common.toLab );
531 common.Ptarget.transform( common.toLab );
532
533 // Calculation of the creation time
534 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
535 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
536 projectile->SetPosition( target->GetPosition() );
537 AdditionalString->SetTimeOfCreation( target->GetTimeOfCreation() );
538 AdditionalString->SetPosition( target->GetPosition() );
539
540 projectile->Set4Momentum( common.Pprojectile );
541 AdditionalString->Set4Momentum( LeftString );
542 target->Set4Momentum( common.Ptarget );
543 projectile->IncrementCollisionCount( 1 );
544 AdditionalString->IncrementCollisionCount( 1 );
545 target->IncrementCollisionCount( 1 );
546
547 return true;
548}
549
550//-----------------------------------------------------------------------
551
552G4int G4FTFAnnihilation::
553Create1DiquarkAntiDiquarkString( G4VSplitableHadron* projectile,
554 G4VSplitableHadron* target,
555 G4FTFAnnihilation::CommonVariables& common ) const {
556 // Simulation of anti-diquark-diquark string creation.
557 // This method returns an integer code - instead of a boolean, with the following meaning:
558 // "0" : successfully ended and nothing else needs to be done;
559 // "1" : successfully completed, but the work needs to be continued;
560 // "99" : unsuccessfully ended, nothing else can be done.
561
562 #ifdef debugFTFannih
563 G4cout << "Process b, quark - anti-quark annihilation, di-q - anti-di-q string" << G4endl;
564 #endif
565
566 G4int CandidatsN = 0, CandAQ[9][2] = {}, CandQ[9][2] = {};
567 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { // index of the 3 constituent anti-quarks of the antibaryon projectile
568 for ( G4int iQ = 0; iQ < 3; ++iQ ) { // index of the 3 constituent quarks of the nucleon target
569 if ( -common.AQ[iAQ] == common.Q[iQ] ) { // antiquark - quark that can annihilate
570 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
571 // of the (anti-baryon) projectile or (nucleon) target.
572 if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; }
573 if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; }
574 if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; }
575 if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; }
576 if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; }
577 if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; }
578 ++CandidatsN;
579 }
580 }
581 }
582
583 // Remaining two (anti-)quarks that form the (anti-)diquark
584 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, LeftQ2 = 0;
585 if ( CandidatsN != 0 ) {
586 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
587 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0] ];
588 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1] ];
589 LeftQ1 = common.Q[ CandQ[SampledCase][0] ];
590 LeftQ2 = common.Q[ CandQ[SampledCase][1] ];
591
592 // Build anti-diquark and diquark : the last digit can be either 3 - for all combinations
593 // of anti-quark - anti-quark and quark - quark - or 1 - only when the two anti-quarks
594 // or quarks are different. For simplicity, only 3 is considered.
595 G4int Anti_DQ = 0, DQ = 0;
596 if ( std::abs( LeftAQ1 ) > std::abs( LeftAQ2 ) ) {
597 Anti_DQ = 1000*LeftAQ1 + 100*LeftAQ2 - 3;
598 } else {
599 Anti_DQ = 1000*LeftAQ2 + 100*LeftAQ1 - 3;
600 }
601 if ( std::abs( LeftQ1 ) > std::abs( LeftQ2 ) ) {
602 DQ = 1000*LeftQ1 + 100*LeftQ2 + 3;
603 } else {
604 DQ = 1000*LeftQ2 + 100*LeftQ1 + 3;
605 }
606
607 // Set the string properties
608 projectile->SplitUp();
609 projectile->SetFirstParton( DQ );
610 projectile->SetSecondParton( Anti_DQ );
611
612 if ( common.RotateStrings ) {
613 G4LorentzVector Pquark = G4LorentzVector( 0.0, 0.0, common.SqrtS/2.0, common.SqrtS/2.0 );
614 Pquark *= common.RandomRotation;
615 G4LorentzVector Paquark = G4LorentzVector( 0.0, 0.0, -common.SqrtS/2.0, common.SqrtS/2.0 );
616 Paquark *= common.RandomRotation;
617 Pquark.transform( common.toLab ); projectile->GetNextParton()->Set4Momentum( Pquark );
618 Paquark.transform( common.toLab ); projectile->GetNextAntiParton()->Set4Momentum( Paquark );
619 }
620
621 projectile->SetStatus( 0 );
622 target->SetStatus( 4 ); // The target nucleon has annihilated 3->4
623 common.Pprojectile.setPx( 0.0 );
624 common.Pprojectile.setPy( 0.0 );
625 common.Pprojectile.setPz( 0.0 );
626 common.Pprojectile.setE( common.SqrtS );
627 common.Pprojectile.transform( common.toLab );
628
629 // Calculation of the creation time
630 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
631 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
632 projectile->SetPosition( target->GetPosition() );
633 projectile->Set4Momentum( common.Pprojectile );
634 projectile->IncrementCollisionCount( 1 );
635 target->IncrementCollisionCount( 1 );
636
637 return 0; // Completed successfully: nothing else to be done
638 } // End of if ( CandidatsN != 0 )
639
640 return 1; // Successfully ended, but the work is not over
641}
642
643//-----------------------------------------------------------------------
644
645G4int G4FTFAnnihilation::
646Create2QuarkAntiQuarkStrings( G4VSplitableHadron* projectile,
647 G4VSplitableHadron* target,
648 G4FTFParameters* theParameters,
649 G4FTFAnnihilation::CommonVariables& common ) const {
650 // Simulation of 2 anti-quark-quark strings creation.
651 // This method returns an integer code - instead of a boolean, with the following meaning:
652 // "0" : successfully ended and nothing else needs to be done;
653 // "1" : successfully completed, but the work needs to be continued;
654 // "99" : unsuccessfully ended, nothing else can be done.
655
656 #ifdef debugFTFannih
657 G4cout << "Process c, quark - anti-quark and string junctions annihilation, 2 strings left."
658 << G4endl;
659 #endif
660
661 G4int CandidatsN = 0, CandAQ[9][2] = {}, CandQ[9][2] = {};
662 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, LeftQ2 = 0;
663 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { // index of the 3 constituent anti-quarks of the antibaryon projectile
664 for ( G4int iQ = 0; iQ < 3; ++iQ ) { // index of the 3 constituent quarks of the nucleon target
665 if ( -common.AQ[iAQ] == common.Q[iQ] ) { // antiquark - quark that can annihilate
666 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
667 // of the (anti-baryon) projectile or (nucleon) target.
668 if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; }
669 if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; }
670 if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; }
671 if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; }
672 if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; }
673 if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; }
674 ++CandidatsN;
675 }
676 }
677 }
678
679 if ( CandidatsN != 0 ) {
680 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
681 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0] ];
682 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1] ];
683 if ( G4UniformRand() < 0.5 ) {
684 LeftQ1 = common.Q[ CandQ[SampledCase][0] ];
685 LeftQ2 = common.Q[ CandQ[SampledCase][1] ];
686 } else {
687 LeftQ2 = common.Q[ CandQ[SampledCase][0] ];
688 LeftQ1 = common.Q[ CandQ[SampledCase][1] ];
689 }
690
691 // Set the string properties
692 // An anti quark - quark pair can have the quantum number of either a scalar meson
693 // or a vector meson: the last digit of the PDG code is, respectively, 1 and 3.
694 // For simplicity only scalar is considered here.
695 G4int NewCode = 0, antiQuark = 0, quark = 0;
696 G4ParticleDefinition* TestParticle = nullptr;
697 for ( G4int iString = 0; iString < 2; ++iString ) { // Loop over the 2 string cases
698 if ( iString == 0 ) {
699 antiQuark = LeftAQ1; quark = LeftQ1;
700 projectile->SplitUp();
701 projectile->SetFirstParton( antiQuark );
702 projectile->SetSecondParton( quark );
703 projectile->SetStatus( 0 );
704 } else { // iString == 1
705 quark = LeftQ2; antiQuark = LeftAQ2;
706 target->SplitUp();
707 target->SetFirstParton( quark );
708 target->SetSecondParton( antiQuark );
709 target->SetStatus( 0 );
710 }
711 G4int absAntiQuark = std::abs( antiQuark ), absQuark = std::abs( quark );
712 G4double aKsi = G4UniformRand();
713 if ( absAntiQuark == absQuark ) {
714 if ( absAntiQuark != 3 ) {
715 NewCode = 111; // Pi0-meson
716 if ( aKsi < 0.5 ) {
717 NewCode = 221; // Eta -meson
718 if ( aKsi < 0.25 ) {
719 NewCode = 331; // Eta'-meson
720 }
721 }
722 } else {
723 NewCode = 221; // Eta -meson
724 if ( aKsi < 0.5 ) {
725 NewCode = 331; // Eta'-meson
726 }
727 }
728 } else {
729 if ( absAntiQuark > absQuark ) {
730 NewCode = absAntiQuark*100 + absQuark*10 + 1; NewCode *= absAntiQuark/antiQuark;
731 } else {
732 NewCode = absQuark*100 + absAntiQuark*10 + 1; NewCode *= absQuark/quark;
733 }
734 }
735 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
736 if ( ! TestParticle ) return 99; // unsuccessfully ended, nothing else can be done
737 if ( iString == 0 ) {
738 projectile->SetDefinition( TestParticle );
739 theParameters->SetProjMinDiffMass( 0.5 );
740 theParameters->SetProjMinNonDiffMass( 0.5 );
741 } else { // iString == 1
742 target->SetDefinition( TestParticle );
743 theParameters->SetTarMinDiffMass( 0.5 );
744 theParameters->SetTarMinNonDiffMass( 0.5 );
745 }
746 } // End of loop over the 2 string cases
747
748 // Sampling kinematical properties: 1st string LeftAQ1-LeftQ1, 2nd string LeftAQ2-LeftQ2
749 G4ThreeVector Quark_Mom[4];
750 G4double ModMom2[4];
751 G4double AveragePt2 = 200.0*200.0, maxPtSquare = common.S, SumMt = 0.0, MassQ2 = 0.0,
752 ScaleFactor = 1.0;
753 G4int NumberOfTries = 0, loopCounter = 0;
754 const G4int maxNumberOfLoops = 1000;
755 do {
756 ++NumberOfTries;
757 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
758 // At large number of tries it would be better to reduce the values of <Pt^2>
759 ScaleFactor /= 2.0;
760 AveragePt2 *= ScaleFactor;
761 }
762 G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
763 for( G4int i = 0; i < 4; ++i ) {
764 Quark_Mom[i] = GaussianPt( AveragePt2, maxPtSquare );
765 PtSum += Quark_Mom[i];
766 }
767 PtSum /= 4.0;
768 SumMt = 0.0;
769 for ( G4int i = 0; i < 4; ++i ) {
770 Quark_Mom[i] -= PtSum;
771 ModMom2[i] = Quark_Mom[i].mag2();
772 SumMt += std::sqrt( ModMom2[i] + MassQ2 );
773 }
774 } while ( ( SumMt > common.SqrtS ) &&
775 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
776 if ( loopCounter >= maxNumberOfLoops ) {
777 return 99; // unsuccessfully ended, nothing else can be done
778 }
779
780 // Sampling X's of the two strings
781 G4double WminusTarget = 0.0, WplusProjectile = 0.0, Alfa_R = 0.5; ScaleFactor = 1.0;
782 G4bool Success = true;
783 NumberOfTries = 0, loopCounter = 0;
784 do {
785 Success = true;
786 ++NumberOfTries;
787 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
788 // At large number of tries it would be better to reduce the values of Pt's
789 ScaleFactor /= 2.0;
790 }
791 G4double Alfa = 0.0, Beta = 0.0;
792 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // Loop over the two strings
793 G4double x = 0.0, r = G4UniformRand();
794 if ( Alfa_R == 1.0 ) {
795 if ( iCase == 0 ) { // first string
796 x = std::sqrt( r );
797 } else { // second string
798 x = 1.0 - std::sqrt( r);
799 }
800 } else {
801 x = sqr( std::sin( pi/2.0*r ) );
802 }
803 G4int index = iCase*2; // 0 for the first string, 2 for the second string
804 Quark_Mom[index].setZ( x ); Quark_Mom[index+1].setZ( 1.0 - x );
805 for ( G4int i = 0; i < 2; ++i ) {
806 if ( Quark_Mom[i].getZ() != 0.0 ) {
807 G4double val = ( ScaleFactor * ModMom2[index+i] + MassQ2 ) / Quark_Mom[index+i].getZ();
808 if ( iCase == 0 ) { // first string
809 Alfa += val;
810 } else { // second string
811 Beta += val;
812 }
813 } else {
814 Success = false;
815 }
816 }
817 }
818 if ( ! Success ) continue;
819 if ( std::sqrt( Alfa ) + std::sqrt( Beta ) > common.SqrtS ) {
820 Success = false;
821 continue;
822 }
823 G4double DecayMomentum2 = sqr(common.S) + sqr(Alfa) + sqr(Beta)
824 - 2.0*( common.S*(Alfa + Beta) + Alfa*Beta );
825 WminusTarget = ( common.S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.SqrtS;
826 WplusProjectile = common.SqrtS - Beta/WminusTarget;
827 } while ( ( ! Success ) &&
828 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
829 if ( loopCounter >= maxNumberOfLoops ) {
830 return 99; // unsuccessfully ended, nothing else can be done
831 }
832
833 G4double SqrtScaleF = std::sqrt( ScaleFactor );
834 G4LorentzVector Pstring1, Pstring2;
835 G4double Ystring1 = 0.0, Ystring2 = 0.0;
836 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // Loop over the two strings
837 G4int index = iCase*2; // 0 for the first string, 2 for the second string
838 for ( G4int i = 0; i < 2; ++i ) {
839 G4double w = WplusProjectile; // For the first string
840 if ( iCase == 1 ) w = - WminusTarget; // For the second string
841 G4double Pz = w * Quark_Mom[index+i].getZ() / 2.0
842 - ( ScaleFactor * ModMom2[index+i] + MassQ2 ) /
843 ( 2.0 * w * Quark_Mom[index+i].getZ() );
844 Quark_Mom[index+i].setZ( Pz );
845 if ( ScaleFactor != 1.0 ) {
846 Quark_Mom[index+i].setX( SqrtScaleF * Quark_Mom[index+i].getX() );
847 Quark_Mom[index+i].setY( SqrtScaleF * Quark_Mom[index+i].getY() );
848 }
849 }
850 }
851 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // Loop over the two strings
852 G4ThreeVector tmp = Quark_Mom[iCase] + Quark_Mom[iCase+2];
853 G4LorentzVector Pstring( tmp, std::sqrt( Quark_Mom[iCase].mag2() + MassQ2 ) +
854 std::sqrt( Quark_Mom[iCase+2].mag2() + MassQ2 ) );
855 //AR-Jun2018 if ( common.RotateStrings ) Pstring *= common.RandomRotation;
856 // Add protection for rapidity = 0.5*ln( (E+Pz)/(E-Pz) )
857 G4double Ystring = 0.0;
858 if ( Pstring.e() > 1.0e-30 ) {
859 if ( Pstring.e() + Pstring.pz() < 1.0e-30 ) { // Very small numerator in the logarithm
860 Ystring = -1.0e30; // A very large negative value (E ~= -Pz)
861 if ( Pstring.e() - Pstring.pz() < 1.0e-30 ) { // Very small denominator in the logarithm
862 Ystring = 1.0e30; // A very large positive value (E ~= Pz)
863 } else { // Normal case
864 Ystring = Pstring.rapidity();
865 }
866 }
867 }
868 if ( iCase == 0 ) { // For the first string
869 Pstring1 = Pstring; Ystring1 = Ystring;
870 } else { // For the second string
871 Pstring2 = Pstring; Ystring2 = Ystring;
872 }
873 }
874 if ( Ystring1 > Ystring2 ) {
875 common.Pprojectile = Pstring1; common.Ptarget = Pstring2;
876 } else {
877 common.Pprojectile = Pstring2; common.Ptarget = Pstring1;
878 }
879
880 if ( common.RotateStrings ) { //AR-Jun2018
881 common.Pprojectile *= common.RandomRotation;
882 common.Ptarget *= common.RandomRotation;
883 }
884
885 common.Pprojectile.transform( common.toLab );
886 common.Ptarget.transform( common.toLab );
887
888 // Calculation of the creation time
889 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
890 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
891 projectile->SetPosition( target->GetPosition() );
892 projectile->Set4Momentum( common.Pprojectile );
893 target->Set4Momentum( common.Ptarget );
894 projectile->IncrementCollisionCount( 1 );
895 target->IncrementCollisionCount( 1 );
896
897 return 0; // Completed successfully: nothing else to be done
898 } // End of if ( CandidatsN != 0 )
899
900 return 1; // Successfully ended, but the work is not over
901}
902
903//-----------------------------------------------------------------------
904
905G4bool G4FTFAnnihilation::
906Create1QuarkAntiQuarkString( G4VSplitableHadron* projectile,
907 G4VSplitableHadron* target,
908 G4FTFParameters* theParameters,
909 G4FTFAnnihilation::CommonVariables& common ) const {
910 // Simulation of anti-quark - quark string creation
911
912 #ifdef debugFTFannih
913 G4cout << "Process d, only 1 quark - anti-quark string" << G4endl;
914 #endif
915
916 // Determine the set of candidates anti-quark - quark pairs that do not annihilate.
917 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
918 // of the (anti-baryon) projectile or (nucleon) target.
919 G4int CandidatsN = 0, CandAQ[36], CandQ[36];
920 G4int LeftAQ = 0, LeftQ = 0;
921 for ( G4int iAQ1 = 0; iAQ1 < 3; ++iAQ1 ) {
922 for ( G4int iAQ2 = 0; iAQ2 < 3; ++iAQ2 ) {
923 if ( iAQ1 != iAQ2 ) {
924 for ( G4int iQ1 = 0; iQ1 < 3; ++iQ1 ) {
925 for ( G4int iQ2 = 0; iQ2 < 3; ++iQ2 ) {
926 if ( iQ1 != iQ2 ) {
927 if ( -common.AQ[iAQ1] == common.Q[iQ1] && -common.AQ[iAQ2] == common.Q[iQ2] ) {
928 if ( ( iAQ1 == 0 && iAQ2 == 1 ) || ( iAQ1 == 1 && iAQ2 == 0 ) ) {
929 CandAQ[CandidatsN] = 2;
930 } else if ( ( iAQ1 == 0 && iAQ2 == 2 ) || ( iAQ1 == 2 && iAQ2 == 0 ) ) {
931 CandAQ[CandidatsN] = 1;
932 } else if ( ( iAQ1 == 1 && iAQ2 == 2 ) || ( iAQ1 == 2 && iAQ2 == 1 ) ) {
933 CandAQ[CandidatsN] = 0;
934 }
935 if ( ( iQ1 == 0 && iQ2 == 1 ) || ( iQ1 == 1 && iQ2 == 0 ) ) {
936 CandQ[CandidatsN] = 2;
937 } else if ( ( iQ1 == 0 && iQ2 == 2 ) || ( iQ1 == 2 && iQ2 == 0 ) ) {
938 CandQ[CandidatsN] = 1;
939 } else if ( ( iQ1 == 1 && iQ2 == 2 ) || ( iQ1 == 2 && iQ2 == 1 ) ) {
940 CandQ[CandidatsN] = 0;
941 }
942 ++CandidatsN;
943 }
944 }
945 }
946 }
947 }
948 }
949 }
950
951 if ( CandidatsN != 0 ) {
952 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
953 LeftAQ = common.AQ[ CandAQ[SampledCase] ];
954 LeftQ = common.Q[ CandQ[SampledCase] ];
955
956 // Set the string properties
957 projectile->SplitUp();
958 projectile->SetFirstParton( LeftQ );
959 projectile->SetSecondParton( LeftAQ );
960 projectile->SetStatus( 0 );
961 G4int aAQ = std::abs( LeftAQ ), aQ = std::abs( LeftQ );
962 G4int NewCode = 0;
963 G4double aKsi = G4UniformRand();
964 // The string can have the quantum number of either a scalar or a vector (whose last digit
965 // of the PDG code is, respectively, 1 and 3). For simplicity only scalar is considered here.
966 if ( aAQ == aQ ) {
967 if ( aAQ != 3 ) {
968 NewCode = 111; // Pi0-meson
969 if ( aKsi < 0.5 ) {
970 NewCode = 221; // Eta -meson
971 if ( aKsi < 0.25 ) {
972 NewCode = 331; // Eta'-meson
973 }
974 }
975 } else {
976 NewCode = 221; // Eta -meson
977 if ( aKsi < 0.5 ) {
978 NewCode = 331; // Eta'-meson
979 }
980 }
981 } else {
982 if ( aAQ > aQ ) {
983 NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ;
984 } else {
985 NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ;
986 }
987 }
988
989 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
990 if ( ! TestParticle ) return false;
991 projectile->SetDefinition( TestParticle );
992 theParameters->SetProjMinDiffMass( 0.5 );
993 theParameters->SetProjMinNonDiffMass( 0.5 );
994
995 target->SetStatus( 4 ); // The target nucleon has annihilated 3->4
996 common.Pprojectile.setPx( 0.0 );
997 common.Pprojectile.setPy( 0.0 );
998 common.Pprojectile.setPz( 0.0 );
999 common.Pprojectile.setE( common.SqrtS );
1000 common.Pprojectile.transform( common.toLab );
1001
1002 G4LorentzVector Pquark = G4LorentzVector( 0.0, 0.0, common.SqrtS/2.0, common.SqrtS/2.0 );
1003 G4LorentzVector Paquark = G4LorentzVector( 0.0, 0.0, -common.SqrtS/2.0, common.SqrtS/2.0 );
1004 if ( common.RotateStrings ) {
1005 Pquark *= common.RandomRotation; Paquark *= common.RandomRotation;
1006 }
1007 Pquark.transform(common.toLab); projectile->GetNextParton()->Set4Momentum(Pquark);
1008 Paquark.transform(common.toLab); projectile->GetNextAntiParton()->Set4Momentum(Paquark);
1009
1010 // Calculation of the creation time
1011 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
1012 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
1013 projectile->SetPosition( target->GetPosition() );
1014 projectile->Set4Momentum( common.Pprojectile );
1015 projectile->IncrementCollisionCount( 1 );
1016 target->IncrementCollisionCount( 1 );
1017
1018 return true;
1019 } // End of if ( CandidatsN != 0 )
1020
1021 return true;
1022}
1023
1024
1025//============================================================================
1026
1027G4double G4FTFAnnihilation::ChooseX( G4double /* Alpha */, G4double /* Beta */ ) const {
1028 // If for sampling Xs other values of Alfa and Beta instead of 0.5 will be
1029 // chosen the method will be implemented
1030 //G4double tmp = Alpha*Beta;
1031 //tmp *= 1.0;
1032 return 0.5;
1033}
1034
1035
1036//============================================================================
1037
1038G4ThreeVector G4FTFAnnihilation::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1039 // @@ this method is used in FTFModel as well. Should go somewhere common!
1040 G4double Pt2 = 0.0;
1041 if ( AveragePt2 <= 0.0 ) {
1042 Pt2 = 0.0;
1043 } else {
1044 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1045 ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
1046 }
1047 G4double Pt = std::sqrt( Pt2 );
1048 G4double phi = G4UniformRand() * twopi;
1049 return G4ThreeVector ( Pt*std::cos( phi ), Pt*std::sin( phi ), 0.0 );
1050}
1051
1052
1053//============================================================================
1054
1055void G4FTFAnnihilation::UnpackBaryon( G4int IdPDG, G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1056 G4int AbsId = std::abs( IdPDG );
1057 Q1 = AbsId / 1000;
1058 Q2 = ( AbsId % 1000 ) / 100;
1059 Q3 = ( AbsId % 100 ) / 10;
1060 if ( IdPDG < 0 ) { Q1 = -Q1; Q2 = -Q2; Q3 = -Q3; } // Anti-baryon
1061 return;
1062}
1063
1064
1065//============================================================================
1066
1067G4FTFAnnihilation::G4FTFAnnihilation( const G4FTFAnnihilation& ) {
1068 throw G4HadronicException( __FILE__, __LINE__,
1069 "G4FTFAnnihilation copy constructor not meant to be called" );
1070}
1071
1072
1073//============================================================================
1074
1075const G4FTFAnnihilation & G4FTFAnnihilation::operator=( const G4FTFAnnihilation& ) {
1076 throw G4HadronicException( __FILE__, __LINE__,
1077 "G4FTFAnnihilation = operator not meant to be called" );
1078}
1079
1080
1081//============================================================================
1082
1083G4bool G4FTFAnnihilation::operator==( const G4FTFAnnihilation& ) const {
1084 throw G4HadronicException( __FILE__, __LINE__,
1085 "G4FTFAnnihilation == operator not meant to be called" );
1086}
1087
1088
1089//============================================================================
1090
1091G4bool G4FTFAnnihilation::operator!=( const G4FTFAnnihilation& ) const {
1092 throw G4HadronicException( __FILE__, __LINE__,
1093 "G4DiffractiveExcitation != operator not meant to be called" );
1094}
1095
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:226
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:34
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4double
double G4double
Definition: G4Types.hh:83
G4long
long G4long
Definition: G4Types.hh:87
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::getZ
double getZ() const
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::setX
void setX(double)
CLHEP::HepLorentzRotation::rotateY
HepLorentzRotation & rotateY(double delta)
Definition: LorentzRotation.cc:158
CLHEP::HepLorentzRotation::rotateZ
HepLorentzRotation & rotateZ(double delta)
Definition: LorentzRotation.cc:170
CLHEP::HepLorentzRotation::inverse
HepLorentzRotation inverse() const
CLHEP::HepLorentzVector::theta
double theta() const
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:173
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
G4FTFAnnihilation::~G4FTFAnnihilation
virtual ~G4FTFAnnihilation()
Definition: G4FTFAnnihilation.cc:81
G4FTFAnnihilation::Annihilate
virtual G4bool Annihilate(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4VSplitableHadron *&AdditionalString, G4FTFParameters *theParameters) const
Definition: G4FTFAnnihilation.cc:86
G4FTFParameters::SetTarMinNonDiffMass
void SetTarMinNonDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:576
G4FTFParameters::SetProjMinDiffMass
void SetProjMinDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:564
G4FTFParameters::SetTarMinDiffMass
void SetTarMinDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:572
G4FTFParameters::SetProjMinNonDiffMass
void SetProjMinNonDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:568
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:586
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
G4Parton::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
G4VSplitableHadron::SetTimeOfCreation
void SetTimeOfCreation(G4double aTime)
Definition: G4VSplitableHadron.hh:147
G4VSplitableHadron::SetStatus
void SetStatus(const G4int aStatus)
Definition: G4VSplitableHadron.hh:167
G4VSplitableHadron::SplitUp
virtual void SplitUp()=0
Definition: G4VSplitableHadron.cc:112
G4VSplitableHadron::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4VSplitableHadron.hh:137
G4VSplitableHadron::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4VSplitableHadron.hh:122
G4VSplitableHadron::GetNextParton
virtual G4Parton * GetNextParton()=0
G4VSplitableHadron::SetSecondParton
virtual void SetSecondParton(G4int PDGcode)=0
G4VSplitableHadron::GetNextAntiParton
virtual G4Parton * GetNextAntiParton()=0
G4VSplitableHadron::SetDefinition
void SetDefinition(const G4ParticleDefinition *aDefinition)
Definition: G4VSplitableHadron.hh:132
G4VSplitableHadron::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4VSplitableHadron.hh:127
G4VSplitableHadron::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4VSplitableHadron.hh:162
G4VSplitableHadron::IncrementCollisionCount
void IncrementCollisionCount(G4int aCount)
Definition: G4VSplitableHadron.hh:142
G4VSplitableHadron::SetFirstParton
virtual void SetFirstParton(G4int PDGcode)=0
G4VSplitableHadron::SetPosition
void SetPosition(const G4ThreeVector &aPosition)
Definition: G4VSplitableHadron.hh:157
sqr
T sqr(const T &x)
Definition: templates.hh:128