Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/Geant4/G4NuclNuclDiffuseElastic.hh was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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 //
0029 // G4 Model: optical elastic scattering with 4-momentum balance 
0030 //
0031 // Class Description
0032 // Final state production model for nucleus-nucleus elastic scattering;
0033 // Coulomb amplitude is not considered as correction 
0034 // (as in G4DiffuseElastic)
0035 // Class Description - End
0036 //
0037 //
0038 // 17.03.09 V. Grichine implementation for Coulomb elastic scattering
0039 
0040 
0041 #ifndef G4NuclNuclDiffuseElastic_h
0042 #define G4NuclNuclDiffuseElastic_h 1
0043  
0044 #include <complex>
0045 #include <CLHEP/Units/PhysicalConstants.h>
0046 #include "globals.hh"
0047 #include "G4Integrator.hh"
0048 #include "G4HadronElastic.hh"
0049 #include "G4HadProjectile.hh"
0050 #include "G4Nucleus.hh"
0051 
0052 #include "G4Exp.hh"
0053 #include "G4Log.hh"
0054 #include "G4Pow.hh"
0055 
0056 
0057 class G4ParticleDefinition;
0058 class G4PhysicsTable;
0059 class G4PhysicsLogVector;
0060 
0061 class G4NuclNuclDiffuseElastic : public G4HadronElastic   // G4HadronicInteraction
0062 {
0063 public:
0064 
0065   G4NuclNuclDiffuseElastic();
0066 
0067   // G4NuclNuclDiffuseElastic(const G4ParticleDefinition* aParticle);
0068 
0069   virtual ~G4NuclNuclDiffuseElastic();
0070 
0071   void Initialise();
0072 
0073   void InitialiseOnFly(G4double Z, G4double A);
0074 
0075   void BuildAngleTable();
0076 
0077   virtual G4double SampleInvariantT(const G4ParticleDefinition* p, 
0078                     G4double plab,
0079                     G4int Z, G4int A);
0080 
0081   void SetPlabLowLimit(G4double value);
0082 
0083   void SetHEModelLowLimit(G4double value);
0084 
0085   void SetQModelLowLimit(G4double value);
0086 
0087   void SetLowestEnergyLimit(G4double value);
0088 
0089   void SetRecoilKinEnergyLimit(G4double value);
0090 
0091   G4double SampleT(const G4ParticleDefinition* aParticle, 
0092                          G4double p, G4double A);
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 SampleCoulombMuCMS( const G4ParticleDefinition* aParticle, G4double p);
0100 
0101   G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p, 
0102                                      G4double Z, G4double A);
0103 
0104   G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
0105 
0106   G4double SampleThetaLab(const G4HadProjectile* aParticle, 
0107                                 G4double tmass, G4double A);
0108 
0109   G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 
0110                                  G4double theta, 
0111                      G4double momentum, 
0112                  G4double A         );
0113 
0114   G4double GetInvElasticXsc( const G4ParticleDefinition* particle, 
0115                                  G4double theta, 
0116                      G4double momentum, 
0117                  G4double A, G4double Z );
0118 
0119   G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 
0120                                  G4double theta, 
0121                      G4double momentum, 
0122                  G4double A, G4double Z );
0123 
0124   G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle, 
0125                                  G4double tMand, 
0126                      G4double momentum, 
0127                  G4double A, G4double Z );
0128 
0129   G4double IntegralElasticProb( const G4ParticleDefinition* particle, 
0130                                  G4double theta, 
0131                      G4double momentum, 
0132                  G4double A            );
0133   
0134 
0135   G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
0136                                  G4double theta, 
0137                      G4double momentum, 
0138                  G4double Z         );
0139 
0140   G4double GetRutherfordXsc(     G4double theta       );
0141 
0142   G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 
0143                                  G4double tMand, 
0144                      G4double momentum, 
0145                  G4double A, G4double Z         );
0146 
0147   G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
0148                      G4double momentum, G4double Z       );
0149 
0150   G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
0151                      G4double momentum, G4double Z, 
0152                                  G4double theta1, G4double theta2         );
0153 
0154 
0155   G4double CalculateParticleBeta( const G4ParticleDefinition* particle, 
0156                                     G4double momentum    );
0157 
0158   G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
0159 
0160   G4double CalculateAm( G4double momentum, G4double n, G4double Z);
0161 
0162   G4double CalculateNuclearRad( G4double A);
0163 
0164   G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle, 
0165                                 G4double tmass, G4double thetaCMS);
0166 
0167   G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle, 
0168                                 G4double tmass, G4double thetaLab);
0169 
0170   void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
0171               G4double Z, G4double A);
0172 
0173 
0174 
0175   G4double BesselJzero(G4double z);
0176   G4double BesselJone(G4double z);
0177   G4double DampFactor(G4double z);
0178   G4double BesselOneByArg(G4double z);
0179 
0180   G4double GetDiffElasticProb(G4double theta);
0181   G4double GetDiffElasticSumProb(G4double theta);
0182   G4double GetDiffElasticSumProbA(G4double alpha);
0183   G4double GetIntegrandFunction(G4double theta);
0184 
0185   G4double GetNuclearRadius(){return fNuclearRadius;};
0186 
0187 
0188   // Technical math functions for strong Coulomb contribution
0189 
0190   G4complex GammaLogarithm(G4complex xx);
0191   G4complex GammaLogB2n(G4complex xx);
0192 
0193   G4double  GetErf(G4double x);
0194 
0195   G4double GetCosHaPit2(G4double t){return std::cos(CLHEP::halfpi*t*t);};
0196   G4double GetSinHaPit2(G4double t){return std::sin(CLHEP::halfpi*t*t);};
0197 
0198   G4double  GetCint(G4double x);
0199   G4double  GetSint(G4double x);
0200 
0201 
0202   G4complex GetErfcComp(G4complex z, G4int nMax);
0203   G4complex GetErfcSer(G4complex z, G4int nMax);
0204   G4complex GetErfcInt(G4complex z); // , G4int nMax);
0205 
0206   G4complex GetErfComp(G4complex z, G4int nMax);  // AandS algorithm != Ser, Int
0207   G4complex GetErfSer(G4complex z, G4int nMax);
0208 
0209   G4double GetExpCos(G4double x);
0210   G4double GetExpSin(G4double x);
0211   G4complex GetErfInt(G4complex z); // , G4int nMax);
0212 
0213   G4double GetLegendrePol(G4int n, G4double x);
0214 
0215   G4complex TestErfcComp(G4complex z, G4int nMax);
0216   G4complex TestErfcSer(G4complex z, G4int nMax);
0217   G4complex TestErfcInt(G4complex z); // , G4int nMax);
0218 
0219   G4complex CoulombAmplitude(G4double theta);
0220   G4double  CoulombAmplitudeMod2(G4double theta);
0221 
0222   void CalculateCoulombPhaseZero();
0223   G4double CalculateCoulombPhase(G4int n);
0224   void CalculateRutherfordAnglePar();
0225 
0226   G4double ProfileNear(G4double theta);
0227   G4double ProfileFar(G4double theta);
0228   G4double Profile(G4double theta);
0229 
0230   G4complex PhaseNear(G4double theta);
0231   G4complex PhaseFar(G4double theta);
0232 
0233   G4complex GammaLess(G4double theta);
0234   G4complex GammaMore(G4double theta);
0235 
0236   G4complex AmplitudeNear(G4double theta);
0237   G4complex AmplitudeFar(G4double theta);
0238 
0239   G4complex Amplitude(G4double theta);
0240   G4double  AmplitudeMod2(G4double theta);
0241 
0242   G4complex AmplitudeSim(G4double theta);
0243   G4double  AmplitudeSimMod2(G4double theta);
0244 
0245   G4double  GetRatioSim(G4double theta);
0246   G4double  GetRatioGen(G4double theta);
0247   
0248   G4double  GetFresnelDiffuseXsc(G4double theta);
0249   G4double  GetFresnelIntegrandXsc(G4double alpha);
0250   
0251 
0252   G4complex AmplitudeGla(G4double theta);
0253   G4double  AmplitudeGlaMod2(G4double theta);
0254 
0255   G4complex AmplitudeGG(G4double theta);
0256   G4double  AmplitudeGGMod2(G4double theta);
0257 
0258   void      InitParameters(const G4ParticleDefinition* theParticle,  
0259                   G4double partMom, G4double Z, G4double A);
0260  
0261   void      InitDynParameters(const G4ParticleDefinition* theParticle,  
0262                   G4double partMom); 
0263 
0264   void      InitParametersGla(const G4DynamicParticle* aParticle,  
0265                   G4double partMom, G4double Z, G4double A);
0266 
0267   G4double GetHadronNucleonXscNS( G4ParticleDefinition* pParticle, 
0268                                                  G4double pTkin, 
0269                   G4ParticleDefinition* tParticle);
0270 
0271   G4double CalcMandelstamS( const G4double mp , const G4double mt , 
0272                 const G4double Plab );
0273 
0274   G4double GetProfileLambda(){return fProfileLambda;};
0275 
0276   inline void SetProfileLambda(G4double pl) {fProfileLambda = pl;};
0277   inline void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
0278   inline void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
0279   inline void SetCofLambda(G4double pa){fCofLambda = pa;};
0280 
0281   inline void SetCofAlpha(G4double pa){fCofAlpha = pa;};
0282   inline void SetCofAlphaMax(G4double pa){fCofAlphaMax = pa;};
0283   inline void SetCofAlphaCoulomb(G4double pa){fCofAlphaCoulomb = pa;};
0284 
0285   inline void SetCofDelta(G4double pa){fCofDelta = pa;};
0286   inline void SetCofPhase(G4double pa){fCofPhase = pa;};
0287   inline void SetCofFar(G4double pa){fCofFar = pa;};
0288   inline void SetEtaRatio(G4double pa){fEtaRatio = pa;};
0289   inline void SetMaxL(G4int l){fMaxL = l;};
0290   inline void SetNuclearRadiusCof(G4double r){fNuclearRadiusCof = r;};
0291 
0292   inline G4double GetCofAlphaMax(){return fCofAlphaMax;};
0293   inline G4double GetCofAlphaCoulomb(){return fCofAlphaCoulomb;};
0294 
0295 private:
0296 
0297   G4ParticleDefinition* theProton;
0298   G4ParticleDefinition* theNeutron;
0299   G4ParticleDefinition* theDeuteron;
0300   G4ParticleDefinition* theAlpha;
0301 
0302   const G4ParticleDefinition* thePionPlus;
0303   const G4ParticleDefinition* thePionMinus;
0304 
0305   G4double lowEnergyRecoilLimit;  
0306   G4double lowEnergyLimitHE;  
0307   G4double lowEnergyLimitQ;  
0308   G4double lowestEnergyLimit;  
0309   G4double plabLowLimit;
0310 
0311   G4int fEnergyBin;
0312   G4int fAngleBin;
0313 
0314   G4PhysicsLogVector*           fEnergyVector;
0315   G4PhysicsTable*               fAngleTable;
0316   std::vector<G4PhysicsTable*>  fAngleBank;
0317 
0318   std::vector<G4double> fElementNumberVector;
0319   std::vector<G4String> fElementNameVector;
0320 
0321   const G4ParticleDefinition* fParticle;
0322 
0323   G4double fWaveVector;
0324   G4double fAtomicWeight;
0325   G4double fAtomicNumber;
0326 
0327   G4double fNuclearRadius1;
0328   G4double fNuclearRadius2;
0329   G4double fNuclearRadius;
0330   G4double fNuclearRadiusSquare;
0331   G4double fNuclearRadiusCof;
0332 
0333   G4double fBeta;
0334   G4double fZommerfeld;
0335   G4double fRutherfordRatio;
0336   G4double fAm;
0337   G4bool   fAddCoulomb;
0338 
0339   G4double fCoulombPhase0;
0340   G4double fHalfRutThetaTg;
0341   G4double fHalfRutThetaTg2;
0342   G4double fRutherfordTheta;
0343 
0344   G4double fProfileLambda;
0345   G4double fProfileDelta;
0346   G4double fProfileAlpha;
0347 
0348   G4double fCofLambda;
0349   G4double fCofAlpha;
0350   G4double fCofDelta;
0351   G4double fCofPhase;
0352   G4double fCofFar;
0353 
0354   G4double fCofAlphaMax;
0355   G4double fCofAlphaCoulomb;
0356 
0357   G4int    fMaxL;
0358   G4double fSumSigma;
0359   G4double fEtaRatio;
0360 
0361   G4double fReZ;
0362   G4double fCoulombMuC;
0363 
0364 };
0365 
0366 
0367 inline void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
0368 {
0369   lowEnergyRecoilLimit = value;
0370 }
0371 
0372 inline void G4NuclNuclDiffuseElastic::SetPlabLowLimit(G4double value)
0373 {
0374   plabLowLimit = value;
0375 }
0376 
0377 inline void G4NuclNuclDiffuseElastic::SetHEModelLowLimit(G4double value)
0378 {
0379   lowEnergyLimitHE = value;
0380 }
0381 
0382 inline void G4NuclNuclDiffuseElastic::SetQModelLowLimit(G4double value)
0383 {
0384   lowEnergyLimitQ = value;
0385 }
0386 
0387 inline void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit(G4double value)
0388 {
0389   lowestEnergyLimit = value;
0390 }
0391 
0392 ////////////////////////////////////////////////////////////////////
0393 //
0394 // damp factor in diffraction x/sh(x), x was already *pi
0395 
0396 inline G4double G4NuclNuclDiffuseElastic::DampFactor(G4double x)
0397 {
0398   G4double df;
0399   G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
0400 
0401   // x *= pi;
0402 
0403   if( std::fabs(x) < 0.01 )
0404   { 
0405     df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
0406   }
0407   else
0408   {
0409     df = x/std::sinh(x); 
0410   }
0411   return df;
0412 }
0413 
0414 
0415 ////////////////////////////////////////////////////////////////////
0416 //
0417 // return J1(x)/x with special case for small x
0418 
0419 inline G4double G4NuclNuclDiffuseElastic::BesselOneByArg(G4double x)
0420 {
0421   G4double x2, result;
0422   
0423   if( std::fabs(x) < 0.01 )
0424   { 
0425    x     *= 0.5;
0426    x2     = x*x;
0427    result = 2. - x2 + x2*x2/6.;
0428   }
0429   else
0430   {
0431     result = BesselJone(x)/x; 
0432   }
0433   return result;
0434 }
0435 
0436 ////////////////////////////////////////////////////////////////////
0437 //
0438 // return particle beta
0439 
0440 inline  G4double 
0441 G4NuclNuclDiffuseElastic::CalculateParticleBeta(const G4ParticleDefinition* particle, 
0442                         G4double momentum)
0443 {
0444   G4double mass = particle->GetPDGMass();
0445   G4double a    = momentum/mass;
0446   fBeta         = a/std::sqrt(1+a*a);
0447 
0448   return fBeta; 
0449 }
0450 
0451 ////////////////////////////////////////////////////////////////////
0452 //
0453 // return Zommerfeld parameter for Coulomb scattering
0454 
0455 inline  G4double 
0456 G4NuclNuclDiffuseElastic::CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2 )
0457 {
0458   fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
0459 
0460   return fZommerfeld; 
0461 }
0462 
0463 ////////////////////////////////////////////////////////////////////
0464 //
0465 // return Wentzel correction for Coulomb scattering
0466 
0467 inline  G4double 
0468 G4NuclNuclDiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
0469 {
0470   G4double k   = momentum/CLHEP::hbarc;
0471   G4double ch  = 1.13 + 3.76*n*n;
0472   G4double zn  = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
0473   G4double zn2 = zn*zn;
0474   fAm          = ch/zn2;
0475 
0476   return fAm;
0477 }
0478 
0479 ////////////////////////////////////////////////////////////////////
0480 //
0481 // calculate nuclear radius for different atomic weights using different approximations
0482 
0483 inline  G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad( G4double A)
0484 {
0485   G4double r0 = 1.*CLHEP::fermi, radius;
0486   // r0 *= 1.12;
0487   // r0 *= 1.44;
0488   r0 *= fNuclearRadiusCof;
0489 
0490   /*
0491   if( A < 50. )
0492   {
0493     if( A > 10. ) r0  = 1.16*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;   // 1.08*fermi;
0494     else          r0  = 1.1*CLHEP::fermi;
0495 
0496     radius = r0*G4Pow::GetInstance()->A13(A);
0497   }
0498   else
0499   {
0500     r0 = 1.7*CLHEP::fermi;   // 1.7*fermi;
0501 
0502     radius = r0*G4Pow::GetInstance()->powA(A, 0.27); // 0.27);
0503   }
0504   */
0505   radius = r0*G4Pow::GetInstance()->A13(A);
0506 
0507   return radius;
0508 }
0509 
0510 ////////////////////////////////////////////////////////////////////
0511 //
0512 // return Coulomb scattering differential xsc with Wentzel correction. Test function  
0513 
0514 inline  G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
0515                                  G4double theta, 
0516                      G4double momentum, 
0517                  G4double Z         )
0518 {
0519   G4double sinHalfTheta  = std::sin(0.5*theta);
0520   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
0521   G4double beta          = CalculateParticleBeta( particle, momentum);
0522   G4double z             = particle->GetPDGCharge();
0523   G4double n             = CalculateZommerfeld( beta, z, Z );
0524   G4double am            = CalculateAm( momentum, n, Z);
0525   G4double k             = momentum/CLHEP::hbarc;
0526   G4double ch            = 0.5*n/k;
0527   G4double ch2           = ch*ch;
0528   G4double xsc           = ch2/((sinHalfTheta2+am)*(sinHalfTheta2+am));
0529 
0530   return xsc;
0531 }
0532 
0533 ////////////////////////////////////////////////////////////////////
0534 //
0535 // return Rutherford scattering differential xsc with Wentzel correction. For Sampling.  
0536 
0537 inline  G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc(G4double theta)
0538 {
0539   G4double sinHalfTheta  = std::sin(0.5*theta);
0540   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
0541 
0542   G4double ch2           = fRutherfordRatio*fRutherfordRatio;
0543 
0544   G4double xsc           = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);
0545 
0546   return xsc;
0547 }
0548 
0549 ////////////////////////////////////////////////////////////////////
0550 //
0551 // return Coulomb scattering total xsc with Wentzel correction  
0552 
0553 inline  G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
0554                                                  G4double momentum, G4double Z  )
0555 {
0556   G4double beta          = CalculateParticleBeta( particle, momentum);
0557   G4cout<<"beta = "<<beta<<G4endl;
0558   G4double z             = particle->GetPDGCharge();
0559   G4double n             = CalculateZommerfeld( beta, z, Z );
0560   G4cout<<"fZomerfeld = "<<n<<G4endl;
0561   G4double am            = CalculateAm( momentum, n, Z);
0562   G4cout<<"cof Am = "<<am<<G4endl;
0563   G4double k             = momentum/CLHEP::hbarc;
0564   G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
0565   G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
0566   G4double ch            = n/k;
0567   G4double ch2           = ch*ch;
0568   G4double xsc           = ch2*CLHEP::pi/(am +am*am);
0569 
0570   return xsc;
0571 }
0572 
0573 ////////////////////////////////////////////////////////////////////
0574 //
0575 // return Coulomb scattering xsc with Wentzel correction  integrated between
0576 // theta1 and < theta2
0577 
0578 inline  G4double 
0579 G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc(const G4ParticleDefinition* particle,  
0580                         G4double momentum, G4double Z, 
0581                         G4double theta1, G4double theta2 )
0582 {
0583   G4double c1 = std::cos(theta1);
0584   //G4cout<<"c1 = "<<c1<<G4endl;
0585   G4double c2 = std::cos(theta2);
0586   // G4cout<<"c2 = "<<c2<<G4endl;
0587   G4double beta          = CalculateParticleBeta( particle, momentum);
0588   // G4cout<<"beta = "<<beta<<G4endl;
0589   G4double z             = particle->GetPDGCharge();
0590   G4double n             = CalculateZommerfeld( beta, z, Z );
0591   // G4cout<<"fZomerfeld = "<<n<<G4endl;
0592   G4double am            = CalculateAm( momentum, n, Z);
0593   // G4cout<<"cof Am = "<<am<<G4endl;
0594   G4double k             = momentum/CLHEP::hbarc;
0595   // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
0596   // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
0597   G4double ch            = n/k;
0598   G4double ch2           = ch*ch;
0599   am *= 2.;
0600   G4double xsc = ch2*CLHEP::twopi*(c1-c2)/((1 - c1 + am)*(1 - c2 + am));
0601 
0602   return xsc;
0603 }
0604 
0605 ///////////////////////////////////////////////////////////////////
0606 //
0607 // For the calculation of arg Gamma(z) one needs complex extension 
0608 // of ln(Gamma(z)) here is approximate algorithm
0609 
0610 inline G4complex G4NuclNuclDiffuseElastic::GammaLogB2n(G4complex z)
0611 {
0612   G4complex z1 = 12.*z;
0613   G4complex z2 = z*z;
0614   G4complex z3 = z2*z;
0615   G4complex z5 = z2*z3;
0616   G4complex z7 = z2*z5;
0617 
0618   z3 *= 360.;
0619   z5 *= 1260.;
0620   z7 *= 1680.;
0621 
0622   G4complex result  = (z-0.5)*std::log(z)-z+0.5*G4Log(CLHEP::twopi);
0623             result += 1./z1 - 1./z3 +1./z5 -1./z7;
0624   return result;
0625 }
0626 
0627 /////////////////////////////////////////////////////////////////
0628 //
0629 //
0630 
0631 inline G4double  G4NuclNuclDiffuseElastic::GetErf(G4double x)
0632 {
0633   G4double t, z, tmp, result;
0634 
0635   z   = std::fabs(x);
0636   t   = 1.0/(1.0+0.5*z);
0637 
0638   tmp = t*std::exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
0639         t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
0640         t*(-0.82215223+t*0.17087277)))))))));
0641 
0642   if( x >= 0.) result = 1. - tmp;
0643   else         result = 1. + tmp; 
0644     
0645   return result;
0646 }
0647 
0648 /////////////////////////////////////////////////////////////////
0649 //
0650 //
0651 
0652 inline G4complex G4NuclNuclDiffuseElastic::GetErfcComp(G4complex z, G4int nMax)
0653 {
0654   G4complex erfcz = 1. - GetErfComp( z, nMax);
0655   return erfcz;
0656 }
0657 
0658 /////////////////////////////////////////////////////////////////
0659 //
0660 //
0661 
0662 inline G4complex G4NuclNuclDiffuseElastic::GetErfcSer(G4complex z, G4int nMax)
0663 {
0664   G4complex erfcz = 1. - GetErfSer( z, nMax);
0665   return erfcz;
0666 }
0667 
0668 /////////////////////////////////////////////////////////////////
0669 //
0670 //
0671 
0672 inline G4complex G4NuclNuclDiffuseElastic::GetErfcInt(G4complex z) // , G4int nMax)
0673 {
0674   G4complex erfcz = 1. - GetErfInt( z); // , nMax);
0675   return erfcz;
0676 }
0677 
0678 
0679 /////////////////////////////////////////////////////////////////
0680 //
0681 //
0682 
0683 inline G4complex G4NuclNuclDiffuseElastic::TestErfcComp(G4complex z, G4int nMax)
0684 {
0685   G4complex miz = G4complex( z.imag(), -z.real() ); 
0686   G4complex erfcz = 1. - GetErfComp( miz, nMax);
0687   G4complex w = std::exp(-z*z)*erfcz;
0688   return w;
0689 }
0690 
0691 /////////////////////////////////////////////////////////////////
0692 //
0693 //
0694 
0695 inline G4complex G4NuclNuclDiffuseElastic::TestErfcSer(G4complex z, G4int nMax)
0696 {
0697   G4complex miz = G4complex( z.imag(), -z.real() ); 
0698   G4complex erfcz = 1. - GetErfSer( miz, nMax);
0699   G4complex w = std::exp(-z*z)*erfcz;
0700   return w;
0701 }
0702 
0703 /////////////////////////////////////////////////////////////////
0704 //
0705 //
0706 
0707 inline G4complex G4NuclNuclDiffuseElastic::TestErfcInt(G4complex z) // , G4int nMax)
0708 {
0709   G4complex miz = G4complex( z.imag(), -z.real() ); 
0710   G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
0711   G4complex w = std::exp(-z*z)*erfcz;
0712   return w;
0713 }
0714 
0715 /////////////////////////////////////////////////////////////////
0716 //
0717 //
0718 
0719 inline G4complex G4NuclNuclDiffuseElastic::GetErfSer(G4complex z, G4int nMax)
0720 {
0721   G4int n;
0722   G4double a =1., b = 1., tmp;
0723   G4complex sum = z, d = z;
0724 
0725   for( n = 1; n <= nMax; n++)
0726   {
0727     a *= 2.;
0728     b *= 2.*n +1.;
0729     d *= z*z;
0730 
0731     tmp = a/b;
0732 
0733     sum += tmp*d;
0734   }
0735   sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);
0736 
0737   return sum;
0738 }
0739 
0740 /////////////////////////////////////////////////////////////////////
0741 
0742 inline  G4double  G4NuclNuclDiffuseElastic::GetExpCos(G4double x)
0743 {
0744   G4double result;
0745 
0746   result = G4Exp(x*x-fReZ*fReZ);
0747   result *= std::cos(2.*x*fReZ);
0748   return result;
0749 }
0750 
0751 /////////////////////////////////////////////////////////////////////
0752 
0753 inline  G4double  G4NuclNuclDiffuseElastic::GetExpSin(G4double x)
0754 {
0755   G4double result;
0756 
0757   result = G4Exp(x*x-fReZ*fReZ);
0758   result *= std::sin(2.*x*fReZ);
0759   return result;
0760 }
0761 
0762 
0763 /////////////////////////////////////////////////////////////////
0764 //
0765 //
0766 
0767 inline G4double G4NuclNuclDiffuseElastic::GetCint(G4double x)
0768 {
0769   G4double out;
0770 
0771   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
0772 
0773   out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetCosHaPit2, 0., x );
0774 
0775   return out;
0776 }
0777 
0778 /////////////////////////////////////////////////////////////////
0779 //
0780 //
0781 
0782 inline G4double G4NuclNuclDiffuseElastic::GetSint(G4double x)
0783 {
0784   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
0785 
0786   G4double out = 
0787     integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetSinHaPit2, 0., x );
0788 
0789   return out;
0790 }
0791 
0792 
0793 /////////////////////////////////////////////////////////////////
0794 //
0795 //
0796 
0797 inline  G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude(G4double theta)
0798 {
0799   G4complex ca;
0800   
0801   G4double sinHalfTheta  = std::sin(0.5*theta);
0802   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 
0803   sinHalfTheta2         += fAm;
0804 
0805   G4double order         = 2.*fCoulombPhase0 - fZommerfeld*G4Log(sinHalfTheta2);
0806   G4complex z            = G4complex(0., order);
0807   ca                     = std::exp(z);
0808 
0809   ca                    *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
0810 
0811   return ca; 
0812 }
0813 
0814 /////////////////////////////////////////////////////////////////
0815 //
0816 //
0817 
0818 inline  G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2(G4double theta)
0819 {
0820   G4complex ca = CoulombAmplitude(theta);
0821   G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();
0822 
0823   return  out;
0824 }
0825 
0826 /////////////////////////////////////////////////////////////////
0827 //
0828 //
0829 
0830 
0831 inline  void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero()
0832 {
0833   G4complex z        = G4complex(1,fZommerfeld); 
0834   // G4complex gammalog = GammaLogarithm(z);
0835   G4complex gammalog = GammaLogB2n(z);
0836   fCoulombPhase0     = gammalog.imag();
0837 }
0838 
0839 /////////////////////////////////////////////////////////////////
0840 //
0841 //
0842 
0843 
0844 inline  G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase(G4int n)
0845 {
0846   G4complex z        = G4complex(1. + n, fZommerfeld); 
0847   // G4complex gammalog = GammaLogarithm(z);
0848   G4complex gammalog = GammaLogB2n(z);
0849   return gammalog.imag();
0850 }
0851 
0852 
0853 /////////////////////////////////////////////////////////////////
0854 //
0855 //
0856 
0857 
0858 inline  void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
0859 {
0860   fHalfRutThetaTg   = fZommerfeld/fProfileLambda;  // (fWaveVector*fNuclearRadius);
0861   fRutherfordTheta  = 2.*std::atan(fHalfRutThetaTg);
0862   fHalfRutThetaTg2  = fHalfRutThetaTg*fHalfRutThetaTg;
0863   // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
0864 
0865 }
0866 
0867 /////////////////////////////////////////////////////////////////
0868 //
0869 //
0870 
0871 inline   G4double G4NuclNuclDiffuseElastic::ProfileNear(G4double theta)
0872 {
0873   G4double dTheta = fRutherfordTheta - theta;
0874   G4double result = 0., argument = 0.;
0875 
0876   if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
0877   else
0878   {
0879     argument = fProfileDelta*dTheta;
0880     result   = CLHEP::pi*argument*G4Exp(fProfileAlpha*argument);
0881     result  /= std::sinh(CLHEP::pi*argument);
0882     result  -= 1.;
0883     result  /= dTheta;
0884   }
0885   return result;
0886 }
0887 
0888 /////////////////////////////////////////////////////////////////
0889 //
0890 //
0891 
0892 inline   G4double G4NuclNuclDiffuseElastic::ProfileFar(G4double theta)
0893 {
0894   G4double dTheta   = fRutherfordTheta + theta;
0895   G4double argument = fProfileDelta*dTheta;
0896 
0897   G4double result   = CLHEP::pi*argument*G4Exp(fProfileAlpha*argument);
0898   result           /= std::sinh(CLHEP::pi*argument);
0899   result           /= dTheta;
0900 
0901   return result;
0902 }
0903 
0904 /////////////////////////////////////////////////////////////////
0905 //
0906 //
0907 
0908 inline   G4double G4NuclNuclDiffuseElastic::Profile(G4double theta)
0909 {
0910   G4double dTheta = fRutherfordTheta - theta;
0911   G4double result = 0., argument = 0.;
0912 
0913   if(std::abs(dTheta) < 0.001) result = 1.;
0914   else
0915   {
0916     argument = fProfileDelta*dTheta;
0917     result   = CLHEP::pi*argument;
0918     result  /= std::sinh(result);
0919   }
0920   return result;
0921 }
0922 
0923 /////////////////////////////////////////////////////////////////
0924 //
0925 //
0926 
0927 inline   G4complex G4NuclNuclDiffuseElastic::PhaseNear(G4double theta)
0928 {
0929   G4double twosigma = 2.*fCoulombPhase0; 
0930   twosigma -= fZommerfeld*G4Log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
0931   twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
0932   twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;
0933 
0934   twosigma *= fCofPhase;
0935 
0936   G4complex z = G4complex(0., twosigma);
0937 
0938   return std::exp(z);
0939 }
0940 
0941 /////////////////////////////////////////////////////////////////
0942 //
0943 //
0944 
0945 inline   G4complex G4NuclNuclDiffuseElastic::PhaseFar(G4double theta)
0946 {
0947   G4double twosigma = 2.*fCoulombPhase0; 
0948   twosigma -= fZommerfeld*G4Log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
0949   twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
0950   twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;
0951 
0952   twosigma *= fCofPhase;
0953 
0954   G4complex z = G4complex(0., twosigma);
0955 
0956   return std::exp(z);
0957 }
0958 
0959 
0960 /////////////////////////////////////////////////////////////////
0961 //
0962 //
0963 
0964 inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeFar(G4double theta)
0965 {
0966   G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
0967   G4complex out = G4complex(kappa/fWaveVector,0.);
0968   out *= ProfileFar(theta);
0969   out *= PhaseFar(theta);
0970   return out;
0971 }
0972 
0973 
0974 /////////////////////////////////////////////////////////////////
0975 //
0976 //
0977 
0978 inline  G4complex G4NuclNuclDiffuseElastic::Amplitude(G4double theta)
0979 {
0980  
0981   G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
0982   // G4complex out = AmplitudeNear(theta);
0983   // G4complex out = AmplitudeFar(theta);
0984   return    out;
0985 }
0986 
0987 /////////////////////////////////////////////////////////////////
0988 //
0989 //
0990 
0991 inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeMod2(G4double theta)
0992 {
0993   G4complex out = Amplitude(theta);
0994   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
0995   return   mod2;
0996 }
0997 
0998 
0999 /////////////////////////////////////////////////////////////////
1000 //
1001 //
1002 
1003 inline G4double G4NuclNuclDiffuseElastic::GetRatioSim(G4double theta)
1004 {
1005   G4double sinThetaR  = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1006   G4double dTheta     = 0.5*(theta - fRutherfordTheta);
1007   G4double sindTheta  = std::sin(dTheta);
1008 
1009   G4double order      = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1010   // G4cout<<"order = "<<order<<G4endl;  
1011   G4double cosFresnel = 0.5 - GetCint(order);  
1012   G4double sinFresnel = 0.5 - GetSint(order);  
1013 
1014   G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
1015 
1016   return out;
1017 }
1018 
1019 /////////////////////////////////////////////////////////////////
1020 //
1021 // The xsc for Fresnel smooth nucleus profile
1022 
1023 inline G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc(G4double theta)
1024 {
1025   G4double ratio   = GetRatioGen(theta);
1026   G4double ruthXsc = GetRutherfordXsc(theta);
1027   G4double xsc     = ratio*ruthXsc;
1028   return xsc;
1029 }
1030 
1031 /////////////////////////////////////////////////////////////////
1032 //
1033 // The xsc for Fresnel smooth nucleus profile for integration
1034 
1035 inline G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc(G4double alpha)
1036 {
1037   G4double theta = std::sqrt(alpha);
1038   G4double xsc     = GetFresnelDiffuseXsc(theta);
1039   return xsc;
1040 }
1041 
1042 /////////////////////////////////////////////////////////////////
1043 //
1044 //
1045 
1046 inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeSimMod2(G4double theta)
1047 {
1048   G4complex out = AmplitudeSim(theta);
1049   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1050   return   mod2;
1051 }
1052 
1053 
1054 /////////////////////////////////////////////////////////////////
1055 //
1056 //
1057 
1058 inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeGlaMod2(G4double theta)
1059 {
1060   G4complex out = AmplitudeGla(theta);
1061   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1062   return   mod2;
1063 }
1064 
1065 /////////////////////////////////////////////////////////////////
1066 //
1067 //
1068 
1069 inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeGGMod2(G4double theta)
1070 {
1071   G4complex out = AmplitudeGG(theta);
1072   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1073   return   mod2;
1074 }
1075 
1076 ////////////////////////////////////////////////////////////////////////////////////
1077 //
1078 //
1079 
1080 inline G4double G4NuclNuclDiffuseElastic::CalcMandelstamS( const G4double mp , 
1081                                                        const G4double mt , 
1082                                                        const G4double Plab )
1083 {
1084   G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1085   G4double sMand  = mp*mp + mt*mt + 2*Elab*mt ;
1086 
1087   return sMand;
1088 }
1089 
1090 
1091 
1092 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!