G4DiffuseElastic.hh

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // Author: V. Grichine (Vladimir,Grichine@cern.ch)
00030 //
00031 //
00032 // G4 Model: diffuse optical elastic scattering with 4-momentum balance 
00033 //
00034 // Class Description
00035 // Final state production model for hadron nuclear elastic scattering; 
00036 // Class Description - End
00037 //
00038 //
00039 // 24.05.07 V. Grichine, first implementation for hadron (no Coulomb) elastic scattering
00040 // 04.09.07 V. Grichine, implementation for Coulomb elastic scattering
00041 // 12.06.11 V. Grichine, new interface to G4hadronElastic
00042 
00043 #ifndef G4DiffuseElastic_h
00044 #define G4DiffuseElastic_h 1
00045 
00046 #include <CLHEP/Units/PhysicalConstants.h>
00047 #include "globals.hh"
00048 #include "G4HadronElastic.hh"
00049 #include "G4HadProjectile.hh"
00050 #include "G4Nucleus.hh"
00051 
00052 using namespace std;
00053 
00054 class G4ParticleDefinition;
00055 class G4PhysicsTable;
00056 class G4PhysicsLogVector;
00057 
00058 class G4DiffuseElastic : public G4HadronElastic // G4HadronicInteraction
00059 {
00060 public:
00061 
00062   G4DiffuseElastic();
00063 
00064   // G4DiffuseElastic(const G4ParticleDefinition* aParticle);
00065 
00066 
00067 
00068 
00069 
00070   virtual ~G4DiffuseElastic();
00071 
00072   void Initialise();
00073 
00074   void InitialiseOnFly(G4double Z, G4double A);
00075 
00076   void BuildAngleTable();
00077 
00078  
00079   // G4HadFinalState* ApplyYourself(const G4HadProjectile & aTrack, G4Nucleus & targetNucleus);
00080 
00081   virtual G4double SampleInvariantT(const G4ParticleDefinition* p, 
00082                                     G4double plab,
00083                                     G4int Z, G4int A);
00084 
00085   void SetPlabLowLimit(G4double value);
00086 
00087   void SetHEModelLowLimit(G4double value);
00088 
00089   void SetQModelLowLimit(G4double value);
00090 
00091   void SetLowestEnergyLimit(G4double value);
00092 
00093   void SetRecoilKinEnergyLimit(G4double value);
00094 
00095   G4double SampleT(const G4ParticleDefinition* aParticle, 
00096                          G4double p, G4double A);
00097 
00098   G4double SampleTableT(const G4ParticleDefinition* aParticle, 
00099                          G4double p, G4double Z, G4double A);
00100 
00101   G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
00102 
00103   G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p, 
00104                                      G4double Z, G4double A);
00105 
00106   G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
00107 
00108   G4double SampleThetaLab(const G4HadProjectile* aParticle, 
00109                                 G4double tmass, G4double A);
00110 
00111   G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 
00112                                  G4double theta, 
00113                                  G4double momentum, 
00114                                  G4double A         );
00115 
00116   G4double GetInvElasticXsc( const G4ParticleDefinition* particle, 
00117                                  G4double theta, 
00118                                  G4double momentum, 
00119                                  G4double A, G4double Z );
00120 
00121   G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 
00122                                  G4double theta, 
00123                                  G4double momentum, 
00124                                  G4double A, G4double Z );
00125 
00126   G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle, 
00127                                  G4double tMand, 
00128                                  G4double momentum, 
00129                                  G4double A, G4double Z );
00130 
00131   G4double IntegralElasticProb( const G4ParticleDefinition* particle, 
00132                                  G4double theta, 
00133                                  G4double momentum, 
00134                                  G4double A            );
00135   
00136 
00137   G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
00138                                  G4double theta, 
00139                                  G4double momentum, 
00140                                  G4double Z         );
00141 
00142   G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 
00143                                  G4double tMand, 
00144                                  G4double momentum, 
00145                                  G4double A, G4double Z         );
00146 
00147   G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
00148                                  G4double momentum, G4double Z       );
00149 
00150   G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
00151                                  G4double momentum, G4double Z, 
00152                                  G4double theta1, G4double theta2         );
00153 
00154 
00155   G4double CalculateParticleBeta( const G4ParticleDefinition* particle, 
00156                                         G4double momentum    );
00157 
00158   G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
00159 
00160   G4double CalculateAm( G4double momentum, G4double n, G4double Z);
00161 
00162   G4double CalculateNuclearRad( G4double A);
00163 
00164   G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle, 
00165                                 G4double tmass, G4double thetaCMS);
00166 
00167   G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle, 
00168                                 G4double tmass, G4double thetaLab);
00169 
00170   void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
00171                       G4double Z, G4double A);
00172 
00173 
00174 
00175   G4double BesselJzero(G4double z);
00176   G4double BesselJone(G4double z);
00177   G4double DampFactor(G4double z);
00178   G4double BesselOneByArg(G4double z);
00179 
00180   G4double GetDiffElasticProb(G4double theta);
00181   G4double GetDiffElasticSumProb(G4double theta);
00182   G4double GetDiffElasticSumProbA(G4double alpha);
00183   G4double GetIntegrandFunction(G4double theta);
00184 
00185 
00186   G4double GetNuclearRadius(){return fNuclearRadius;};
00187 
00188 private:
00189 
00190 
00191   G4ParticleDefinition* theProton;
00192   G4ParticleDefinition* theNeutron;
00193   G4ParticleDefinition* theDeuteron;
00194   G4ParticleDefinition* theAlpha;
00195 
00196   const G4ParticleDefinition* thePionPlus;
00197   const G4ParticleDefinition* thePionMinus;
00198 
00199   G4double lowEnergyRecoilLimit;  
00200   G4double lowEnergyLimitHE;  
00201   G4double lowEnergyLimitQ;  
00202   G4double lowestEnergyLimit;  
00203   G4double plabLowLimit;
00204 
00205   G4int fEnergyBin;
00206   G4int fAngleBin;
00207 
00208   G4PhysicsLogVector*           fEnergyVector;
00209   G4PhysicsTable*               fAngleTable;
00210   std::vector<G4PhysicsTable*>  fAngleBank;
00211 
00212   std::vector<G4double> fElementNumberVector;
00213   std::vector<G4String> fElementNameVector;
00214 
00215   const G4ParticleDefinition* fParticle;
00216   G4double fWaveVector;
00217   G4double fAtomicWeight;
00218   G4double fAtomicNumber;
00219   G4double fNuclearRadius;
00220   G4double fBeta;
00221   G4double fZommerfeld;
00222   G4double fAm;
00223   G4bool fAddCoulomb;
00224 
00225 };
00226 
00227 
00228 inline void G4DiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
00229 {
00230   lowEnergyRecoilLimit = value;
00231 }
00232 
00233 inline void G4DiffuseElastic::SetPlabLowLimit(G4double value)
00234 {
00235   plabLowLimit = value;
00236 }
00237 
00238 inline void G4DiffuseElastic::SetHEModelLowLimit(G4double value)
00239 {
00240   lowEnergyLimitHE = value;
00241 }
00242 
00243 inline void G4DiffuseElastic::SetQModelLowLimit(G4double value)
00244 {
00245   lowEnergyLimitQ = value;
00246 }
00247 
00248 inline void G4DiffuseElastic::SetLowestEnergyLimit(G4double value)
00249 {
00250   lowestEnergyLimit = value;
00251 }
00252 
00253 
00255 //
00256 // Bessel J0 function based on rational approximation from 
00257 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 
00258 
00259 inline G4double G4DiffuseElastic::BesselJzero(G4double value)
00260 {
00261   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00262 
00263   modvalue = fabs(value);
00264 
00265   if ( value < 8.0 && value > -8.0 )
00266   {
00267     value2 = value*value;
00268 
00269     fact1  = 57568490574.0 + value2*(-13362590354.0 
00270                            + value2*( 651619640.7 
00271                            + value2*(-11214424.18 
00272                            + value2*( 77392.33017 
00273                            + value2*(-184.9052456   ) ) ) ) );
00274 
00275     fact2  = 57568490411.0 + value2*( 1029532985.0 
00276                            + value2*( 9494680.718
00277                            + value2*(59272.64853
00278                            + value2*(267.8532712 
00279                            + value2*1.0               ) ) ) );
00280 
00281     bessel = fact1/fact2;
00282   } 
00283   else 
00284   {
00285     arg    = 8.0/modvalue;
00286 
00287     value2 = arg*arg;
00288 
00289     shift  = modvalue-0.785398164;
00290 
00291     fact1  = 1.0 + value2*(-0.1098628627e-2 
00292                  + value2*(0.2734510407e-4
00293                  + value2*(-0.2073370639e-5 
00294                  + value2*0.2093887211e-6    ) ) );
00295 
00296     fact2  = -0.1562499995e-1 + value2*(0.1430488765e-3
00297                               + value2*(-0.6911147651e-5 
00298                               + value2*(0.7621095161e-6
00299                               - value2*0.934945152e-7    ) ) );
00300 
00301     bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
00302   }
00303   return bessel;
00304 }
00305 
00307 //
00308 // Bessel J1 function based on rational approximation from 
00309 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 
00310 
00311 inline G4double G4DiffuseElastic::BesselJone(G4double value)
00312 {
00313   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00314 
00315   modvalue = fabs(value);
00316 
00317   if ( modvalue < 8.0 ) 
00318   {
00319     value2 = value*value;
00320 
00321     fact1  = value*(72362614232.0 + value2*(-7895059235.0 
00322                                   + value2*( 242396853.1
00323                                   + value2*(-2972611.439 
00324                                   + value2*( 15704.48260 
00325                                   + value2*(-30.16036606  ) ) ) ) ) );
00326 
00327     fact2  = 144725228442.0 + value2*(2300535178.0 
00328                             + value2*(18583304.74
00329                             + value2*(99447.43394 
00330                             + value2*(376.9991397 
00331                             + value2*1.0             ) ) ) );
00332     bessel = fact1/fact2;
00333   } 
00334   else 
00335   {
00336     arg    = 8.0/modvalue;
00337 
00338     value2 = arg*arg;
00339 
00340     shift  = modvalue - 2.356194491;
00341 
00342     fact1  = 1.0 + value2*( 0.183105e-2 
00343                  + value2*(-0.3516396496e-4
00344                  + value2*(0.2457520174e-5 
00345                  + value2*(-0.240337019e-6          ) ) ) );
00346 
00347     fact2 = 0.04687499995 + value2*(-0.2002690873e-3
00348                           + value2*( 0.8449199096e-5
00349                           + value2*(-0.88228987e-6
00350                           + value2*0.105787412e-6       ) ) );
00351 
00352     bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
00353 
00354     if (value < 0.0) bessel = -bessel;
00355   }
00356   return bessel;
00357 }
00358 
00360 //
00361 // damp factor in diffraction x/sh(x), x was already *pi
00362 
00363 inline G4double G4DiffuseElastic::DampFactor(G4double x)
00364 {
00365   G4double df;
00366   G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
00367 
00368   // x *= pi;
00369 
00370   if( std::fabs(x) < 0.01 )
00371   { 
00372     df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
00373   }
00374   else
00375   {
00376     df = x/std::sinh(x); 
00377   }
00378   return df;
00379 }
00380 
00381 
00383 //
00384 // return J1(x)/x with special case for small x
00385 
00386 inline G4double G4DiffuseElastic::BesselOneByArg(G4double x)
00387 {
00388   G4double x2, result;
00389   
00390   if( std::fabs(x) < 0.01 )
00391   { 
00392    x     *= 0.5;
00393    x2     = x*x;
00394    result = 2. - x2 + x2*x2/6.;
00395   }
00396   else
00397   {
00398     result = BesselJone(x)/x; 
00399   }
00400   return result;
00401 }
00402 
00404 //
00405 // return particle beta
00406 
00407 inline  G4double G4DiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle, 
00408                                         G4double momentum    )
00409 {
00410   G4double mass = particle->GetPDGMass();
00411   G4double a    = momentum/mass;
00412   fBeta         = a/std::sqrt(1+a*a);
00413 
00414   return fBeta; 
00415 }
00416 
00418 //
00419 // return Zommerfeld parameter for Coulomb scattering
00420 
00421 inline  G4double G4DiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
00422 {
00423   fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
00424 
00425   return fZommerfeld; 
00426 }
00427 
00429 //
00430 // return Wentzel correction for Coulomb scattering
00431 
00432 inline  G4double G4DiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
00433 {
00434   G4double k   = momentum/CLHEP::hbarc;
00435   G4double ch  = 1.13 + 3.76*n*n;
00436   G4double zn  = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
00437   G4double zn2 = zn*zn;
00438   fAm          = ch/zn2;
00439 
00440   return fAm;
00441 }
00442 
00444 //
00445 // calculate nuclear radius for different atomic weights using different approximations
00446 
00447 inline  G4double G4DiffuseElastic::CalculateNuclearRad( G4double A)
00448 {
00449   G4double r0;
00450 
00451   if( A < 50. )
00452   {
00453     if( A > 10. ) r0  = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi;   // 1.08*fermi;
00454     else          r0  = 1.1*CLHEP::fermi;
00455 
00456     fNuclearRadius = r0*std::pow(A, 1./3.);
00457   }
00458   else
00459   {
00460     r0 = 1.7*CLHEP::fermi;   // 1.7*fermi;
00461 
00462     fNuclearRadius = r0*std::pow(A, 0.27); // 0.27);
00463   }
00464   return fNuclearRadius;
00465 }
00466 
00468 //
00469 // return Coulomb scattering differential xsc with Wentzel correction  
00470 
00471 inline  G4double G4DiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
00472                                  G4double theta, 
00473                                  G4double momentum, 
00474                                  G4double Z         )
00475 {
00476   G4double sinHalfTheta  = std::sin(0.5*theta);
00477   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00478   G4double beta          = CalculateParticleBeta( particle, momentum);
00479   G4double z             = particle->GetPDGCharge();
00480   G4double n             = CalculateZommerfeld( beta, z, Z );
00481   G4double am            = CalculateAm( momentum, n, Z);
00482   G4double k             = momentum/CLHEP::hbarc;
00483   G4double ch            = 0.5*n/k;
00484   G4double ch2           = ch*ch;
00485   G4double xsc           = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
00486 
00487   return xsc;
00488 }
00489 
00490 
00492 //
00493 // return Coulomb scattering total xsc with Wentzel correction  
00494 
00495 inline  G4double G4DiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
00496                                                              G4double momentum, G4double Z  )
00497 {
00498   G4double beta          = CalculateParticleBeta( particle, momentum);
00499   G4cout<<"beta = "<<beta<<G4endl;
00500   G4double z             = particle->GetPDGCharge();
00501   G4double n             = CalculateZommerfeld( beta, z, Z );
00502   G4cout<<"fZomerfeld = "<<n<<G4endl;
00503   G4double am            = CalculateAm( momentum, n, Z);
00504   G4cout<<"cof Am = "<<am<<G4endl;
00505   G4double k             = momentum/CLHEP::hbarc;
00506   G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00507   G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
00508   G4double ch            = n/k;
00509   G4double ch2           = ch*ch;
00510   G4double xsc           = ch2*CLHEP::pi/(am +am*am);
00511 
00512   return xsc;
00513 }
00514 
00516 //
00517 // return Coulomb scattering xsc with Wentzel correction  integrated between
00518 // theta1 and < theta2
00519 
00520 inline  G4double G4DiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
00521                                  G4double momentum, G4double Z, 
00522                                  G4double theta1, G4double theta2 )
00523 {
00524   G4double c1 = std::cos(theta1);
00525   G4cout<<"c1 = "<<c1<<G4endl;
00526   G4double c2 = std::cos(theta2);
00527   G4cout<<"c2 = "<<c2<<G4endl;
00528   G4double beta          = CalculateParticleBeta( particle, momentum);
00529   // G4cout<<"beta = "<<beta<<G4endl;
00530   G4double z             = particle->GetPDGCharge();
00531   G4double n             = CalculateZommerfeld( beta, z, Z );
00532   // G4cout<<"fZomerfeld = "<<n<<G4endl;
00533   G4double am            = CalculateAm( momentum, n, Z);
00534   // G4cout<<"cof Am = "<<am<<G4endl;
00535   G4double k             = momentum/CLHEP::hbarc;
00536   // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00537   // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
00538   G4double ch            = n/k;
00539   G4double ch2           = ch*ch;
00540   am *= 2.;
00541   G4double xsc           = ch2*CLHEP::twopi*(c1-c2);
00542            xsc          /= (1 - c1 + am)*(1 - c2 + am);
00543 
00544   return xsc;
00545 }
00546 
00547 #endif

Generated on Mon May 27 17:48:01 2013 for Geant4 by  doxygen 1.4.7