ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DiffuseElastic.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DiffuseElastic.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 
42 #ifndef G4DiffuseElastic_h
43 #define G4DiffuseElastic_h 1
44 
45 #include <CLHEP/Units/PhysicalConstants.h>
46 #include "globals.hh"
47 #include "G4HadronElastic.hh"
48 #include "G4HadProjectile.hh"
49 #include "G4Nucleus.hh"
50 
51 #include "G4Pow.hh"
52 
53 
54 class G4ParticleDefinition;
55 class G4PhysicsTable;
56 class G4PhysicsLogVector;
57 
58 class G4DiffuseElastic : public G4HadronElastic // G4HadronicInteraction
59 {
60 public:
61 
62  G4DiffuseElastic();
63 
64  // G4DiffuseElastic(const G4ParticleDefinition* aParticle);
65 
66 
67 
68 
69 
70  virtual ~G4DiffuseElastic();
71 
72  virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/,
73  G4Nucleus & /*targetNucleus*/);
74 
75  void Initialise();
76 
77  void InitialiseOnFly(G4double Z, G4double A);
78 
79  void BuildAngleTable();
80 
81 
82  // G4HadFinalState* ApplyYourself(const G4HadProjectile & aTrack, G4Nucleus & targetNucleus);
83 
84  virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
85  G4double plab,
86  G4int Z, G4int A);
87 
88  G4double NeutronTuniform(G4int Z);
89 
90  void SetPlabLowLimit(G4double value);
91 
92  void SetHEModelLowLimit(G4double value);
93 
94  void SetQModelLowLimit(G4double value);
95 
96  void SetLowestEnergyLimit(G4double value);
97 
98  void SetRecoilKinEnergyLimit(G4double value);
99 
100  G4double SampleT(const G4ParticleDefinition* aParticle,
101  G4double p, G4double A);
102 
103  G4double SampleTableT(const G4ParticleDefinition* aParticle,
104  G4double p, G4double Z, G4double A);
105 
106  G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
107 
108  G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
109  G4double Z, G4double A);
110 
111  G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
112 
113  G4double SampleThetaLab(const G4HadProjectile* aParticle,
114  G4double tmass, G4double A);
115 
116  G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
117  G4double theta,
118  G4double momentum,
119  G4double A );
120 
121  G4double GetInvElasticXsc( const G4ParticleDefinition* particle,
122  G4double theta,
123  G4double momentum,
124  G4double A, G4double Z );
125 
126  G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
127  G4double theta,
128  G4double momentum,
129  G4double A, G4double Z );
130 
131  G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle,
132  G4double tMand,
133  G4double momentum,
134  G4double A, G4double Z );
135 
136  G4double IntegralElasticProb( const G4ParticleDefinition* particle,
137  G4double theta,
138  G4double momentum,
139  G4double A );
140 
141 
142  G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle,
143  G4double theta,
144  G4double momentum,
145  G4double Z );
146 
147  G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
148  G4double tMand,
149  G4double momentum,
150  G4double A, G4double Z );
151 
152  G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,
153  G4double momentum, G4double Z );
154 
155  G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
156  G4double momentum, G4double Z,
157  G4double theta1, G4double theta2 );
158 
159 
160  G4double CalculateParticleBeta( const G4ParticleDefinition* particle,
161  G4double momentum );
162 
163  G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
164 
165  G4double CalculateAm( G4double momentum, G4double n, G4double Z);
166 
167  G4double CalculateNuclearRad( G4double A);
168 
169  G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
170  G4double tmass, G4double thetaCMS);
171 
172  G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
173  G4double tmass, G4double thetaLab);
174 
175  void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
176  G4double Z, G4double A);
177 
178 
179 
180  G4double BesselJzero(G4double z);
181  G4double BesselJone(G4double z);
182  G4double DampFactor(G4double z);
183  G4double BesselOneByArg(G4double z);
184 
185  G4double GetDiffElasticProb(G4double theta);
186  G4double GetDiffElasticSumProb(G4double theta);
187  G4double GetDiffElasticSumProbA(G4double alpha);
188  G4double GetIntegrandFunction(G4double theta);
189 
190 
191  G4double GetNuclearRadius(){return fNuclearRadius;};
192 
193 private:
194 
195 
196  G4ParticleDefinition* theProton;
197  G4ParticleDefinition* theNeutron;
198  G4ParticleDefinition* theDeuteron;
199  G4ParticleDefinition* theAlpha;
200 
201  const G4ParticleDefinition* thePionPlus;
202  const G4ParticleDefinition* thePionMinus;
203 
204  G4double lowEnergyRecoilLimit;
205  G4double lowEnergyLimitHE;
206  G4double lowEnergyLimitQ;
207  G4double lowestEnergyLimit;
208  G4double plabLowLimit;
209 
210  G4int fEnergyBin;
211  G4int fAngleBin;
212 
213  G4PhysicsLogVector* fEnergyVector;
214  G4PhysicsTable* fAngleTable;
215  std::vector<G4PhysicsTable*> fAngleBank;
216 
217  std::vector<G4double> fElementNumberVector;
218  std::vector<G4String> fElementNameVector;
219 
220  const G4ParticleDefinition* fParticle;
221  G4double fWaveVector;
222  G4double fAtomicWeight;
223  G4double fAtomicNumber;
224  G4double fNuclearRadius;
225  G4double fBeta;
226  G4double fZommerfeld;
227  G4double fAm;
228  G4bool fAddCoulomb;
229 
230 };
231 
232 inline G4bool G4DiffuseElastic::IsApplicable(const G4HadProjectile & projectile,
233  G4Nucleus & nucleus)
234 {
235  if( ( projectile.GetDefinition() == G4Proton::Proton() ||
236  projectile.GetDefinition() == G4Neutron::Neutron() ||
237  projectile.GetDefinition() == G4PionPlus::PionPlus() ||
238  projectile.GetDefinition() == G4PionMinus::PionMinus() ||
239  projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
240  projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
241 
242  nucleus.GetZ_asInt() >= 2 ) return true;
243  else return false;
244 }
245 
246 inline void G4DiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
247 {
248  lowEnergyRecoilLimit = value;
249 }
250 
251 inline void G4DiffuseElastic::SetPlabLowLimit(G4double value)
252 {
253  plabLowLimit = value;
254 }
255 
256 inline void G4DiffuseElastic::SetHEModelLowLimit(G4double value)
257 {
258  lowEnergyLimitHE = value;
259 }
260 
261 inline void G4DiffuseElastic::SetQModelLowLimit(G4double value)
262 {
263  lowEnergyLimitQ = value;
264 }
265 
266 inline void G4DiffuseElastic::SetLowestEnergyLimit(G4double value)
267 {
268  lowestEnergyLimit = value;
269 }
270 
271 
273 //
274 // Bessel J0 function based on rational approximation from
275 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
276 
277 inline G4double G4DiffuseElastic::BesselJzero(G4double value)
278 {
279  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
280 
281  modvalue = std::fabs(value);
282 
283  if ( value < 8.0 && value > -8.0 )
284  {
285  value2 = value*value;
286 
287  fact1 = 57568490574.0 + value2*(-13362590354.0
288  + value2*( 651619640.7
289  + value2*(-11214424.18
290  + value2*( 77392.33017
291  + value2*(-184.9052456 ) ) ) ) );
292 
293  fact2 = 57568490411.0 + value2*( 1029532985.0
294  + value2*( 9494680.718
295  + value2*(59272.64853
296  + value2*(267.8532712
297  + value2*1.0 ) ) ) );
298 
299  bessel = fact1/fact2;
300  }
301  else
302  {
303  arg = 8.0/modvalue;
304 
305  value2 = arg*arg;
306 
307  shift = modvalue-0.785398164;
308 
309  fact1 = 1.0 + value2*(-0.1098628627e-2
310  + value2*(0.2734510407e-4
311  + value2*(-0.2073370639e-5
312  + value2*0.2093887211e-6 ) ) );
313 
314  fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
315  + value2*(-0.6911147651e-5
316  + value2*(0.7621095161e-6
317  - value2*0.934945152e-7 ) ) );
318 
319  bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
320  }
321  return bessel;
322 }
323 
325 //
326 // Bessel J1 function based on rational approximation from
327 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
328 
329 inline G4double G4DiffuseElastic::BesselJone(G4double value)
330 {
331  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
332 
333  modvalue = std::fabs(value);
334 
335  if ( modvalue < 8.0 )
336  {
337  value2 = value*value;
338 
339  fact1 = value*(72362614232.0 + value2*(-7895059235.0
340  + value2*( 242396853.1
341  + value2*(-2972611.439
342  + value2*( 15704.48260
343  + value2*(-30.16036606 ) ) ) ) ) );
344 
345  fact2 = 144725228442.0 + value2*(2300535178.0
346  + value2*(18583304.74
347  + value2*(99447.43394
348  + value2*(376.9991397
349  + value2*1.0 ) ) ) );
350  bessel = fact1/fact2;
351  }
352  else
353  {
354  arg = 8.0/modvalue;
355 
356  value2 = arg*arg;
357 
358  shift = modvalue - 2.356194491;
359 
360  fact1 = 1.0 + value2*( 0.183105e-2
361  + value2*(-0.3516396496e-4
362  + value2*(0.2457520174e-5
363  + value2*(-0.240337019e-6 ) ) ) );
364 
365  fact2 = 0.04687499995 + value2*(-0.2002690873e-3
366  + value2*( 0.8449199096e-5
367  + value2*(-0.88228987e-6
368  + value2*0.105787412e-6 ) ) );
369 
370  bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
371 
372  if (value < 0.0) bessel = -bessel;
373  }
374  return bessel;
375 }
376 
378 //
379 // damp factor in diffraction x/sh(x), x was already *pi
380 
381 inline G4double G4DiffuseElastic::DampFactor(G4double x)
382 {
383  G4double df;
384  G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
385 
386  // x *= pi;
387 
388  if( std::fabs(x) < 0.01 )
389  {
390  df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
391  }
392  else
393  {
394  df = x/std::sinh(x);
395  }
396  return df;
397 }
398 
399 
401 //
402 // return J1(x)/x with special case for small x
403 
404 inline G4double G4DiffuseElastic::BesselOneByArg(G4double x)
405 {
406  G4double x2, result;
407 
408  if( std::fabs(x) < 0.01 )
409  {
410  x *= 0.5;
411  x2 = x*x;
412  result = 2. - x2 + x2*x2/6.;
413  }
414  else
415  {
416  result = BesselJone(x)/x;
417  }
418  return result;
419 }
420 
422 //
423 // return particle beta
424 
425 inline G4double G4DiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle,
426  G4double momentum )
427 {
428  G4double mass = particle->GetPDGMass();
429  G4double a = momentum/mass;
430  fBeta = a/std::sqrt(1+a*a);
431 
432  return fBeta;
433 }
434 
436 //
437 // return Zommerfeld parameter for Coulomb scattering
438 
439 inline G4double G4DiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
440 {
441  fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
442 
443  return fZommerfeld;
444 }
445 
447 //
448 // return Wentzel correction for Coulomb scattering
449 
450 inline G4double G4DiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
451 {
452  G4double k = momentum/CLHEP::hbarc;
453  G4double ch = 1.13 + 3.76*n*n;
454  G4double zn = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
455  G4double zn2 = zn*zn;
456  fAm = ch/zn2;
457 
458  return fAm;
459 }
460 
462 //
463 // calculate nuclear radius for different atomic weights using different approximations
464 
465 inline G4double G4DiffuseElastic::CalculateNuclearRad( G4double A)
466 {
467  G4double R, r0, a11, a12, a13, a2, a3;
468 
469  a11 = 1.26; // 1.08, 1.16
470  a12 = 1.; // 1.08, 1.16
471  a13 = 1.12; // 1.08, 1.16
472  a2 = 1.1;
473  a3 = 1.;
474 
475  // Special rms radii for light nucleii
476 
477  if (A < 50.)
478  {
479  if (std::abs(A-1.) < 0.5) return 0.89*CLHEP::fermi; // p
480  else if(std::abs(A-2.) < 0.5) return 2.13*CLHEP::fermi; // d
481  else if( // std::abs(Z-1.) < 0.5 &&
482 std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi; // t
483 
484  // else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96CLHEP::fermi; // He3
485  else if( // std::abs(Z-2.) < 0.5 &&
486 std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi; // He4
487 
488  else if( // std::abs(Z-3.) < 0.5
489  std::abs(A-7.) < 0.5 ) return 2.40*CLHEP::fermi; // Li7
490  else if( // std::abs(Z-4.) < 0.5
491 std::abs(A-9.) < 0.5) return 2.51*CLHEP::fermi; // Be9
492 
493  else if( 10. < A && A <= 16. ) r0 = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; // 1.08CLHEP::fermi;
494  else if( 15. < A && A <= 20. ) r0 = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
495  else if( 20. < A && A <= 30. ) r0 = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
496  else r0 = a2*CLHEP::fermi;
497 
498  R = r0*G4Pow::GetInstance()->A13(A);
499  }
500  else
501  {
502  r0 = a3*CLHEP::fermi;
503 
504  R = r0*G4Pow::GetInstance()->powA(A, 0.27);
505  }
506  fNuclearRadius = R;
507  return R;
508  /*
509  G4double r0;
510  if( A < 50. )
511  {
512  if( A > 10. ) r0 = 1.16*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; // 1.08CLHEP::fermi;
513  else r0 = 1.1*CLHEP::fermi;
514  fNuclearRadius = r0*G4Pow::GetInstance()->A13(A);
515  }
516  else
517  {
518  r0 = 1.7*CLHEP::fermi; // 1.7*CLHEP::fermi;
519  fNuclearRadius = r0*G4Pow::GetInstance()->powA(A, 0.27); // 0.27);
520  }
521  return fNuclearRadius;
522  */
523 }
524 
526 //
527 // return Coulomb scattering differential xsc with Wentzel correction
528 
529 inline G4double G4DiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle,
530  G4double theta,
531  G4double momentum,
532  G4double Z )
533 {
534  G4double sinHalfTheta = std::sin(0.5*theta);
535  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
536  G4double beta = CalculateParticleBeta( particle, momentum);
537  G4double z = particle->GetPDGCharge();
538  G4double n = CalculateZommerfeld( beta, z, Z );
539  G4double am = CalculateAm( momentum, n, Z);
540  G4double k = momentum/CLHEP::hbarc;
541  G4double ch = 0.5*n/k;
542  G4double ch2 = ch*ch;
543  G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
544 
545  return xsc;
546 }
547 
548 
550 //
551 // return Coulomb scattering total xsc with Wentzel correction
552 
553 inline G4double G4DiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,
554  G4double momentum, G4double Z )
555 {
556  G4double beta = CalculateParticleBeta( particle, momentum);
557  G4cout<<"beta = "<<beta<<G4endl;
558  G4double z = particle->GetPDGCharge();
559  G4double n = CalculateZommerfeld( beta, z, Z );
560  G4cout<<"fZomerfeld = "<<n<<G4endl;
561  G4double am = CalculateAm( momentum, n, Z);
562  G4cout<<"cof Am = "<<am<<G4endl;
563  G4double k = momentum/CLHEP::hbarc;
564  G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
565  G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
566  G4double ch = n/k;
567  G4double ch2 = ch*ch;
568  G4double xsc = ch2*CLHEP::pi/(am +am*am);
569 
570  return xsc;
571 }
572 
574 //
575 // return Coulomb scattering xsc with Wentzel correction integrated between
576 // theta1 and < theta2
577 
578 inline G4double G4DiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
579  G4double momentum, G4double Z,
580  G4double theta1, G4double theta2 )
581 {
582  G4double c1 = std::cos(theta1);
583  G4cout<<"c1 = "<<c1<<G4endl;
584  G4double c2 = std::cos(theta2);
585  G4cout<<"c2 = "<<c2<<G4endl;
586  G4double beta = CalculateParticleBeta( particle, momentum);
587  // G4cout<<"beta = "<<beta<<G4endl;
588  G4double z = particle->GetPDGCharge();
589  G4double n = CalculateZommerfeld( beta, z, Z );
590  // G4cout<<"fZomerfeld = "<<n<<G4endl;
591  G4double am = CalculateAm( momentum, n, Z);
592  // G4cout<<"cof Am = "<<am<<G4endl;
593  G4double k = momentum/CLHEP::hbarc;
594  // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
595  // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
596  G4double ch = n/k;
597  G4double ch2 = ch*ch;
598  am *= 2.;
599  G4double xsc = ch2*CLHEP::twopi*(c1-c2);
600  xsc /= (1 - c1 + am)*(1 - c2 + am);
601 
602  return xsc;
603 }
604 
605 #endif