ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HadronElastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4HadronElastic.cc
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 // Geant4 Header : G4HadronElastic
28 //
29 // Author : V.Ivanchenko 29 June 2009 (redesign old elastic model)
30 //
31 
32 #include "G4HadronElastic.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4ParticleTable.hh"
35 #include "G4ParticleDefinition.hh"
36 #include "G4IonTable.hh"
37 #include "Randomize.hh"
38 #include "G4Proton.hh"
39 #include "G4Neutron.hh"
40 #include "G4Deuteron.hh"
41 #include "G4Alpha.hh"
42 #include "G4Pow.hh"
43 #include "G4Exp.hh"
44 #include "G4Log.hh"
45 #include "G4HadronicParameters.hh"
46 
47 
49  : G4HadronicInteraction(name)
50 {
51  SetMinEnergy( 0.0*GeV );
53  lowestEnergyLimit= 1.e-6*eV;
54  pLocalTmax = 0.0;
55  nwarn = 0;
56 
61 }
62 
64 {}
65 
66 
67 void G4HadronElastic::ModelDescription(std::ostream& outFile) const
68 {
69  outFile << "G4HadronElastic is the base class for all hadron-nucleus\n"
70  << "elastic scattering models except HP.\n"
71  << "By default it uses the Gheisha two-exponential momentum\n"
72  << "transfer parameterization. The model is fully relativistic\n"
73  << "as opposed to the original Gheisha model which was not.\n"
74  << "This model may be used for all long-lived hadrons at all\n"
75  << "incident energies but fit the data only for relativistic scattering.\n";
76 }
77 
79  const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
80 {
82 
83  const G4HadProjectile* aParticle = &aTrack;
84  G4double ekin = aParticle->GetKineticEnergy();
85 
86  // no scattering below the limit
87  if(ekin <= lowestEnergyLimit) {
90  return &theParticleChange;
91  }
92 
93  G4int A = targetNucleus.GetA_asInt();
94  G4int Z = targetNucleus.GetZ_asInt();
95 
96  // Scattered particle referred to axis of incident particle
97  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
98  G4double m1 = theParticle->GetPDGMass();
99  G4double plab = std::sqrt(ekin*(ekin + 2.0*m1));
100 
101  if (verboseLevel>1) {
102  G4cout << "G4HadronElastic: "
103  << aParticle->GetDefinition()->GetParticleName()
104  << " Plab(GeV/c)= " << plab/GeV
105  << " Ekin(MeV) = " << ekin/MeV
106  << " scattered off Z= " << Z
107  << " A= " << A
108  << G4endl;
109  }
110 
112  G4double e1 = m1 + ekin;
113  G4LorentzVector lv(0.0,0.0,plab,e1+mass2);
114  G4ThreeVector bst = lv.boostVector();
115  G4double momentumCMS = plab*mass2/std::sqrt(m1*m1 + mass2*mass2 + 2.*mass2*e1);
116 
117  pLocalTmax = 4.0*momentumCMS*momentumCMS;
118 
119  // Sampling in CM system
120  G4double t = SampleInvariantT(theParticle, plab, Z, A);
121 
122  if(t < 0.0 || t > pLocalTmax) {
123  // For the very rare cases where cos(theta) is greater than 1 or smaller than -1,
124  // print some debugging information via a "JustWarning" exception, and resample
125  // using the default algorithm
126 #ifdef G4VERBOSE
127  if(nwarn < 2) {
129  ed << GetModelName() << " wrong sampling t= " << t << " tmax= " << pLocalTmax
130  << " for " << aParticle->GetDefinition()->GetParticleName()
131  << " ekin=" << ekin << " MeV"
132  << " off (Z,A)=(" << Z << "," << A << ") - will be resampled" << G4endl;
133  G4Exception( "G4HadronElastic::ApplyYourself", "hadEla001", JustWarning, ed);
134  ++nwarn;
135  }
136 #endif
137  t = G4HadronElastic::SampleInvariantT(theParticle, plab, Z, A);
138  }
139 
141  G4double cost = 1. - 2.0*t/pLocalTmax;
142 
143  if (cost > 1.0) { cost = 1.0; }
144  else if(cost < -1.0) { cost = -1.0; }
145 
146  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
147 
148  if (verboseLevel>1) {
149  G4cout << " t= " << t << " tmax(GeV^2)= " << pLocalTmax/(GeV*GeV)
150  << " Pcms(GeV)= " << momentumCMS/GeV << " cos(t)=" << cost
151  << " sin(t)=" << sint << G4endl;
152  }
153  G4LorentzVector nlv1(momentumCMS*sint*std::cos(phi),
154  momentumCMS*sint*std::sin(phi),
155  momentumCMS*cost,
156  std::sqrt(momentumCMS*momentumCMS + m1*m1));
157 
158  nlv1.boost(bst);
159 
160  G4double eFinal = nlv1.e() - m1;
161  if (verboseLevel > 1) {
162  G4cout <<"G4HadronElastic: m= " << m1 << " Efin(MeV)= " << eFinal
163  << " 4-M Final: " << nlv1
164  << G4endl;
165  }
166 
167  if(eFinal <= 0.0) {
170  } else {
173  }
174  lv -= nlv1;
175  G4double erec = std::max(lv.e() - mass2, 0.0);
176  if (verboseLevel > 1) {
177  G4cout << "Recoil: " <<" m= " << mass2 << " Erec(MeV)= " << erec
178  << " 4-mom: " << lv
179  << G4endl;
180  }
181 
182  // the recoil is created if kinetic energy above the threshold
183  if(erec > GetRecoilEnergyThreshold()) {
184  G4ParticleDefinition * theDef = nullptr;
185  if(Z == 1 && A == 1) { theDef = theProton; }
186  else if (Z == 1 && A == 2) { theDef = theDeuteron; }
187  else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
188  else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
189  else if (Z == 2 && A == 4) { theDef = theAlpha; }
190  else {
191  theDef =
193  }
194  G4DynamicParticle * aSec = new G4DynamicParticle(theDef, lv.vect().unit(), erec);
196  } else {
198  }
199 
200  return &theParticleChange;
201 }
202 
203 // sample momentum transfer in the CMS system
204 G4double
207 {
208  const G4double plabLowLimit = 400.0*CLHEP::MeV;
209  const G4double GeV2 = GeV*GeV;
210  const G4double z07in13 = std::pow(0.7, 0.3333333333);
211 
212  G4int pdg = std::abs(part->GetPDGEncoding());
213  G4double tmax = pLocalTmax/GeV2;
214 
215  G4double aa, bb, cc, dd;
216  G4Pow* g4pow = G4Pow::GetInstance();
217  if (A <= 62) {
218  if (pdg == 211){ //Pions
219  if(mom >= plabLowLimit){ //High energy
220  bb = 14.5*g4pow->Z23(A);/*14.5*/
221  dd = 10.;
222  cc = 0.075*g4pow->Z13(A)/dd;//1.4
223  //aa = g4pow->powZ(A, 1.93)/bb;//1.63
224  aa = (A*A)/bb;//1.63
225  } else { //Low energy
226  bb = 29.*z07in13*z07in13*g4pow->Z23(A);
227  dd = 15.;
228  cc = 0.04*g4pow->Z13(A)/dd;//1.4
229  aa = g4pow->powZ(A, 1.63)/bb;//1.63
230  }
231  } else { //Other particles
232  bb = 14.5*g4pow->Z23(A);
233  dd = 20.;
234  aa = (A*A)/bb;//1.63
235  cc = 1.4*g4pow->Z13(A)/dd;
236  }
237  //===========================
238  } else { //(A>62)
239  if (pdg == 211) {
240  if(mom >= plabLowLimit){ //high
241  bb = 60.*z07in13*g4pow->Z13(A);//60
242  dd = 30.;
243  aa = 0.5*(A*A)/bb;//1.33
244  cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4 --- 2: 0.4
245  } else { //low
246  bb = 120.*z07in13*g4pow->Z13(A);//60
247  dd = 30.;
248  aa = 2.*g4pow->powZ(A,1.33)/bb;
249  cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4 --- 2: 0.4
250  }
251  } else {
252  bb = 60.*g4pow->Z13(A);
253  dd = 25.;
254  aa = g4pow->powZ(A,1.33)/bb;//1.33
255  cc = 0.2*g4pow->powZ(A,0.4)/dd;//1:0.4 --- 2: 0.4
256  }
257  }
258  G4double q1 = 1.0 - G4Exp(-bb*tmax);
259  G4double q2 = 1.0 - G4Exp(-dd*tmax);
260  G4double s1 = q1*aa;
261  G4double s2 = q2*cc;
262  if((s1 + s2)*G4UniformRand() < s2) {
263  q1 = q2;
264  bb = dd;
265  }
266  return -GeV2*G4Log(1.0 - G4UniformRand()*q1)/bb;
267 }
268 
270 //
271 // Cofs for s-,c-,b-particles ds/dt slopes
272 
274 {
275  // The input parameter "pdg" should be the absolute value of the PDG code
276  // (i.e. the same value for a particle and its antiparticle).
277 
278  G4double coeff = 1.0;
279 
280  // heavy barions
281 
282  static const G4double lBarCof1S = 0.88;
283  static const G4double lBarCof2S = 0.76;
284  static const G4double lBarCof3S = 0.64;
285  static const G4double lBarCof1C = 0.784378;
286  static const G4double lBarCofSC = 0.664378;
287  static const G4double lBarCof2SC = 0.544378;
288  static const G4double lBarCof1B = 0.740659;
289  static const G4double lBarCofSB = 0.620659;
290  static const G4double lBarCof2SB = 0.500659;
291 
292  if( pdg == 3122 || pdg == 3222 || pdg == 3112 || pdg == 3212 )
293  {
294  coeff = lBarCof1S; // Lambda, Sigma+, Sigma-, Sigma0
295 
296  } else if( pdg == 3322 || pdg == 3312 )
297  {
298  coeff = lBarCof2S; // Xi-, Xi0
299  }
300  else if( pdg == 3324)
301  {
302  coeff = lBarCof3S; // Omega
303  }
304  else if( pdg == 4122 || pdg == 4212 || pdg == 4222 || pdg == 4112 )
305  {
306  coeff = lBarCof1C; // LambdaC+, SigmaC+, SigmaC++, SigmaC0
307  }
308  else if( pdg == 4332 )
309  {
310  coeff = lBarCof2SC; // OmegaC
311  }
312  else if( pdg == 4232 || pdg == 4132 )
313  {
314  coeff = lBarCofSC; // XiC+, XiC0
315  }
316  else if( pdg == 5122 || pdg == 5222 || pdg == 5112 || pdg == 5212 )
317  {
318  coeff = lBarCof1B; // LambdaB, SigmaB+, SigmaB-, SigmaB0
319  }
320  else if( pdg == 5332 )
321  {
322  coeff = lBarCof2SB; // OmegaB-
323  }
324  else if( pdg == 5132 || pdg == 5232 ) // XiB-, XiB0
325  {
326  coeff = lBarCofSB;
327  }
328  // heavy mesons Kaons?
329  static const G4double lMesCof1S = 0.82; // Kp/piP kaons?
330  static const G4double llMesCof1C = 0.676568;
331  static const G4double llMesCof1B = 0.610989;
332  static const G4double llMesCof2C = 0.353135;
333  static const G4double llMesCof2B = 0.221978;
334  static const G4double llMesCofSC = 0.496568;
335  static const G4double llMesCofSB = 0.430989;
336  static const G4double llMesCofCB = 0.287557;
337  static const G4double llMesCofEtaP = 0.88;
338  static const G4double llMesCofEta = 0.76;
339 
340  if( pdg == 321 || pdg == 311 || pdg == 310 )
341  {
342  coeff = lMesCof1S; //K+-0
343  }
344  else if( pdg == 511 || pdg == 521 )
345  {
346  coeff = llMesCof1B; // BMeson0, BMeson+
347  }
348  else if(pdg == 421 || pdg == 411 )
349  {
350  coeff = llMesCof1C; // DMeson+, DMeson0
351  }
352  else if( pdg == 531 )
353  {
354  coeff = llMesCofSB; // BSMeson0
355  }
356  else if( pdg == 541 )
357  {
358  coeff = llMesCofCB; // BCMeson+-
359  }
360  else if(pdg == 431 )
361  {
362  coeff = llMesCofSC; // DSMeson+-
363  }
364  else if(pdg == 441 || pdg == 443 )
365  {
366  coeff = llMesCof2C; // Etac, JPsi
367  }
368  else if(pdg == 553 )
369  {
370  coeff = llMesCof2B; // Upsilon
371  }
372  else if(pdg == 221 )
373  {
374  coeff = llMesCofEta; // Eta
375  }
376  else if(pdg == 331 )
377  {
378  coeff = llMesCofEtaP; // Eta'
379  }
380  return coeff;
381 }
382 
383