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