include/Geant4/G4DiffuseElasticV2.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // Author: V. Grichine (Vladimir,Grichine@cern.ch)
0029 //
0030 //
0031 // G4 Model: diffuse optical elastic scattering with 4-momentum balance
0032 //
0033 // Class Description
0034 // Final state production model for hadron nuclear elastic scattering;
0035 // Class Description - End
0036 //
0037 //
0038 // 24.05.07 V. Grichine, first implementation for hadron (no Coulomb) elastic scattering
0039 // 04.09.07 V. Grichine, implementation for Coulomb elastic scattering
0040 // 12.06.11 V. Grichine, new interface to G4hadronElastic
0041 // 24.11.17 W. Pokorski, code cleanup and performance improvements
0042
0043
0044 #ifndef G4DiffuseElasticV2_h
0045 #define G4DiffuseElasticV2_h 1
0046
0047 #include <CLHEP/Units/PhysicalConstants.h>
0048 #include "globals.hh"
0049 #include "G4HadronElastic.hh"
0050 #include "G4HadProjectile.hh"
0051 #include "G4Nucleus.hh"
0052
0053 #include "G4Pow.hh"
0054
0055 #include <vector>
0056
0057 class G4ParticleDefinition;
0058 class G4PhysicsTable;
0059 class G4PhysicsLogVector;
0060
0061 class G4DiffuseElasticV2 : public G4HadronElastic // G4HadronicInteraction
0062 {
0063 public:
0064
0065   G4DiffuseElasticV2();
0066
0067   virtual ~G4DiffuseElasticV2();
0068
0069   virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/,
0070                   G4Nucleus & /*targetNucleus*/);
0071
0072   void Initialise();
0073
0074   void InitialiseOnFly(G4double Z, G4double A);
0075
0076   void BuildAngleTable();
0077
0078   virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
0079                     G4double plab,
0080                     G4int Z, G4int A);
0081
0082   G4double NeutronTuniform(G4int Z);
0083
0084   void SetPlabLowLimit(G4double value);
0085
0086   void SetHEModelLowLimit(G4double value);
0087
0088   void SetQModelLowLimit(G4double value);
0089
0090   void SetLowestEnergyLimit(G4double value);
0091
0092   void SetRecoilKinEnergyLimit(G4double value);
0093
0094   G4double SampleTableT(const G4ParticleDefinition* aParticle,
0095                          G4double p, G4double Z, G4double A);
0096
0097   G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
0098
0099   G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
0100                                      G4double Z, G4double A);
0101
0102   G4double GetScatteringAngle(G4int iMomentum, unsigned long iAngle, G4double position);
0103
0104   G4double SampleThetaLab(const G4HadProjectile* aParticle,
0105                                 G4double tmass, G4double A);
0106
0107   G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
0108
0109   G4double CalculateAm( G4double momentum, G4double n, G4double Z);
0110
0111   G4double CalculateNuclearRad( G4double A);
0112
0113   G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
0114                                 G4double tmass, G4double thetaCMS);
0115
0116   G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
0117                                 G4double tmass, G4double thetaLab);
0118
0119
0120   G4double BesselJzero(G4double z);
0121   G4double BesselJone(G4double z);
0122   G4double DampFactor(G4double z);
0123   G4double BesselOneByArg(G4double z);
0124
0125   G4double GetDiffElasticSumProbA(G4double alpha);
0126   G4double GetIntegrandFunction(G4double theta);
0127
0128
0129   G4double GetNuclearRadius(){return fNuclearRadius;};
0130
0131 private:
0132
0133
0134   G4ParticleDefinition* theProton;
0135   G4ParticleDefinition* theNeutron;
0136
0137   G4double lowEnergyRecoilLimit;
0138   G4double lowEnergyLimitHE;
0139   G4double lowEnergyLimitQ;
0140   G4double lowestEnergyLimit;
0141   G4double plabLowLimit;
0142
0143   G4int fEnergyBin;
0144   unsigned long fAngleBin;
0145
0146   G4PhysicsLogVector* fEnergyVector;
0147
0148   std::vector<std::vector<std::vector<double>*>*>  fEnergyAngleVectorBank;
0149   std::vector<std::vector<std::vector<double>*>*>  fEnergySumVectorBank;
0150
0151   std::vector<std::vector<double>*>*  fEnergyAngleVector;
0152   std::vector<std::vector<double>*>*  fEnergySumVector;
0153
0154
0155   std::vector<G4double> fElementNumberVector;
0156   std::vector<G4String> fElementNameVector;
0157
0158   const G4ParticleDefinition* fParticle;
0159   G4double fWaveVector;
0160   G4double fAtomicWeight;
0161   G4double fAtomicNumber;
0162   G4double fNuclearRadius;
0163   G4double fBeta;
0164   G4double fZommerfeld;
0165   G4double fAm;
0166   G4bool fAddCoulomb;
0167
0168 };
0169
0170 inline G4bool G4DiffuseElasticV2::IsApplicable(const G4HadProjectile & projectile,
0171                   G4Nucleus & nucleus)
0172 {
0173   if( ( projectile.GetDefinition() == G4Proton::Proton() ||
0174         projectile.GetDefinition() == G4Neutron::Neutron() ||
0175         projectile.GetDefinition() == G4PionPlus::PionPlus() ||
0176         projectile.GetDefinition() == G4PionMinus::PionMinus() ||
0177         projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
0178         projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
0179
0180         nucleus.GetZ_asInt() >= 2 ) return true;
0181   else                              return false;
0182 }
0183
0184 inline void G4DiffuseElasticV2::SetRecoilKinEnergyLimit(G4double value)
0185 {
0186   lowEnergyRecoilLimit = value;
0187 }
0188
0189 inline void G4DiffuseElasticV2::SetPlabLowLimit(G4double value)
0190 {
0191   plabLowLimit = value;
0192 }
0193
0194 inline void G4DiffuseElasticV2::SetHEModelLowLimit(G4double value)
0195 {
0196   lowEnergyLimitHE = value;
0197 }
0198
0199 inline void G4DiffuseElasticV2::SetQModelLowLimit(G4double value)
0200 {
0201   lowEnergyLimitQ = value;
0202 }
0203
0204 inline void G4DiffuseElasticV2::SetLowestEnergyLimit(G4double value)
0205 {
0206   lowestEnergyLimit = value;
0207 }
0208
0209
0210 /////////////////////////////////////////////////////////////
0211 //
0212 // Bessel J0 function based on rational approximation from
0213 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
0214
0215 inline G4double G4DiffuseElasticV2::BesselJzero(G4double value)
0216 {
0217   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
0218
0219   modvalue = std::fabs(value);
0220
0221   if ( value < 8.0 && value > -8.0 )
0222   {
0223     value2 = value*value;
0224
0225     fact1  = 57568490574.0 + value2*(-13362590354.0
0226                            + value2*( 651619640.7
0227                            + value2*(-11214424.18
0228                            + value2*( 77392.33017
0229                            + value2*(-184.9052456   ) ) ) ) );
0230
0231     fact2  = 57568490411.0 + value2*( 1029532985.0
0232                            + value2*( 9494680.718
0233                            + value2*(59272.64853
0234                            + value2*(267.8532712
0235                            + value2*1.0               ) ) ) );
0236
0237     bessel = fact1/fact2;
0238   }
0239   else
0240   {
0241     arg    = 8.0/modvalue;
0242
0243     value2 = arg*arg;
0244
0245     shift  = modvalue-0.785398164;
0246
0247     fact1  = 1.0 + value2*(-0.1098628627e-2
0248                  + value2*(0.2734510407e-4
0249                  + value2*(-0.2073370639e-5
0250                  + value2*0.2093887211e-6    ) ) );
0251
0252     fact2  = -0.1562499995e-1 + value2*(0.1430488765e-3
0253                               + value2*(-0.6911147651e-5
0254                               + value2*(0.7621095161e-6
0255                               - value2*0.934945152e-7    ) ) );
0256
0257     bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
0258   }
0259   return bessel;
0260 }
0261
0262 /////////////////////////////////////////////////////////////
0263 //
0264 // Bessel J1 function based on rational approximation from
0265 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
0266
0267 inline G4double G4DiffuseElasticV2::BesselJone(G4double value)
0268 {
0269   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
0270
0271   modvalue = std::fabs(value);
0272
0273   if ( modvalue < 8.0 )
0274   {
0275     value2 = value*value;
0276
0277     fact1  = value*(72362614232.0 + value2*(-7895059235.0
0278                                   + value2*( 242396853.1
0279                                   + value2*(-2972611.439
0280                                   + value2*( 15704.48260
0281                                   + value2*(-30.16036606  ) ) ) ) ) );
0282
0283     fact2  = 144725228442.0 + value2*(2300535178.0
0284                             + value2*(18583304.74
0285                             + value2*(99447.43394
0286                             + value2*(376.9991397
0287                             + value2*1.0             ) ) ) );
0288     bessel = fact1/fact2;
0289   }
0290   else
0291   {
0292     arg    = 8.0/modvalue;
0293
0294     value2 = arg*arg;
0295
0296     shift  = modvalue - 2.356194491;
0297
0298     fact1  = 1.0 + value2*( 0.183105e-2
0299                  + value2*(-0.3516396496e-4
0300                  + value2*(0.2457520174e-5
0301                  + value2*(-0.240337019e-6          ) ) ) );
0302
0303     fact2 = 0.04687499995 + value2*(-0.2002690873e-3
0304                           + value2*( 0.8449199096e-5
0305                           + value2*(-0.88228987e-6
0306                           + value2*0.105787412e-6       ) ) );
0307
0308     bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
0309
0310     if (value < 0.0) bessel = -bessel;
0311   }
0312   return bessel;
0313 }
0314
0315 ////////////////////////////////////////////////////////////////////
0316 //
0317 // damp factor in diffraction x/sh(x), x was already *pi
0318
0319 inline G4double G4DiffuseElasticV2::DampFactor(G4double x)
0320 {
0321   G4double df;
0322   G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
0323
0324   // x *= pi;
0325
0326   if( std::fabs(x) < 0.01 )
0327   {
0328     df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
0329   }
0330   else
0331   {
0332     df = x/std::sinh(x);
0333   }
0334   return df;
0335 }
0336
0337
0338 ////////////////////////////////////////////////////////////////////
0339 //
0340 // return J1(x)/x with special case for small x
0341
0342 inline G4double G4DiffuseElasticV2::BesselOneByArg(G4double x)
0343 {
0344   G4double x2, result;
0345
0346   if( std::fabs(x) < 0.01 )
0347   {
0348    x     *= 0.5;
0349    x2     = x*x;
0350    result = 2. - x2 + x2*x2/6.;
0351   }
0352   else
0353   {
0354     result = BesselJone(x)/x;
0355   }
0356   return result;
0357 }
0358
0359
0360 ////////////////////////////////////////////////////////////////////
0361 //
0362 // return Zommerfeld parameter for Coulomb scattering
0363
0364 inline  G4double G4DiffuseElasticV2::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
0365 {
0366   fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
0367
0368   return fZommerfeld;
0369 }
0370
0371 ////////////////////////////////////////////////////////////////////
0372 //
0373 // return Wentzel correction for Coulomb scattering
0374
0375 inline  G4double G4DiffuseElasticV2::CalculateAm( G4double momentum, G4double n, G4double Z)
0376 {
0377   G4double k   = momentum/CLHEP::hbarc;
0378   G4double ch  = 1.13 + 3.76*n*n;
0379   G4double zn  = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
0380   G4double zn2 = zn*zn;
0381   fAm          = ch/zn2;
0382
0383   return fAm;
0384 }
0385
0386 ////////////////////////////////////////////////////////////////////
0387 //
0388 // calculate nuclear radius for different atomic weights using different approximations
0389
0390 inline  G4double G4DiffuseElasticV2::CalculateNuclearRad( G4double A)
0391 {
0392   G4double R, r0, a11, a12, a13, a2, a3;
0393
0394   a11 = 1.26;  // 1.08, 1.16
0395   a12 = 1.;  // 1.08, 1.16
0396   a13 = 1.12;  // 1.08, 1.16
0397   a2 = 1.1;
0398   a3 = 1.;
0399
0400   // Special rms radii for light nucleii
0401
0402   if (A < 50.)
0403     {
0404       if     (std::abs(A-1.) < 0.5)                         return 0.89*CLHEP::fermi; // p
0405       else if(std::abs(A-2.) < 0.5)                         return 2.13*CLHEP::fermi; // d
0406       else if(  // std::abs(Z-1.) < 0.5 &&
0407           std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi; // t
0408
0409       // else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96CLHEP::fermi; // He3
0410       else if( // std::abs(Z-2.) < 0.5 &&
0411           std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi; // He4
0412
0413       else if(  // std::abs(Z-3.) < 0.5
0414           std::abs(A-7.) < 0.5   )                         return 2.40*CLHEP::fermi; // Li7
0415       else if(  // std::abs(Z-4.) < 0.5
0416           std::abs(A-9.) < 0.5)                         return 2.51*CLHEP::fermi; // Be9
0417
0418       else if( 10.  < A && A <= 16. ) r0  = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;   // 1.08CLHEP::fermi;
0419       else if( 15.  < A && A <= 20. ) r0  = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
0420       else if( 20.  < A && A <= 30. ) r0  = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
0421       else                            r0  = a2*CLHEP::fermi;
0422
0423       R = r0*G4Pow::GetInstance()->A13(A);
0424     }
0425   else
0426     {
0427       r0 = a3*CLHEP::fermi;
0428
0429       R  = r0*G4Pow::GetInstance()->powA(A, 0.27);
0430     }
0431   fNuclearRadius = R;
0432
0433   return R;
0434 }
0435
0436
0437 #endif