ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DiffuseElasticV2.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DiffuseElasticV2.hh
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 //
28 // Author: V. Grichine (Vladimir,Grichine@cern.ch)
29 //
30 //
31 // G4 Model: diffuse optical elastic scattering with 4-momentum balance
32 //
33 // Class Description
34 // Final state production model for hadron nuclear elastic scattering;
35 // Class Description - End
36 //
37 //
38 // 24.05.07 V. Grichine, first implementation for hadron (no Coulomb) elastic scattering
39 // 04.09.07 V. Grichine, implementation for Coulomb elastic scattering
40 // 12.06.11 V. Grichine, new interface to G4hadronElastic
41 // 24.11.17 W. Pokorski, code cleanup and performance improvements
42 
43 
44 #ifndef G4DiffuseElasticV2_h
45 #define G4DiffuseElasticV2_h 1
46 
47 #include <CLHEP/Units/PhysicalConstants.h>
48 #include "globals.hh"
49 #include "G4HadronElastic.hh"
50 #include "G4HadProjectile.hh"
51 #include "G4Nucleus.hh"
52 
53 #include "G4Pow.hh"
54 
55 #include <vector>
56 
57 class G4ParticleDefinition;
58 class G4PhysicsTable;
59 class G4PhysicsLogVector;
60 
61 class G4DiffuseElasticV2 : public G4HadronElastic // G4HadronicInteraction
62 {
63 public:
64 
65  G4DiffuseElasticV2();
66 
67  virtual ~G4DiffuseElasticV2();
68 
69  virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/,
70  G4Nucleus & /*targetNucleus*/);
71 
72  void Initialise();
73 
74  void InitialiseOnFly(G4double Z, G4double A);
75 
76  void BuildAngleTable();
77 
78  virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
79  G4double plab,
80  G4int Z, G4int A);
81 
82  G4double NeutronTuniform(G4int Z);
83 
84  void SetPlabLowLimit(G4double value);
85 
86  void SetHEModelLowLimit(G4double value);
87 
88  void SetQModelLowLimit(G4double value);
89 
90  void SetLowestEnergyLimit(G4double value);
91 
92  void SetRecoilKinEnergyLimit(G4double value);
93 
94  G4double SampleTableT(const G4ParticleDefinition* aParticle,
95  G4double p, G4double Z, G4double A);
96 
97  G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
98 
99  G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
100  G4double Z, G4double A);
101 
102  G4double GetScatteringAngle(G4int iMomentum, unsigned long iAngle, G4double position);
103 
104  G4double SampleThetaLab(const G4HadProjectile* aParticle,
105  G4double tmass, G4double A);
106 
107  G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
108 
109  G4double CalculateAm( G4double momentum, G4double n, G4double Z);
110 
111  G4double CalculateNuclearRad( G4double A);
112 
113  G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
114  G4double tmass, G4double thetaCMS);
115 
116  G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
117  G4double tmass, G4double thetaLab);
118 
119 
120  G4double BesselJzero(G4double z);
121  G4double BesselJone(G4double z);
122  G4double DampFactor(G4double z);
123  G4double BesselOneByArg(G4double z);
124 
125  G4double GetDiffElasticSumProbA(G4double alpha);
126  G4double GetIntegrandFunction(G4double theta);
127 
128 
129  G4double GetNuclearRadius(){return fNuclearRadius;};
130 
131 private:
132 
133 
134  G4ParticleDefinition* theProton;
135  G4ParticleDefinition* theNeutron;
136 
137  G4double lowEnergyRecoilLimit;
138  G4double lowEnergyLimitHE;
139  G4double lowEnergyLimitQ;
140  G4double lowestEnergyLimit;
141  G4double plabLowLimit;
142 
143  G4int fEnergyBin;
144  unsigned long fAngleBin;
145 
146  G4PhysicsLogVector* fEnergyVector;
147 
148  std::vector<std::vector<std::vector<double>*>*> fEnergyAngleVectorBank;
149  std::vector<std::vector<std::vector<double>*>*> fEnergySumVectorBank;
150 
151  std::vector<std::vector<double>*>* fEnergyAngleVector;
152  std::vector<std::vector<double>*>* fEnergySumVector;
153 
154 
155  std::vector<G4double> fElementNumberVector;
156  std::vector<G4String> fElementNameVector;
157 
158  const G4ParticleDefinition* fParticle;
159  G4double fWaveVector;
160  G4double fAtomicWeight;
161  G4double fAtomicNumber;
162  G4double fNuclearRadius;
163  G4double fBeta;
164  G4double fZommerfeld;
165  G4double fAm;
166  G4bool fAddCoulomb;
167 
168 };
169 
170 inline G4bool G4DiffuseElasticV2::IsApplicable(const G4HadProjectile & projectile,
171  G4Nucleus & nucleus)
172 {
173  if( ( projectile.GetDefinition() == G4Proton::Proton() ||
174  projectile.GetDefinition() == G4Neutron::Neutron() ||
175  projectile.GetDefinition() == G4PionPlus::PionPlus() ||
176  projectile.GetDefinition() == G4PionMinus::PionMinus() ||
177  projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
178  projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
179 
180  nucleus.GetZ_asInt() >= 2 ) return true;
181  else return false;
182 }
183 
184 inline void G4DiffuseElasticV2::SetRecoilKinEnergyLimit(G4double value)
185 {
186  lowEnergyRecoilLimit = value;
187 }
188 
189 inline void G4DiffuseElasticV2::SetPlabLowLimit(G4double value)
190 {
191  plabLowLimit = value;
192 }
193 
194 inline void G4DiffuseElasticV2::SetHEModelLowLimit(G4double value)
195 {
196  lowEnergyLimitHE = value;
197 }
198 
199 inline void G4DiffuseElasticV2::SetQModelLowLimit(G4double value)
200 {
201  lowEnergyLimitQ = value;
202 }
203 
204 inline void G4DiffuseElasticV2::SetLowestEnergyLimit(G4double value)
205 {
206  lowestEnergyLimit = value;
207 }
208 
209 
211 //
212 // Bessel J0 function based on rational approximation from
213 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
214 
215 inline G4double G4DiffuseElasticV2::BesselJzero(G4double value)
216 {
217  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
218 
219  modvalue = std::fabs(value);
220 
221  if ( value < 8.0 && value > -8.0 )
222  {
223  value2 = value*value;
224 
225  fact1 = 57568490574.0 + value2*(-13362590354.0
226  + value2*( 651619640.7
227  + value2*(-11214424.18
228  + value2*( 77392.33017
229  + value2*(-184.9052456 ) ) ) ) );
230 
231  fact2 = 57568490411.0 + value2*( 1029532985.0
232  + value2*( 9494680.718
233  + value2*(59272.64853
234  + value2*(267.8532712
235  + value2*1.0 ) ) ) );
236 
237  bessel = fact1/fact2;
238  }
239  else
240  {
241  arg = 8.0/modvalue;
242 
243  value2 = arg*arg;
244 
245  shift = modvalue-0.785398164;
246 
247  fact1 = 1.0 + value2*(-0.1098628627e-2
248  + value2*(0.2734510407e-4
249  + value2*(-0.2073370639e-5
250  + value2*0.2093887211e-6 ) ) );
251 
252  fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
253  + value2*(-0.6911147651e-5
254  + value2*(0.7621095161e-6
255  - value2*0.934945152e-7 ) ) );
256 
257  bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
258  }
259  return bessel;
260 }
261 
263 //
264 // Bessel J1 function based on rational approximation from
265 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
266 
267 inline G4double G4DiffuseElasticV2::BesselJone(G4double value)
268 {
269  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
270 
271  modvalue = std::fabs(value);
272 
273  if ( modvalue < 8.0 )
274  {
275  value2 = value*value;
276 
277  fact1 = value*(72362614232.0 + value2*(-7895059235.0
278  + value2*( 242396853.1
279  + value2*(-2972611.439
280  + value2*( 15704.48260
281  + value2*(-30.16036606 ) ) ) ) ) );
282 
283  fact2 = 144725228442.0 + value2*(2300535178.0
284  + value2*(18583304.74
285  + value2*(99447.43394
286  + value2*(376.9991397
287  + value2*1.0 ) ) ) );
288  bessel = fact1/fact2;
289  }
290  else
291  {
292  arg = 8.0/modvalue;
293 
294  value2 = arg*arg;
295 
296  shift = modvalue - 2.356194491;
297 
298  fact1 = 1.0 + value2*( 0.183105e-2
299  + value2*(-0.3516396496e-4
300  + value2*(0.2457520174e-5
301  + value2*(-0.240337019e-6 ) ) ) );
302 
303  fact2 = 0.04687499995 + value2*(-0.2002690873e-3
304  + value2*( 0.8449199096e-5
305  + value2*(-0.88228987e-6
306  + value2*0.105787412e-6 ) ) );
307 
308  bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
309 
310  if (value < 0.0) bessel = -bessel;
311  }
312  return bessel;
313 }
314 
316 //
317 // damp factor in diffraction x/sh(x), x was already *pi
318 
319 inline G4double G4DiffuseElasticV2::DampFactor(G4double x)
320 {
321  G4double df;
322  G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
323 
324  // x *= pi;
325 
326  if( std::fabs(x) < 0.01 )
327  {
328  df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
329  }
330  else
331  {
332  df = x/std::sinh(x);
333  }
334  return df;
335 }
336 
337 
339 //
340 // return J1(x)/x with special case for small x
341 
342 inline G4double G4DiffuseElasticV2::BesselOneByArg(G4double x)
343 {
344  G4double x2, result;
345 
346  if( std::fabs(x) < 0.01 )
347  {
348  x *= 0.5;
349  x2 = x*x;
350  result = 2. - x2 + x2*x2/6.;
351  }
352  else
353  {
354  result = BesselJone(x)/x;
355  }
356  return result;
357 }
358 
359 
361 //
362 // return Zommerfeld parameter for Coulomb scattering
363 
364 inline G4double G4DiffuseElasticV2::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
365 {
366  fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
367 
368  return fZommerfeld;
369 }
370 
372 //
373 // return Wentzel correction for Coulomb scattering
374 
375 inline G4double G4DiffuseElasticV2::CalculateAm( G4double momentum, G4double n, G4double Z)
376 {
377  G4double k = momentum/CLHEP::hbarc;
378  G4double ch = 1.13 + 3.76*n*n;
379  G4double zn = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
380  G4double zn2 = zn*zn;
381  fAm = ch/zn2;
382 
383  return fAm;
384 }
385 
387 //
388 // calculate nuclear radius for different atomic weights using different approximations
389 
390 inline G4double G4DiffuseElasticV2::CalculateNuclearRad( G4double A)
391 {
392  G4double R, r0, a11, a12, a13, a2, a3;
393 
394  a11 = 1.26; // 1.08, 1.16
395  a12 = 1.; // 1.08, 1.16
396  a13 = 1.12; // 1.08, 1.16
397  a2 = 1.1;
398  a3 = 1.;
399 
400  // Special rms radii for light nucleii
401 
402  if (A < 50.)
403  {
404  if (std::abs(A-1.) < 0.5) return 0.89*CLHEP::fermi; // p
405  else if(std::abs(A-2.) < 0.5) return 2.13*CLHEP::fermi; // d
406  else if( // std::abs(Z-1.) < 0.5 &&
407  std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi; // t
408 
409  // else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96CLHEP::fermi; // He3
410  else if( // std::abs(Z-2.) < 0.5 &&
411  std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi; // He4
412 
413  else if( // std::abs(Z-3.) < 0.5
414  std::abs(A-7.) < 0.5 ) return 2.40*CLHEP::fermi; // Li7
415  else if( // std::abs(Z-4.) < 0.5
416  std::abs(A-9.) < 0.5) return 2.51*CLHEP::fermi; // Be9
417 
418  else if( 10. < A && A <= 16. ) r0 = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; // 1.08CLHEP::fermi;
419  else if( 15. < A && A <= 20. ) r0 = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
420  else if( 20. < A && A <= 30. ) r0 = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
421  else r0 = a2*CLHEP::fermi;
422 
423  R = r0*G4Pow::GetInstance()->A13(A);
424  }
425  else
426  {
427  r0 = a3*CLHEP::fermi;
428 
429  R = r0*G4Pow::GetInstance()->powA(A, 0.27);
430  }
431  fNuclearRadius = R;
432 
433  return R;
434 }
435 
436 
437 #endif