ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FTFAnnihilation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4FTFAnnihilation.cc
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 
76 G4FTFAnnihilation::G4FTFAnnihilation() {}
77 
78 
79 //============================================================================
80 
81 G4FTFAnnihilation::~G4FTFAnnihilation() {}
82 
83 
84 //============================================================================
85 
86 G4bool 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 ) {
186  X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // Pbar P
187  } else if ( ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114 ) {
188  X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // NeutrBar P
189  } else if ( ProjectilePDGcode == -3122 || ProjectilePDGcode == -3124 ) {
190  X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar P
191  } else if ( ProjectilePDGcode == -3112 || ProjectilePDGcode == -3114 ) {
192  X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma-Bar P
193  } else if ( ProjectilePDGcode == -3212 || ProjectilePDGcode == -3214 ) {
194  X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar P
195  } else if ( ProjectilePDGcode == -3222 || ProjectilePDGcode == -3224 ) {
196  X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma+Bar P
197  } else if ( ProjectilePDGcode == -3312 || ProjectilePDGcode == -3314 ) {
198  X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi-Bar P
199  } else if ( ProjectilePDGcode == -3322 || ProjectilePDGcode == -3324 ) {
200  X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi0Bar P
201  } else if ( ProjectilePDGcode == -3334 ) {
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 ) {
208  X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // Pbar N
209  } else if ( ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114 ) {
210  X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // NeutrBar N
211  } else if ( ProjectilePDGcode == -3122 || ProjectilePDGcode == -3124 ) {
212  X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar N
213  } else if ( ProjectilePDGcode == -3112 || ProjectilePDGcode == -3114 ) {
214  X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma-Bar N
215  } else if ( ProjectilePDGcode == -3212 || ProjectilePDGcode == -3214 ) {
216  X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar N
217  } else if ( ProjectilePDGcode == -3222 || ProjectilePDGcode == -3224 ) {
218  X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma+Bar N
219  } else if ( ProjectilePDGcode == -3312 || ProjectilePDGcode == -3314 ) {
220  X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi-Bar N
221  } else if ( ProjectilePDGcode == -3322 || ProjectilePDGcode == -3324 ) {
222  X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi0Bar N
223  } else if ( ProjectilePDGcode == -3334 ) {
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 
281 G4bool G4FTFAnnihilation::
282 Create3QuarkAntiQuarkStrings( 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 
552 G4int G4FTFAnnihilation::
553 Create1DiquarkAntiDiquarkString( 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 
645 G4int G4FTFAnnihilation::
646 Create2QuarkAntiQuarkStrings( 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 
905 G4bool G4FTFAnnihilation::
906 Create1QuarkAntiQuarkString( 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 
1027 G4double 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 
1038 G4ThreeVector 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 
1055 void 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 
1067 G4FTFAnnihilation::G4FTFAnnihilation( const G4FTFAnnihilation& ) {
1068  throw G4HadronicException( __FILE__, __LINE__,
1069  "G4FTFAnnihilation copy constructor not meant to be called" );
1070 }
1071 
1072 
1073 //============================================================================
1074 
1075 const G4FTFAnnihilation & G4FTFAnnihilation::operator=( const G4FTFAnnihilation& ) {
1076  throw G4HadronicException( __FILE__, __LINE__,
1077  "G4FTFAnnihilation = operator not meant to be called" );
1078 }
1079 
1080 
1081 //============================================================================
1082 
1083 G4bool G4FTFAnnihilation::operator==( const G4FTFAnnihilation& ) const {
1084  throw G4HadronicException( __FILE__, __LINE__,
1085  "G4FTFAnnihilation == operator not meant to be called" );
1086 }
1087 
1088 
1089 //============================================================================
1090 
1091 G4bool G4FTFAnnihilation::operator!=( const G4FTFAnnihilation& ) const {
1092  throw G4HadronicException( __FILE__, __LINE__,
1093  "G4DiffractiveExcitation != operator not meant to be called" );
1094 }
1095