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