Geant4  10.02.p02
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 97625 2016-06-06 13:35:58Z 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 #include "G4ParticleTable.hh" // Uzhi March 2016
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 //theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016 ??? for other Anti_bar annih.
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( 3 ); // 2->3 Uzhi Oct 2014
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 
148  G4double MesonProdThreshold = projectile->GetDefinition()->GetPDGMass() +
149  target->GetDefinition()->GetPDGMass() +
150  ( 2.0*140.0 + 16.0 )*MeV; // 2 Mpi + DeltaE
151  G4double Prel2 = S*S + M0projectile2*M0projectile2 + M0target2*M0target2 -
152  2.0*S*M0projectile2 - 2.0*S*M0target2 - 2.0*M0projectile2*M0target2;
153  Prel2 /= S;
154  //G4cout << "Prel2 " << Prel2 << G4endl;
155  if ( Prel2 <= 0.0 ) { // *MeV*MeV 1600.
156  // Annihilation at rest! Values are copied from Parameters
157  X_a = 625.1; // mb // 3-shirt diagram
158  X_b = 0.0; // 9.780 12 Dec. 2012; // mb // anti-quark-quark annihilation
159  X_c = 49.989; // mb
160  X_d = 6.614; // mb
161 
162  #ifdef debugFTFannih
163  G4cout << "Annih at Rest X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d
164  << G4endl;
165  #endif
166 
167  } else { // Annihilation in flight!
168  G4double FlowF = 1.0 / std::sqrt( Prel2 )*GeV;
169 
170  // Process cross sections
171  X_a = 25.0*FlowF; // mb 3-shirt diagram
172  if ( SqrtS < MesonProdThreshold ) {
173  X_b = 3.13 + 140.0*G4Pow::GetInstance()->powA( ( MesonProdThreshold - SqrtS )/GeV, 2.5 );
174  } else {
175  X_b = 6.8*GeV / SqrtS; // mb anti-quark-quark annihilation
176  }
177  if ( projectile->GetDefinition()->GetPDGMass() + target->GetDefinition()->GetPDGMass()
178  > SqrtS ) {
179  X_b = 0.0;
180  }
181  // This can be in an interaction of low energy anti-baryon with off-shell nuclear nucleon
182  X_c = 2.0 * FlowF * sqr( projectile->GetDefinition()->GetPDGMass() +
183  target->GetDefinition()->GetPDGMass() ) / S; // mb re-arrangement of
184  // 2 quarks and 2 anti-quarks
185  X_d = 23.3*GeV*GeV / S; // mb anti-quark-quark string creation
186 
187  #ifdef debugFTFannih
188  G4cout << "Annih in Flight X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d
189  << G4endl << "SqrtS MesonProdThreshold " << SqrtS << " " << MesonProdThreshold
190  << G4endl;
191  #endif
192 
193  }
194 
195  if ((ProjectilePDGcode == -2212 || ProjectilePDGcode == -2214)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
196  X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // Pbar P
197  } else if ((ProjectilePDGcode == -2212 || ProjectilePDGcode == -2214)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
198  X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // Pbar N
199  } else if ((ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
200  X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // NeutrBar P
201  } else if ((ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
202  X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // NeutrBar N
203  } else if ((ProjectilePDGcode == -3122 || ProjectilePDGcode == -3124)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
204  X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar P
205  } else if ((ProjectilePDGcode == -3122 || ProjectilePDGcode == -3124)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
206  X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar N
207  } else if ((ProjectilePDGcode == -3112 || ProjectilePDGcode == -3114)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
208  X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma-Bar P
209  } else if ((ProjectilePDGcode == -3112 || ProjectilePDGcode == -3114)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
210  X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma-Bar N
211  } else if ((ProjectilePDGcode == -3212 || ProjectilePDGcode == -3214)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
212  X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar P
213  } else if ((ProjectilePDGcode == -3212 || ProjectilePDGcode == -3214)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
214  X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar N
215  } else if ((ProjectilePDGcode == -3222 || ProjectilePDGcode == -3224)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
216  X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma+Bar P
217  } else if ((ProjectilePDGcode == -3222 || ProjectilePDGcode == -3224)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
218  X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma+Bar N
219  } else if ((ProjectilePDGcode == -3312 || ProjectilePDGcode == -3314)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
220  X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi-Bar P
221  } else if ((ProjectilePDGcode == -3312 || ProjectilePDGcode == -3314)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
222  X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi-Bar N
223  } else if ((ProjectilePDGcode == -3322 || ProjectilePDGcode == -3324)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
224  X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi0Bar P
225  } else if ((ProjectilePDGcode == -3322 || ProjectilePDGcode == -3324)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
226  X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi0Bar N
227  } else if ( ProjectilePDGcode == -3334 && ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) {
228  X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // Omega-Bar P
229  } else if ( ProjectilePDGcode == -3334 && ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) {
230  X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // Omega-Bar N
231  } else {
232  G4cout << "Unknown anti-baryon for FTF annihilation: PDGcodes - "
233  << ProjectilePDGcode << " " << TargetPDGcode << G4endl;
234  }
235 
236  #ifdef debugFTFannih
237  G4cout << "Annih Actual X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d << G4endl;
238  #endif
239 
240  G4double Xannihilation = X_a + X_b + X_c + X_d;
241 //X_a=0.; // Uzhi
242 //X_b=0.;
243 //X_c=0.;
244 //X_d=0.;
245 //Xannihilation = X_a + X_b + X_c + X_d;
246 
247  // Projectile unpacking
248  G4int AQ[3];
249  UnpackBaryon( ProjectilePDGcode, AQ[0], AQ[1], AQ[2] );
250 
251  // Target unpacking
252  G4int Q[3];
253  UnpackBaryon( TargetPDGcode, Q[0], Q[1], Q[2] );
254 
255  G4double Ksi = G4UniformRand();
256 
257  if ( Ksi < X_a / Xannihilation ) {
258 
259  // Simulation of 3 anti-quark-quark strings creation
260  // Sampling of anti-quark order in projectile
261 
262  #ifdef debugFTFannih
263  G4cout << "Process a, 3 shirt diagram" << G4endl;
264  #endif
265 
266  G4int SampledCase = G4RandFlat::shootInt( G4long( 6 ) );
267 
268  G4int Tmp1( 0 ), Tmp2( 0 );
269  if ( SampledCase == 0 ) {
270  } else if ( SampledCase == 1 ) {
271  Tmp1 = AQ[1]; AQ[1] = AQ[2]; AQ[2] = Tmp1;
272  } else if ( SampledCase == 2 ) {
273  Tmp1 = AQ[0]; AQ[0] = AQ[1]; AQ[1] = Tmp1;
274  } else if ( SampledCase == 3 ) {
275  Tmp1 = AQ[0]; Tmp2 = AQ[1]; AQ[0] = AQ[2]; AQ[1] = Tmp1; AQ[2] = Tmp2;
276  } else if ( SampledCase == 4 ) {
277  Tmp1 = AQ[0]; Tmp2 = AQ[1]; AQ[0] = Tmp2; AQ[1] = AQ[2]; AQ[2] = Tmp1;
278  } else if ( SampledCase == 5 ) {
279  Tmp1 = AQ[0]; Tmp2 = AQ[1]; AQ[0] = AQ[2]; AQ[1] = Tmp2; AQ[2] = Tmp1;
280  }
281 
282  // Set the string properties
283 
284  //G4cout << "String 1 " << AQ[0] << " " << Q[0] << G4endl;
285  projectile->SplitUp();
286  projectile->SetFirstParton( AQ[0] );
287  projectile->SetSecondParton( Q[0] );
288  projectile->SetStatus( 0 );
289 
290 // Uzhi March 2016 start
291 G4int aAQ, aQ;
292 aAQ=std::abs( AQ[0] ); aQ=std::abs( Q[0] );
293 G4int NewCode;
294 G4double aKsi = G4UniformRand();
295 
296 if ( aAQ == aQ )
297 {
298  if ( aAQ != 3 )
299  {
300  NewCode = 111; // Pi0-meson
301  if ( aKsi < 0.5 )
302  {
303  NewCode = 221; // Eta -meson
304  if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson
305  }
306  } else
307  {
308  NewCode = 221; // Eta -meson
309  if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson
310  }
311 } else
312 {
313  if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/AQ[0]; }
314  else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/Q[0]; }
315 }
316 
317 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
318 if(!TestParticle) return false;
319 projectile->SetDefinition( TestParticle );
320 
321 theParameters->SetProjMinDiffMass( 0.5 ); // Uzhi 2016 M+140 ??
322 theParameters->SetProjMinNonDiffMass( 0.5 ); // Uzhi 2016 M+140 ??
323 // Uzhi March 2016 end
324 
325  //G4cout << "String 2 " << Q[1] << " " << AQ[1] << G4endl;
326  target->SplitUp();
327  target->SetFirstParton( Q[1] );
328  target->SetSecondParton( AQ[1] );
329  target->SetStatus( 0 );
330 
331 // Uzhi March 2016 Start
332 aAQ=std::abs( AQ[1] ); aQ=std::abs( Q[1] ); aKsi = G4UniformRand();
333 if ( aAQ == aQ )
334 {
335  if ( aAQ != 3 )
336  {
337  NewCode = 111; // Pi0-meson
338  if ( aKsi < 0.5 )
339  {
340  NewCode = 221; // Eta -meson
341  if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson
342  }
343  } else
344  {
345  NewCode = 221; // Eta -meson
346  if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson
347  }
348 } else
349 {
350  if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/AQ[1]; }
351  else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/Q[1]; }
352 }
353 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
354 if(!TestParticle) return false;
355 target->SetDefinition( TestParticle );
356 
357 theParameters->SetTarMinDiffMass( 0.5 ); // Uzhi 2016 M+140 ??
358 theParameters->SetTarMinNonDiffMass( 0.5 ); // Uzhi 2016 M+140 ??
359 // Uzhi March 2016 end
360 
361 
362  //G4cout << "String 3 " << AQ[2] << " " << Q[2] << G4endl;
363  AdditionalString = new G4DiffractiveSplitableHadron();
364 
365 // Uzhi March 2016 start
366 aAQ=std::abs( AQ[2] ); aQ=std::abs( Q[2] ); aKsi = G4UniformRand();
367 
368 if ( aAQ == aQ )
369 {
370  if ( aAQ != 3 )
371  {
372  NewCode = 111; // Pi0-meson
373  if ( aKsi < 0.5 )
374  {
375  NewCode = 221; // Eta -meson
376  if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson
377  }
378  } else
379  {
380  NewCode = 221; // Eta -meson
381  if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson
382  }
383 } else
384 {
385  if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/AQ[2]; }
386  else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/Q[2]; }
387 }
388 
389 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
390 if(!TestParticle) return false;
391 AdditionalString->SetDefinition( TestParticle );
392 // Uzhi March 2016 end
393 
394  AdditionalString->SplitUp();
395  AdditionalString->SetFirstParton( AQ[2] );
396  AdditionalString->SetSecondParton( Q[2] );
397  AdditionalString->SetStatus( 0 );
398  //G4cout << G4endl << "*AdditionalString in Annih" << AdditionalString << G4endl;
399 
400  // Sampling kinematical properties
401  // 1 string AQ[0]-Q[0]// 2 string AQ[1]-Q[1]// 3 string AQ[2]-Q[2]
402 
403  G4ThreeVector Quark_Mom[6];
404  G4double ModMom2[6]; //ModMom[6]
405 
406  AveragePt2 = 200.0*200.0; maxPtSquare = S;
407 
408  G4double SumMt( 0.0 );
409  G4double MassQ2 = 0.0; // 100.0*100.0*MeV*MeV;
410  G4int NumberOfTries( 0 );
411  G4double ScaleFactor( 1.0 );
412 
413  const G4int maxNumberOfLoops = 1000;
414  G4int loopCounter = 0;
415  do {
416  NumberOfTries++;
417  if ( NumberOfTries == 100*(NumberOfTries/100) ) {
418  // At large number of tries it would be better to reduce the values of <Pt^2>
419  ScaleFactor /= 2.0;
420  AveragePt2 *= ScaleFactor;
421  }
422  G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
423  for ( G4int i = 0; i < 6; i++ ) {
424  Quark_Mom [i] = GaussianPt( AveragePt2, maxPtSquare );
425  PtSum += Quark_Mom[i];
426  }
427  PtSum /= 6.0;
428  SumMt = 0.0;
429  for( G4int i = 0; i < 6; i++ ) {
430  Quark_Mom[i] -= PtSum;
431  //ModMom[i] = Quark_Mom[i].mag();
432  ModMom2[i] = Quark_Mom[i].mag2();
433  SumMt += std::sqrt( ModMom2[i] + MassQ2 );
434  }
435  } while ( ( SumMt > SqrtS ) &&
436  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
437  if ( loopCounter >= maxNumberOfLoops ) {
438  return false;
439  }
440 
441  G4double WminusTarget( 0.0 ), WplusProjectile( 0.0 );
442 
443  // Closed is variant with sampling of Xs at minimum
444  //G4double SumMod_anti = ModMom[0] + ModMom[1] + ModMom[2];
445  //Quark_Mom[0].setZ( ModMom[0]/SumMod_anti );
446  //Quark_Mom[1].setZ( ModMom[1]/SumMod_anti );
447  //Quark_Mom[2].setZ( ModMom[2]/SumMod_anti );
448  //G4double SumMod_bary = ModMom[3] + ModMom[4] + ModMom[5];
449  //Quark_Mom[3].setZ( ModMom[3]/SumMod_bary );
450  //Quark_Mom[4].setZ( ModMom[4]/SumMod_bary );
451  //Quark_Mom[5].setZ( ModMom[5]/SumMod_bary );
452  //G4double Alfa = SumMod_anti*SumMod_anti;
453  //G4double Beta = SumMod_bary*SumMod_bary;
454  //G4double DecayMomentum2 = S*S + Alfa*Alfa + Beta*Beta
455  // - 2.0*S*Alfa - 2.0*S*Beta - 2.0*Alfa*Beta;
456  //WminusTarget = ( S - Alfa + Beta + std::sqrt( DecayMomentum2 ) )/2.0/SqrtS;
457  //WplusProjectile = SqrtS - Beta/WminusTarget;
458  // Closed is variant with sampling of Xs at minimum
459 
460  // Sampling X's of anti-baryon
461  G4double Alfa_R = 0.5;
462  NumberOfTries = 0;
463  ScaleFactor = 1.0;
464  G4bool Succes( true );
465 
466  loopCounter = 0;
467  do {
468 
469  Succes = true;
470  NumberOfTries++;
471  if ( NumberOfTries == 100*(NumberOfTries/100) ) {
472  // At large number of tries it would be better to reduce the values of Pt's
473  ScaleFactor /= 2.0;
474  }
475 
476  if ( Alfa_R == 1.0 ) {
477  G4double Xaq1 = 1.0 - std::sqrt( G4UniformRand() );
478  G4double Xaq2 = (1.0 - Xaq1) * G4UniformRand();
479  G4double Xaq3 = 1.0 - Xaq1 - Xaq2;
480  Quark_Mom[0].setZ( Xaq1 ); Quark_Mom[1].setZ( Xaq2 ); Quark_Mom[2].setZ( Xaq3 );
481  } else {
482  G4double Xaq1 = sqr( G4UniformRand() );
483  G4double Xaq2 = (1.0 - Xaq1)*sqr( std::sin( pi/2.0*G4UniformRand() ) );
484  G4double Xaq3 = 1.0 - Xaq1 - Xaq2;
485  Quark_Mom[0].setZ( Xaq1 ); Quark_Mom[1].setZ( Xaq2 ); Quark_Mom[2].setZ( Xaq3 );
486  }
487 
488  // Sampling X's of baryon
489  if ( Alfa_R == 1.0 ) {
490  G4double Xq1 = 1.0 - std::sqrt( G4UniformRand() );
491  G4double Xq2 = (1.0 - Xq1) * G4UniformRand();
492  G4double Xq3 = 1.0 - Xq1 - Xq2;
493  Quark_Mom[3].setZ( Xq1 ); Quark_Mom[4].setZ( Xq2 ); Quark_Mom[5].setZ( Xq3 );
494  } else {
495  G4double Xq1 = sqr( G4UniformRand() );
496  G4double Xq2 = (1.0 - Xq1) * sqr( std::sin( pi/2.0*G4UniformRand() ) );
497  G4double Xq3 = 1.0 - Xq1 - Xq2;
498  Quark_Mom[3].setZ( Xq1 ); Quark_Mom[4].setZ( Xq2 ); Quark_Mom[5].setZ( Xq3 );
499  }
500 
501  G4double Alfa( 0.0 ), Beta( 0.0 );
502  for ( G4int i = 0; i < 3; i++ ) { // For Anti-baryon
503  if ( Quark_Mom[i].getZ() != 0.0 ) {
504  Alfa += ( ScaleFactor * ModMom2[i] + MassQ2 ) / Quark_Mom[i].getZ();
505  } else {
506  Succes = false;
507  }
508  }
509  for ( G4int i = 3; i < 6; i++ ) { // For baryon
510  if ( Quark_Mom[i].getZ() != 0.0 ) {
511  Beta += ( ScaleFactor * ModMom2[i] + MassQ2 ) / Quark_Mom[i].getZ();
512  } else {
513  Succes = false;
514  }
515  }
516 
517  if ( ! Succes ) continue;
518 
519  if ( std::sqrt( Alfa ) + std::sqrt( Beta ) > SqrtS ) {
520  Succes = false;
521  continue;
522  }
523 
524  G4double DecayMomentum2 = S*S + Alfa*Alfa + Beta*Beta
525  - 2.0*S*Alfa - 2.0*S*Beta - 2.0*Alfa*Beta;
526 
527  WminusTarget = ( S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS;
528  WplusProjectile = SqrtS - Beta/WminusTarget;
529 
530  } while ( ( ! Succes ) &&
531  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
532  if ( loopCounter >= maxNumberOfLoops ) {
533  return false;
534  }
535 
536  G4double SqrtScaleF = std::sqrt( ScaleFactor );
537  for ( G4int i = 0; i < 3; i++ ) {
538  G4double Pz = WplusProjectile * Quark_Mom[i].getZ() / 2.0 -
539  ( ScaleFactor * ModMom2[i] + MassQ2 ) /
540  ( 2.0 * WplusProjectile * Quark_Mom[i].getZ() );
541  Quark_Mom[i].setZ( Pz );
542  if ( ScaleFactor != 1.0 ) {
543  Quark_Mom[i].setX( SqrtScaleF * Quark_Mom[i].getX() );
544  Quark_Mom[i].setY( SqrtScaleF * Quark_Mom[i].getY() );
545  }
546  }
547  for ( G4int i = 3; i < 6; i++ ) {
548  G4double Pz = -WminusTarget * Quark_Mom[i].getZ() / 2.0 +
549  ( ScaleFactor * ModMom2[i] + MassQ2 ) /
550  ( 2.0 * WminusTarget * Quark_Mom[i].getZ() );
551  Quark_Mom[i].setZ( Pz );
552  if ( ScaleFactor != 1.0 ) {
553  Quark_Mom[i].setX( SqrtScaleF * Quark_Mom[i].getX() );
554  Quark_Mom[i].setY( SqrtScaleF * Quark_Mom[i].getY() );
555  }
556  }
557  //G4cout << "Sum AQ " << Quark_Mom[0] + Quark_Mom[1] + Quark_Mom[2] << G4endl
558  // << "Sum Q " << Quark_Mom[3] + Quark_Mom[4] + Quark_Mom[5] << G4endl;
559 
560  G4ThreeVector tmp = Quark_Mom[0] + Quark_Mom[3];
561  G4LorentzVector Pstring1( tmp, std::sqrt( Quark_Mom[0].mag2() + MassQ2 ) +
562  std::sqrt( Quark_Mom[3].mag2() + MassQ2 ) );
563  G4double Ystring1 = Pstring1.rapidity();
564 
565  //G4cout << "Mom 1 string " << G4endl << Quark_Mom[0] << G4endl << Quark_Mom[3] << G4endl
566  // << tmp << " " << tmp.mag() << G4endl;
567  //G4cout << "1 str " << Pstring1 << " " << Pstring1.mag() << " " << Ystring1 << G4endl;
568 
569  tmp = Quark_Mom[1] + Quark_Mom[4];
570  G4LorentzVector Pstring2( tmp, std::sqrt( Quark_Mom[1].mag2() + MassQ2 ) +
571  std::sqrt( Quark_Mom[4].mag2() + MassQ2 ) );
572  G4double Ystring2 = Pstring2.rapidity();
573 
574  //G4cout << "Mom 2 string " << G4endl << Quark_Mom[1] << G4endl << Quark_Mom[4] << G4endl
575  // << tmp << " " << tmp.mag() << G4endl;
576  //G4cout << "2 str " << Pstring2 << " " << Pstring2.mag() << " " << Ystring2 << G4endl;
577 
578  tmp = Quark_Mom[2] + Quark_Mom[5];
579  G4LorentzVector Pstring3( tmp, std::sqrt( Quark_Mom[2].mag2() + MassQ2 ) +
580  std::sqrt( Quark_Mom[5].mag2() + MassQ2 ) );
581  G4double Ystring3 = Pstring3.rapidity();
582 
583  //G4cout << "Mom 3 string " << G4endl << Quark_Mom[2] << G4endl << Quark_Mom[5] << G4endl
584  // << tmp << " " << tmp.mag() << G4endl;
585  //G4cout << "3 str " << Pstring3 << " " << Pstring3.mag() << " " << Ystring3 << G4endl
586  // << "SumE " << Pstring1.e() + Pstring2.e() + Pstring3.e() << G4endl
587  // << Pstring1.mag() << " " <<Pstring2.mag() << " " << Pstring3.mag() << G4endl;
588  //G4int Uzhi; G4cin >> Uzhi;
589 
590  G4LorentzVector LeftString( 0.0, 0.0, 0.0, 0.0 );
591  if ( Ystring1 > Ystring2 && Ystring2 > Ystring3 ) {
592  Pprojectile = Pstring1;
593  LeftString = Pstring2;
594  Ptarget = Pstring3;
595  }
596  if ( Ystring1 > Ystring3 && Ystring3 > Ystring2 ) {
597  Pprojectile = Pstring1;
598  LeftString = Pstring3;
599  Ptarget = Pstring2;
600  }
601 
602  if ( Ystring2 > Ystring1 && Ystring1 > Ystring3 ) {
603  Pprojectile = Pstring2;
604  LeftString = Pstring1;
605  Ptarget = Pstring3;
606  }
607  if ( Ystring2 > Ystring3 && Ystring3 > Ystring1 ) {
608  Pprojectile = Pstring2;
609  LeftString = Pstring3;
610  Ptarget = Pstring1;
611  }
612 
613  if ( Ystring3 > Ystring1 && Ystring1 > Ystring2 ) {
614  Pprojectile = Pstring3;
615  LeftString = Pstring1;
616  Ptarget = Pstring2;
617  }
618  if ( Ystring3 > Ystring2 && Ystring2 > Ystring1 ) {
619  Pprojectile = Pstring3;
620  LeftString = Pstring2;
621  Ptarget = Pstring1;
622  }
623  //G4cout << "SumP " << Pprojectile + LeftString + Ptarget << " " << SqrtS << G4endl;
624 
625  Pprojectile.transform( toLab );
626  LeftString.transform( toLab );
627  Ptarget.transform( toLab );
628  //G4cout << "SumP " << Pprojectile + LeftString + Ptarget << " " << SqrtS << G4endl;
629 
630  // Calculation of the creation time
631  projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
632  projectile->SetPosition( target->GetPosition() );
633  AdditionalString->SetTimeOfCreation( target->GetTimeOfCreation() );
634  AdditionalString->SetPosition( target->GetPosition() );
635  // Creation time and position of target nucleon were determined in
636  // ReggeonCascade() of G4FTFModel
637 
638  //G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl;
639  projectile->Set4Momentum( Pprojectile );
640  AdditionalString->Set4Momentum( LeftString );
641  target->Set4Momentum( Ptarget );
642  projectile->IncrementCollisionCount( 1 );
643  AdditionalString->IncrementCollisionCount( 1 );
644  target->IncrementCollisionCount( 1 );
645 
646 theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016
647 
648  return true;
649 
650  } // End of if ( Ksi < X_a / Xannihilation )
651 
652  // Simulation of anti-diquark-diquark string creation
653 
654  if ( Ksi < (X_a + X_b) / Xannihilation ) {
655 
656  #ifdef debugFTFannih
657  G4cout << "Process b, quark - anti-quark annihilation, di-q - anti-di-q string" << G4endl;
658  #endif
659 
660  G4int CandidatsN( 0 ), CandAQ[9][2], CandQ[9][2];
661  G4int LeftAQ1( 0 ), LeftAQ2( 0 ), LeftQ1( 0 ), LeftQ2( 0 );
662 
663  for ( G4int iAQ = 0; iAQ < 3; iAQ++ ) {
664  for ( G4int iQ = 0; iQ < 3; iQ++ ) {
665  if ( -AQ[iAQ] == Q[iQ] ) {
666  if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; }
667  if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; }
668  if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; }
669  if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; }
670  if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; }
671  if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; }
672  CandidatsN++;
673  }
674  }
675  }
676  //G4cout << "CandidatsN " << CandidatsN << G4endl;
677 
678  if ( CandidatsN != 0 ) {
679  G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
680  LeftAQ1 = AQ[ CandAQ[SampledCase][0] ];
681  LeftAQ2 = AQ[ CandAQ[SampledCase][1] ];
682  LeftQ1 = Q[ CandQ[SampledCase][0] ];
683  LeftQ2 = Q[ CandQ[SampledCase][1] ];
684 
685  // Build anti-diquark and diquark
686  G4int Anti_DQ( 0 ), DQ( 0 );
687  if ( std::abs( LeftAQ1 ) > std::abs( LeftAQ2 ) ) {
688  Anti_DQ = 1000*LeftAQ1 + 100*LeftAQ2 - 3; // 1
689  } else {
690  Anti_DQ = 1000*LeftAQ2 + 100*LeftAQ1 - 3; // 1
691  }
692  //if ( G4UniformRand() > 0.5 ) Anti_DQ -= 2;
693  if ( std::abs( LeftQ1 ) > std::abs( LeftQ2 ) ) {
694  DQ = 1000*LeftQ1 + 100*LeftQ2 + 3; // 1
695  } else {
696  DQ = 1000*LeftQ2 + 100*LeftQ1 + 3; // 1
697  }
698  // if ( G4UniformRand() > 0.5 ) DQ += 2;
699 
700  // Set the string properties
701  //G4cout << "Left ADiQ DiQ " << Anti_DQ << " " << DQ << G4endl;
702  projectile->SplitUp();
703  //projectile->SetFirstParton( Anti_DQ );
704  //projectile->SetSecondParton( DQ );
705  projectile->SetFirstParton( DQ );
706  projectile->SetSecondParton( Anti_DQ );
707  projectile->SetStatus( 0 );
708  target->SetStatus( 4 ); // The target nucleon has annihilated 3->4 Uzhi Oct 2014
709  Pprojectile.setPx( 0.0 );
710  Pprojectile.setPy( 0.0 );
711  Pprojectile.setPz( 0.0 );
712  Pprojectile.setE( SqrtS );
713  Pprojectile.transform( toLab );
714 // Uzhi March 2016 if QQ_QQbar will interact Set Mmin, MdifMin
715 
716  // Calculation of the creation time
717  projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
718  projectile->SetPosition( target->GetPosition() );
719  // Creation time and position of target nucleon were determined in
720  // ReggeonCascade() of G4FTFModel
721 
722  //G4cout << "Mproj " << Pprojectile.mag() << G4endl
723  // << "Mtarg " << Ptarget.mag() << G4endl;
724  projectile->Set4Momentum( Pprojectile );
725 
726  projectile->IncrementCollisionCount( 1 );
727  target->IncrementCollisionCount( 1 );
728 
729 //theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016
730 // In the case baryon and anti-baryon are created. Thus the antibaryon can annihilate latter.
731 
732  return true;
733  }
734 
735  } // End of if ( Ksi < (X_a + X_b) / Xannihilation )
736 
737  if ( Ksi < ( X_a + X_b + X_c ) / Xannihilation ) {
738 
739  // Simulation of 2 anti-quark-quark strings creation
740 
741  #ifdef debugFTFannih
742  G4cout << "Process c, quark - anti-quark and string junctions annihilation, 2 strings left."
743  << G4endl;
744  #endif
745 
746  G4int CandidatsN( 0 ), CandAQ[9][2], CandQ[9][2];
747  G4int LeftAQ1( 0 ), LeftAQ2( 0 ), LeftQ1( 0 ), LeftQ2( 0 );
748 
749  for ( G4int iAQ = 0; iAQ < 3; iAQ++ ) {
750  for ( G4int iQ = 0; iQ < 3; iQ++ ) {
751  if ( -AQ[iAQ] == Q[iQ] ) {
752  if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; }
753  if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; }
754  if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; }
755  if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; }
756  if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; }
757  if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; }
758  CandidatsN++;
759  }
760  }
761  }
762  //G4cout << "CandidatsN " << CandidatsN << G4endl;
763 
764  if ( CandidatsN != 0 ) {
765  G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
766  LeftAQ1 = AQ[ CandAQ[SampledCase][0] ];
767  LeftAQ2 = AQ[ CandAQ[SampledCase][1] ];
768  if ( G4UniformRand() < 0.5 ) {
769  LeftQ1 = Q[ CandQ[SampledCase][0] ];
770  LeftQ2 = Q[ CandQ[SampledCase][1] ];
771  } else {
772  LeftQ2 = Q[ CandQ[SampledCase][0] ];
773  LeftQ1 = Q[ CandQ[SampledCase][1] ];
774  }
775 
776  // Set the string properties
777  //G4cout << "String 1 " << LeftAQ1 << " " << LeftQ1 << G4endl;
778  projectile->SplitUp();
779  projectile->SetFirstParton( LeftAQ1 );
780  projectile->SetSecondParton( LeftQ1 );
781  projectile->SetStatus( 0 );
782 
783 // Uzhi March 2016 start
784 G4int aAQ, aQ;
785 aAQ=std::abs( LeftAQ1 ); aQ=std::abs( LeftQ1 );
786 
787 G4int NewCode;
788 G4double aKsi = G4UniformRand();
789 
790 if ( aAQ == aQ )
791 {
792  if ( aAQ != 3 )
793  {
794  NewCode = 111; // Pi0-meson
795  if ( aKsi < 0.5 )
796  {
797  NewCode = 221; // Eta -meson
798  if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson
799  }
800  } else
801  {
802  NewCode = 221; // Eta -meson
803  if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson
804  }
805 } else
806 {
807  if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ1; }
808  else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ1; }
809 }
810 
811 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
812 if(!TestParticle) return false;
813 projectile->SetDefinition( TestParticle );
814 theParameters->SetProjMinDiffMass( 0.5 ); // (0.5) // GeV Uzhi March 2016 ???
815 theParameters->SetProjMinNonDiffMass( 0.5 );
816 // Uzhi March 2016 end
817 
818  //G4cout << "String 2 " << LeftAQ2 << " " << LeftQ2 << G4endl;
819  target->SplitUp();
820  target->SetFirstParton( LeftQ2 );
821  target->SetSecondParton( LeftAQ2 );
822  target->SetStatus( 0 );
823 
824 // Uzhi March 2016 start
825 aAQ=std::abs( LeftAQ2 ); aQ=std::abs( LeftQ2 ); aKsi = G4UniformRand();
826 
827 if ( aAQ == aQ )
828 {
829  if ( aAQ != 3 )
830  {
831  NewCode = 111; // Pi0-meson
832  if ( aKsi < 0.5 )
833  {
834  NewCode = 221; // Eta -meson
835  if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson
836  }
837  } else
838  {
839  NewCode = 221; // Eta -meson
840  if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson
841  }
842 } else
843 {
844  if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ2; }
845  else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ2; }
846 }
847 
848 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
849 if(!TestParticle) return false;
850 target->SetDefinition( TestParticle );
851 theParameters->SetTarMinDiffMass( 0.5 ); // Uzhi March 2016 ???
852 theParameters->SetTarMinNonDiffMass( 0.5 );
853 // Uzhi March 2016
854 
855  // Sampling kinematical properties
856  // 1 string LeftAQ1-LeftQ1// 2 string LeftAQ2-LeftQ2
857  G4ThreeVector Quark_Mom[4];
858  G4double ModMom2[4]; //ModMom[4],
859 
860  AveragePt2 = 200.0*200.0; maxPtSquare = S;
861 
862  G4double SumMt( 0.0 );
863  G4double MassQ2 = 0.0; //100.0*100.0*MeV*MeV;
864  G4int NumberOfTries( 0 );
865  G4double ScaleFactor( 1.0 );
866 
867  const G4int maxNumberOfLoops = 1000;
868  G4int loopCounter = 0;
869  do {
870  NumberOfTries++;
871  if ( NumberOfTries == 100*(NumberOfTries/100) ) {
872  // At large number of tries it would be better to reduce the values of <Pt^2>
873  ScaleFactor /= 2.0;
874  AveragePt2 *= ScaleFactor;
875  }
876  G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
877  for( G4int i = 0; i < 4; i++ ) {
878  Quark_Mom[i] = GaussianPt( AveragePt2, maxPtSquare );
879  PtSum += Quark_Mom[i];
880  }
881  PtSum /= 4.0;
882  SumMt = 0.0;
883  for ( G4int i = 0; i < 4; i++ ) {
884  Quark_Mom[i] -= PtSum;
885  //ModMom[i] = Quark_Mom[i].mag();
886  ModMom2[i] = Quark_Mom[i].mag2();
887  SumMt += std::sqrt( ModMom2[i] + MassQ2 );
888  }
889  } while ( ( SumMt > SqrtS ) &&
890  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
891  if ( loopCounter >= maxNumberOfLoops ) {
892  return false;
893  }
894 
895  G4double WminusTarget( 0.0 ), WplusProjectile( 0.0 );
896 
897  // Sampling X's of anti-baryon
898  G4double Alfa_R = 0.5;
899  NumberOfTries = 0;
900  ScaleFactor = 1.0;
901  G4bool Succes( true );
902 
903  loopCounter = 0;
904  do {
905 
906  Succes = true;
907  NumberOfTries++;
908  if ( NumberOfTries == 100*(NumberOfTries/100) ) {
909  // At large number of tries it would be better to reduce the values of Pt's
910  ScaleFactor /= 2.0;
911  }
912 
913  if ( Alfa_R == 1.0 ) {
914  G4double Xaq1 = std::sqrt( G4UniformRand() );
915  G4double Xaq2 = 1.0 - Xaq1;
916  Quark_Mom[0].setZ( Xaq1 ); Quark_Mom[1].setZ( Xaq2 );
917  } else {
918  G4double Xaq1 = sqr( std::sin( pi/2.0*G4UniformRand() ) );
919  G4double Xaq2 = 1.0 - Xaq1;
920  Quark_Mom[0].setZ( Xaq1 ); Quark_Mom[1].setZ( Xaq2 );
921  }
922 
923  // Sampling X's of baryon ------------
924  if ( Alfa_R == 1.0 ) {
925  G4double Xq1 = 1.0 - std::sqrt( G4UniformRand() );
926  G4double Xq2 = 1.0 - Xq1;
927  Quark_Mom[2].setZ( Xq1 ); Quark_Mom[3].setZ( Xq2 );
928  } else {
929  G4double Xq1 = sqr( std::sin( pi/2.0*G4UniformRand() ) );
930  G4double Xq2 = 1.0 - Xq1;
931  Quark_Mom[2].setZ( Xq1 ); Quark_Mom[3].setZ( Xq2 );
932  }
933 
934  G4double Alfa( 0.0 ), Beta( 0.0 );
935  for ( G4int i = 0; i < 2; i++ ) { // For Anti-baryon
936  if ( Quark_Mom[i].getZ() != 0.0 ) {
937  Alfa += ( ScaleFactor * ModMom2[i] + MassQ2 ) / Quark_Mom[i].getZ();
938  } else {
939  Succes = false;
940  }
941  }
942  for ( G4int i = 2; i < 4; i++ ) { // For baryon
943  if ( Quark_Mom[i].getZ() != 0.0 ) {
944  Beta += ( ScaleFactor * ModMom2[i] + MassQ2 ) / Quark_Mom[i].getZ();
945  } else {
946  Succes = false;
947  }
948  }
949 
950  if ( ! Succes ) continue;
951 
952  if ( std::sqrt( Alfa ) + std::sqrt( Beta ) > SqrtS ) {
953  Succes = false;
954  continue;
955  }
956 
957  G4double DecayMomentum2 = S*S + Alfa*Alfa + Beta*Beta
958  - 2.0*S*Alfa - 2.0*S*Beta - 2.0*Alfa*Beta;
959  WminusTarget = ( S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS;
960  WplusProjectile = SqrtS - Beta/WminusTarget;
961 
962  } while ( ( ! Succes ) &&
963  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
964  if ( loopCounter >= maxNumberOfLoops ) {
965  return false;
966  }
967 
968  G4double SqrtScaleF = std::sqrt( ScaleFactor );
969 
970  for ( G4int i = 0; i < 2; i++ ) {
971  G4double Pz = WplusProjectile * Quark_Mom[i].getZ() / 2.0 -
972  ( ScaleFactor * ModMom2[i] + MassQ2 ) /
973  ( 2.0 * WplusProjectile * Quark_Mom[i].getZ() );
974  Quark_Mom[i].setZ( Pz );
975  if ( ScaleFactor != 1.0 ) {
976  Quark_Mom[i].setX( SqrtScaleF * Quark_Mom[i].getX() );
977  Quark_Mom[i].setY( SqrtScaleF * Quark_Mom[i].getY() );
978  }
979  //G4cout << "Anti Q " << i << " " << Quark_Mom[i] << G4endl;
980  }
981  for ( G4int i = 2; i < 4; i++ ) {
982  G4double Pz = -WminusTarget * Quark_Mom[i].getZ() / 2.0 +
983  ( ScaleFactor * ModMom2[i] + MassQ2 ) /
984  ( 2.0 * WminusTarget * Quark_Mom[i].getZ() );
985  Quark_Mom[i].setZ( Pz );
986  if ( ScaleFactor != 1.0 ) {
987  Quark_Mom[i].setX( SqrtScaleF * Quark_Mom[i].getX() );
988  Quark_Mom[i].setY( SqrtScaleF * Quark_Mom[i].getY() );
989  }
990  //G4cout << "Bary Q " << i << " " << Quark_Mom[i] << G4endl;
991  }
992  //G4cout << "Sum AQ " << Quark_Mom[0] + Quark_Mom[1] << G4endl
993  // << "Sum Q " << Quark_Mom[2] + Quark_Mom[3] << G4endl;
994 
995  G4ThreeVector tmp = Quark_Mom[0] + Quark_Mom[2];
996  G4LorentzVector Pstring1( tmp, std::sqrt( Quark_Mom[0].mag2() + MassQ2 ) +
997  std::sqrt( Quark_Mom[2].mag2() + MassQ2 ) );
998  G4double Ystring1 = Pstring1.rapidity();
999 
1000  //G4cout << "Mom 1 string " << G4endl << Quark_Mom[0] << G4endl << Quark_Mom[2] << G4endl
1001  // << tmp << " " << tmp.mag() << G4endl;
1002  //G4cout << "1 str " << Pstring1 << " " << Pstring1.mag() << " " << Ystring1 << G4endl;
1003 
1004  tmp = Quark_Mom[1] + Quark_Mom[3];
1005  G4LorentzVector Pstring2( tmp, std::sqrt( Quark_Mom[1].mag2() + MassQ2 ) +
1006  std::sqrt( Quark_Mom[3].mag2() + MassQ2 ) );
1007  G4double Ystring2 = Pstring2.rapidity();
1008 
1009  //G4cout << "Mom 2 string " << G4endl <<Quark_Mom[1] << G4endl << Quark_Mom[3] << G4endl
1010  // << tmp << " " << tmp.mag() << G4endl;
1011  //G4cout << "2 str " << Pstring2 << " " << Pstring2.mag() << " " << Ystring2 << G4endl;
1012 
1013  if ( Ystring1 > Ystring2 ) {
1014  Pprojectile = Pstring1;
1015  Ptarget = Pstring2;
1016  } else {
1017  Pprojectile = Pstring2;
1018  Ptarget = Pstring1;
1019  }
1020 
1021  //G4cout << "SumP CMS " << Pprojectile + Ptarget << " " << SqrtS << G4endl;
1022  Pprojectile.transform( toLab );
1023  Ptarget.transform( toLab );
1024  //G4cout << " SumP Lab " << Pprojectile + Ptarget << " " << SqrtS << G4endl;
1025 
1026  // Calculation of the creation time
1027  projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
1028  projectile->SetPosition( target->GetPosition() );
1029  // Creation time and position of target nucleon were determined in
1030  // ReggeonCascade() of G4FTFModel
1031  //G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl;
1032  projectile->Set4Momentum( Pprojectile );
1033  target->Set4Momentum( Ptarget );
1034  projectile->IncrementCollisionCount( 1 );
1035  target->IncrementCollisionCount( 1 );
1036 
1037 theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016
1038 
1039  return true;
1040 
1041  } // End of if ( CandidatsN != 0 )
1042 
1043  } // End of if ( Ksi < ( X_a + X_b + X_c ) / Xannihilation )
1044 
1045  // Simulation of anti-quark-quark string creation
1046 
1047  if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xannihilation ) {
1048 
1049  #ifdef debugFTFannih
1050  G4cout << "Process d, only 1 quark - anti-quark string" << G4endl;
1051  #endif
1052 
1053  G4int CandidatsN( 0 ), CandAQ[36], CandQ[36];
1054  G4int LeftAQ( 0 ), LeftQ( 0 );
1055 
1056  for ( G4int iAQ1 = 0; iAQ1 < 3; iAQ1++ ) {
1057  for ( G4int iAQ2 = 0; iAQ2 < 3; iAQ2++ ) {
1058  if ( iAQ1 != iAQ2 ) {
1059  for ( G4int iQ1 = 0; iQ1 < 3; iQ1++ ) {
1060  for ( G4int iQ2 = 0; iQ2 < 3; iQ2++ ) {
1061  if ( iQ1 != iQ2 ) {
1062  if ( -AQ[iAQ1] == Q[iQ1] && -AQ[iAQ2] == Q[iQ2] ) {
1063  if ( iAQ1 == 0 && iAQ2 == 1 ) { CandAQ[CandidatsN] = 2; }
1064  if ( iAQ1 == 1 && iAQ2 == 0 ) { CandAQ[CandidatsN] = 2; }
1065 
1066  if ( iAQ1 == 0 && iAQ2 == 2 ) { CandAQ[CandidatsN] = 1; }
1067  if ( iAQ1 == 2 && iAQ2 == 0 ) { CandAQ[CandidatsN] = 1; }
1068 
1069  if ( iAQ1 == 1 && iAQ2 == 2 ) { CandAQ[CandidatsN] = 0; }
1070  if ( iAQ1 == 2 && iAQ2 == 1 ) { CandAQ[CandidatsN] = 0; }
1071 
1072  if ( iQ1 == 0 && iQ2 == 1 ) { CandQ[CandidatsN] = 2; }
1073  if ( iQ1 == 1 && iQ2 == 0 ) { CandQ[CandidatsN] = 2; }
1074 
1075  if ( iQ1 == 0 && iQ2 == 2 ) { CandQ[CandidatsN] = 1; }
1076  if ( iQ1 == 2 && iQ2 == 0 ) { CandQ[CandidatsN] = 1; }
1077 
1078  if ( iQ1 == 1 && iQ2 == 2 ) { CandQ[CandidatsN] = 0; }
1079  if ( iQ1 == 2 && iQ2 == 1 ) { CandQ[CandidatsN] = 0; }
1080  CandidatsN++;
1081  }
1082  }
1083  }
1084  }
1085  }
1086  }
1087  }
1088 
1089  if ( CandidatsN != 0 ) {
1090  G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
1091  LeftAQ = AQ[ CandAQ[SampledCase] ];
1092  LeftQ = Q[ CandQ[SampledCase] ];
1093  //G4cout << "Left Aq Q " << LeftAQ << " " << LeftQ << G4endl;
1094 
1095  // Set the string properties
1096  projectile->SplitUp();
1097  //projectile->SetFirstParton( LeftAQ );
1098  //projectile->SetSecondParton( LeftQ );
1099  projectile->SetFirstParton( LeftQ );
1100  projectile->SetSecondParton( LeftAQ );
1101  projectile->SetStatus( 0 );
1102 
1103 // Uzhi March 2016 start
1104 G4int aAQ, aQ;
1105 aAQ=std::abs( LeftAQ ); aQ=std::abs( LeftQ );
1106 
1107 G4int NewCode;
1108 G4double aKsi = G4UniformRand();
1109 
1110 if ( aAQ == aQ )
1111 {
1112  if ( aAQ != 3 )
1113  {
1114  NewCode = 111; // Pi0-meson
1115  if ( aKsi < 0.5 )
1116  {
1117  NewCode = 221; // Eta -meson
1118  if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson
1119  }
1120  } else
1121  {
1122  NewCode = 221; // Eta -meson
1123  if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson
1124  }
1125 } else
1126 {
1127  if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ; }
1128  else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ; }
1129 }
1130 
1131 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
1132 if(!TestParticle) return false;
1133 projectile->SetDefinition( TestParticle );
1134 theParameters->SetProjMinDiffMass( 0.5 ); // (0.5) // GeV Uzhi March 2016
1135 theParameters->SetProjMinNonDiffMass( 0.5 );
1136 // Uzhi March 2016 end
1137 
1138  target->SetStatus( 4 ); // The target nucleon has annihilated 3->4 Uzhi Oct 2014
1139  Pprojectile.setPx( 0.0 );
1140  Pprojectile.setPy( 0.0 );
1141  Pprojectile.setPz( 0.0 );
1142  Pprojectile.setE( SqrtS );
1143  Pprojectile.transform( toLab );
1144 
1145  // Calculation of the creation time
1146  projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
1147  projectile->SetPosition( target->GetPosition() );
1148  // Creation time and position of target nucleon were determined in
1149  // ReggeonCascade() of G4FTFModel
1150 
1151  //G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl;
1152  projectile->Set4Momentum( Pprojectile );
1153 
1154  projectile->IncrementCollisionCount( 1 );
1155  target->IncrementCollisionCount( 1 );
1156 
1157 theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016
1158 
1159  return true;
1160  }
1161 
1162  } // End of if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xannihilation )
1163 
1164  //G4cout << "Pr Y " << Pprojectile.rapidity() << " Tr Y " << Ptarget.rapidity() << G4endl;
1165  return true;
1166 }
1167 
1168 
1169 //============================================================================
1170 
1171 G4double G4FTFAnnihilation::ChooseX( G4double /* Alpha */, G4double /* Beta */ ) const {
1172  // If for sampling Xs other values of Alfa and Beta instead of 0.5 will be
1173  // chosen the method will be implemented
1174  //G4double tmp = Alpha*Beta;
1175  //tmp *= 1.0;
1176  return 0.5;
1177 }
1178 
1179 
1180 
1181 //============================================================================
1182 
1183 G4ThreeVector G4FTFAnnihilation::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1184  // @@ this method is used in FTFModel as well. Should go somewhere common!
1185  G4double Pt2( 0.0 );
1186  if ( AveragePt2 <= 0.0 ) {
1187  Pt2 = 0.0;
1188  } else {
1189  Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1190  ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
1191  }
1192  G4double Pt = std::sqrt( Pt2 );
1193  G4double phi = G4UniformRand() * twopi;
1194  return G4ThreeVector ( Pt*std::cos( phi ), Pt*std::sin( phi ), 0.0 );
1195 }
1196 
1197 
1198 //============================================================================
1199 
1200 void G4FTFAnnihilation::UnpackBaryon( G4int IdPDG, G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1201  G4int AbsId = std::abs( IdPDG );
1202  Q1 = AbsId / 1000;
1203  Q2 = ( AbsId % 1000 ) / 100;
1204  Q3 = ( AbsId % 100 ) / 10;
1205  if ( IdPDG < 0 ) { Q1 = -Q1; Q2 = -Q2; Q3 = -Q3; } // Anti-baryon
1206  return;
1207 }
1208 
1209 
1210 //============================================================================
1211 
1212 G4FTFAnnihilation::G4FTFAnnihilation( const G4FTFAnnihilation& ) {
1213  throw G4HadronicException( __FILE__, __LINE__,
1214  "G4FTFAnnihilation copy contructor not meant to be called" );
1215 }
1216 
1217 
1218 //============================================================================
1219 
1220 const G4FTFAnnihilation & G4FTFAnnihilation::operator=( const G4FTFAnnihilation& ) {
1221  throw G4HadronicException( __FILE__, __LINE__,
1222  "G4FTFAnnihilation = operator not meant to be called" );
1223 }
1224 
1225 
1226 //============================================================================
1227 
1228 int G4FTFAnnihilation::operator==( const G4FTFAnnihilation& ) const {
1229  throw G4HadronicException( __FILE__, __LINE__,
1230  "G4FTFAnnihilation == operator not meant to be called" );
1231 }
1232 
1233 
1234 //============================================================================
1235 
1236 int G4FTFAnnihilation::operator!=( const G4FTFAnnihilation& ) const {
1237  throw G4HadronicException( __FILE__, __LINE__,
1238  "G4DiffractiveExcitation != operator not meant to be called" );
1239 }
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
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:566
G4FTFAnnihilation::GaussianPt
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
Definition: G4FTFAnnihilation.cc:1183
G4FTFAnnihilation::operator!=
int operator!=(const G4FTFAnnihilation &right) const
Definition: G4FTFAnnihilation.cc:1236
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:1228
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:1200
G4FTFParameters::SetTarMinNonDiffMass
void SetTarMinNonDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:302
G4VSplitableHadron::SetDefinition
void SetDefinition(const G4ParticleDefinition *aDefinition)
Definition: G4VSplitableHadron.hh:137
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
G4FTFParameters::SetProjMinNonDiffMass
void SetProjMinNonDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:294
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
G4FTFParameters::SetProjMinDiffMass
void SetProjMinDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:290
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
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:96
G4FTFAnnihilation::ChooseX
G4double ChooseX(G4double Alpha, G4double Beta) const
Definition: G4FTFAnnihilation.cc:1171
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
G4FTFParameters::SetProbabilityOfAnnihilation
void SetProbabilityOfAnnihilation(const G4double aValue)
Definition: G4FTFParameters.hh:249
G4SystemOfUnits.hh
G4ThreeVector.hh
G4Pow.hh
G4FTFAnnihilation::operator=
const G4FTFAnnihilation & operator=(const G4FTFAnnihilation &right)
Definition: G4FTFAnnihilation.cc:1220
G4FTFParameters::SetTarMinDiffMass
void SetTarMinDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:298
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4LorentzRotation.hh