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