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