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