G4NuclNuclDiffuseElastic.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 //
00030 // G4 Model: optical elastic scattering with 4-momentum balance 
00031 //
00032 // Class Description
00033 // Final state production model for nucleus-nucleus elastic scattering;
00034 // Coulomb amplitude is not considered as correction 
00035 // (as in G4DiffuseElastic)
00036 // Class Description - End
00037 //
00038 //
00039 // 17.03.09 V. Grichine implementation for Coulomb elastic scattering
00040 
00041 
00042 #ifndef G4NuclNuclDiffuseElastic_h
00043 #define G4NuclNuclDiffuseElastic_h 1
00044  
00045 #include <complex>
00046 #include <CLHEP/Units/PhysicalConstants.h>
00047 #include "globals.hh"
00048 #include "G4Integrator.hh"
00049 #include "G4HadronElastic.hh"
00050 #include "G4HadProjectile.hh"
00051 #include "G4Nucleus.hh"
00052 
00053 using namespace std;
00054 
00055 class G4ParticleDefinition;
00056 class G4PhysicsTable;
00057 class G4PhysicsLogVector;
00058 
00059 class G4NuclNuclDiffuseElastic : public G4HadronElastic   // G4HadronicInteraction
00060 {
00061 public:
00062 
00063   G4NuclNuclDiffuseElastic();
00064 
00065   // G4NuclNuclDiffuseElastic(const G4ParticleDefinition* aParticle);
00066 
00067 
00068 
00069 
00070 
00071   virtual ~G4NuclNuclDiffuseElastic();
00072 
00073   void Initialise();
00074 
00075   void InitialiseOnFly(G4double Z, G4double A);
00076 
00077   void BuildAngleTable();
00078 
00079  
00080   // G4HadFinalState * ApplyYourself(const G4HadProjectile & aTrack, G4Nucleus & targetNucleus);
00081 
00082   virtual G4double SampleInvariantT(const G4ParticleDefinition* p, 
00083                                     G4double plab,
00084                                     G4int Z, G4int A);
00085 
00086   void SetPlabLowLimit(G4double value);
00087 
00088   void SetHEModelLowLimit(G4double value);
00089 
00090   void SetQModelLowLimit(G4double value);
00091 
00092   void SetLowestEnergyLimit(G4double value);
00093 
00094   void SetRecoilKinEnergyLimit(G4double value);
00095 
00096   G4double SampleT(const G4ParticleDefinition* aParticle, 
00097                          G4double p, G4double A);
00098 
00099   G4double SampleTableT(const G4ParticleDefinition* aParticle, 
00100                          G4double p, G4double Z, G4double A);
00101 
00102   G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
00103 
00104   G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p, 
00105                                      G4double Z, G4double A);
00106 
00107   G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
00108 
00109   G4double SampleThetaLab(const G4HadProjectile* aParticle, 
00110                                 G4double tmass, G4double A);
00111 
00112   G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 
00113                                  G4double theta, 
00114                                  G4double momentum, 
00115                                  G4double A         );
00116 
00117   G4double GetInvElasticXsc( const G4ParticleDefinition* particle, 
00118                                  G4double theta, 
00119                                  G4double momentum, 
00120                                  G4double A, G4double Z );
00121 
00122   G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 
00123                                  G4double theta, 
00124                                  G4double momentum, 
00125                                  G4double A, G4double Z );
00126 
00127   G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle, 
00128                                  G4double tMand, 
00129                                  G4double momentum, 
00130                                  G4double A, G4double Z );
00131 
00132   G4double IntegralElasticProb( const G4ParticleDefinition* particle, 
00133                                  G4double theta, 
00134                                  G4double momentum, 
00135                                  G4double A            );
00136   
00137 
00138   G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
00139                                  G4double theta, 
00140                                  G4double momentum, 
00141                                  G4double Z         );
00142 
00143   G4double GetRutherfordXsc(     G4double theta       );
00144 
00145   G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 
00146                                  G4double tMand, 
00147                                  G4double momentum, 
00148                                  G4double A, G4double Z         );
00149 
00150   G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
00151                                  G4double momentum, G4double Z       );
00152 
00153   G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
00154                                  G4double momentum, G4double Z, 
00155                                  G4double theta1, G4double theta2         );
00156 
00157 
00158   G4double CalculateParticleBeta( const G4ParticleDefinition* particle, 
00159                                         G4double momentum    );
00160 
00161   G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
00162 
00163   G4double CalculateAm( G4double momentum, G4double n, G4double Z);
00164 
00165   G4double CalculateNuclearRad( G4double A);
00166 
00167   G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle, 
00168                                 G4double tmass, G4double thetaCMS);
00169 
00170   G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle, 
00171                                 G4double tmass, G4double thetaLab);
00172 
00173   void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
00174                       G4double Z, G4double A);
00175 
00176 
00177 
00178   G4double BesselJzero(G4double z);
00179   G4double BesselJone(G4double z);
00180   G4double DampFactor(G4double z);
00181   G4double BesselOneByArg(G4double z);
00182 
00183   G4double GetDiffElasticProb(G4double theta);
00184   G4double GetDiffElasticSumProb(G4double theta);
00185   G4double GetDiffElasticSumProbA(G4double alpha);
00186   G4double GetIntegrandFunction(G4double theta);
00187 
00188   G4double GetNuclearRadius(){return fNuclearRadius;};
00189 
00190 
00191   // Technical math functions for strong Coulomb contribution
00192 
00193   G4complex GammaLogarithm(G4complex xx);
00194   G4complex GammaLogB2n(G4complex xx);
00195 
00196   G4double  GetErf(G4double x);
00197 
00198   G4double GetCosHaPit2(G4double t){return std::cos(CLHEP::halfpi*t*t);};
00199   G4double GetSinHaPit2(G4double t){return std::sin(CLHEP::halfpi*t*t);};
00200 
00201   G4double  GetCint(G4double x);
00202   G4double  GetSint(G4double x);
00203 
00204 
00205   G4complex GetErfcComp(G4complex z, G4int nMax);
00206   G4complex GetErfcSer(G4complex z, G4int nMax);
00207   G4complex GetErfcInt(G4complex z); // , G4int nMax);
00208 
00209   G4complex GetErfComp(G4complex z, G4int nMax);  // AandS algorithm != Ser, Int
00210   G4complex GetErfSer(G4complex z, G4int nMax);
00211 
00212   G4double GetExpCos(G4double x);
00213   G4double GetExpSin(G4double x);
00214   G4complex GetErfInt(G4complex z); // , G4int nMax);
00215 
00216 
00217   
00218 
00219   G4double GetLegendrePol(G4int n, G4double x);
00220 
00221   G4complex TestErfcComp(G4complex z, G4int nMax);
00222   G4complex TestErfcSer(G4complex z, G4int nMax);
00223   G4complex TestErfcInt(G4complex z); // , G4int nMax);
00224 
00225   G4complex CoulombAmplitude(G4double theta);
00226   G4double  CoulombAmplitudeMod2(G4double theta);
00227 
00228 
00229   void CalculateCoulombPhaseZero();
00230   G4double CalculateCoulombPhase(G4int n);
00231   void CalculateRutherfordAnglePar();
00232 
00233   G4double ProfileNear(G4double theta);
00234   G4double ProfileFar(G4double theta);
00235   G4double Profile(G4double theta);
00236 
00237   G4complex PhaseNear(G4double theta);
00238   G4complex PhaseFar(G4double theta);
00239 
00240   G4complex GammaLess(G4double theta);
00241   G4complex GammaMore(G4double theta);
00242 
00243   G4complex AmplitudeNear(G4double theta);
00244   G4complex AmplitudeFar(G4double theta);
00245 
00246   G4complex Amplitude(G4double theta);
00247   G4double  AmplitudeMod2(G4double theta);
00248 
00249   G4complex AmplitudeSim(G4double theta);
00250   G4double  AmplitudeSimMod2(G4double theta);
00251 
00252   G4double  GetRatioSim(G4double theta);
00253   G4double  GetRatioGen(G4double theta);
00254   
00255   G4double  GetFresnelDiffuseXsc(G4double theta);
00256   G4double  GetFresnelIntegrandXsc(G4double alpha);
00257   
00258 
00259   G4complex AmplitudeGla(G4double theta);
00260   G4double  AmplitudeGlaMod2(G4double theta);
00261 
00262   G4complex AmplitudeGG(G4double theta);
00263   G4double  AmplitudeGGMod2(G4double theta);
00264 
00265   void      InitParameters(const G4ParticleDefinition* theParticle,  
00266                               G4double partMom, G4double Z, G4double A);
00267  
00268   void      InitDynParameters(const G4ParticleDefinition* theParticle,  
00269                               G4double partMom); 
00270 
00271   void      InitParametersGla(const G4DynamicParticle* aParticle,  
00272                               G4double partMom, G4double Z, G4double A);
00273 
00274   G4double GetHadronNucleonXscNS( G4ParticleDefinition* pParticle, 
00275                                                  G4double pTkin, 
00276                                   G4ParticleDefinition* tParticle);
00277 
00278   G4double CalcMandelstamS( const G4double mp , 
00279                                                        const G4double mt , 
00280                                                     const G4double Plab );
00281 
00282   G4double GetProfileLambda(){return fProfileLambda;};
00283 
00284   void SetProfileLambda(G4double pl) {fProfileLambda = pl;};
00285   void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
00286   void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
00287   void SetCofLambda(G4double pa){fCofLambda = pa;};
00288 
00289   void SetCofAlpha(G4double pa){fCofAlpha = pa;};
00290   void SetCofAlphaMax(G4double pa){fCofAlphaMax = pa;};
00291   void SetCofAlphaCoulomb(G4double pa){fCofAlphaCoulomb = pa;};
00292 
00293   void SetCofDelta(G4double pa){fCofDelta = pa;};
00294   void SetCofPhase(G4double pa){fCofPhase = pa;};
00295   void SetCofFar(G4double pa){fCofFar = pa;};
00296   void SetEtaRatio(G4double pa){fEtaRatio = pa;};
00297   void SetMaxL(G4int l){fMaxL = l;};
00298   void SetNuclearRadiusCof(G4double r){fNuclearRadiusCof = r;};
00299 
00300   G4double GetCofAlphaMax(){return fCofAlphaMax;};
00301   G4double GetCofAlphaCoulomb(){return fCofAlphaCoulomb;};
00302 
00303 private:
00304 
00305 
00306   G4ParticleDefinition* theProton;
00307   G4ParticleDefinition* theNeutron;
00308   G4ParticleDefinition* theDeuteron;
00309   G4ParticleDefinition* theAlpha;
00310 
00311   const G4ParticleDefinition* thePionPlus;
00312   const G4ParticleDefinition* thePionMinus;
00313 
00314   G4double lowEnergyRecoilLimit;  
00315   G4double lowEnergyLimitHE;  
00316   G4double lowEnergyLimitQ;  
00317   G4double lowestEnergyLimit;  
00318   G4double plabLowLimit;
00319 
00320   G4int fEnergyBin;
00321   G4int fAngleBin;
00322 
00323   G4PhysicsLogVector*           fEnergyVector;
00324   G4PhysicsTable*               fAngleTable;
00325   std::vector<G4PhysicsTable*>  fAngleBank;
00326 
00327   std::vector<G4double> fElementNumberVector;
00328   std::vector<G4String> fElementNameVector;
00329 
00330   const G4ParticleDefinition* fParticle;
00331 
00332   G4double fWaveVector;
00333   G4double fAtomicWeight;
00334   G4double fAtomicNumber;
00335 
00336   G4double fNuclearRadius1;
00337   G4double fNuclearRadius2;
00338   G4double fNuclearRadius;
00339   G4double fNuclearRadiusSquare;
00340   G4double fNuclearRadiusCof;
00341 
00342   G4double fBeta;
00343   G4double fZommerfeld;
00344   G4double fRutherfordRatio;
00345   G4double fAm;
00346   G4bool   fAddCoulomb;
00347 
00348   G4double fCoulombPhase0;
00349   G4double fHalfRutThetaTg;
00350   G4double fHalfRutThetaTg2;
00351   G4double fRutherfordTheta;
00352 
00353   G4double fProfileLambda;
00354   G4double fProfileDelta;
00355   G4double fProfileAlpha;
00356 
00357   G4double fCofLambda;
00358   G4double fCofAlpha;
00359   G4double fCofDelta;
00360   G4double fCofPhase;
00361   G4double fCofFar;
00362 
00363   G4double fCofAlphaMax;
00364   G4double fCofAlphaCoulomb;
00365 
00366   G4int    fMaxL;
00367   G4double fSumSigma;
00368   G4double fEtaRatio;
00369 
00370   G4double fReZ;
00371 
00372 };
00373 
00374 
00375 inline void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
00376 {
00377   lowEnergyRecoilLimit = value;
00378 }
00379 
00380 inline void G4NuclNuclDiffuseElastic::SetPlabLowLimit(G4double value)
00381 {
00382   plabLowLimit = value;
00383 }
00384 
00385 inline void G4NuclNuclDiffuseElastic::SetHEModelLowLimit(G4double value)
00386 {
00387   lowEnergyLimitHE = value;
00388 }
00389 
00390 inline void G4NuclNuclDiffuseElastic::SetQModelLowLimit(G4double value)
00391 {
00392   lowEnergyLimitQ = value;
00393 }
00394 
00395 inline void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit(G4double value)
00396 {
00397   lowestEnergyLimit = value;
00398 }
00399 
00400 
00402 //
00403 // Bessel J0 function based on rational approximation from 
00404 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 
00405 
00406 inline G4double G4NuclNuclDiffuseElastic::BesselJzero(G4double value)
00407 {
00408   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00409 
00410   modvalue = fabs(value);
00411 
00412   if ( value < 8.0 && value > -8.0 )
00413   {
00414     value2 = value*value;
00415 
00416     fact1  = 57568490574.0 + value2*(-13362590354.0 
00417                            + value2*( 651619640.7 
00418                            + value2*(-11214424.18 
00419                            + value2*( 77392.33017 
00420                            + value2*(-184.9052456   ) ) ) ) );
00421 
00422     fact2  = 57568490411.0 + value2*( 1029532985.0 
00423                            + value2*( 9494680.718
00424                            + value2*(59272.64853
00425                            + value2*(267.8532712 
00426                            + value2*1.0               ) ) ) );
00427 
00428     bessel = fact1/fact2;
00429   } 
00430   else 
00431   {
00432     arg    = 8.0/modvalue;
00433 
00434     value2 = arg*arg;
00435 
00436     shift  = modvalue-0.785398164;
00437 
00438     fact1  = 1.0 + value2*(-0.1098628627e-2 
00439                  + value2*(0.2734510407e-4
00440                  + value2*(-0.2073370639e-5 
00441                  + value2*0.2093887211e-6    ) ) );
00442 
00443     fact2  = -0.1562499995e-1 + value2*(0.1430488765e-3
00444                               + value2*(-0.6911147651e-5 
00445                               + value2*(0.7621095161e-6
00446                               - value2*0.934945152e-7    ) ) );
00447 
00448     bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
00449   }
00450   return bessel;
00451 }
00452 
00454 //
00455 // Bessel J1 function based on rational approximation from 
00456 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 
00457 
00458 inline G4double G4NuclNuclDiffuseElastic::BesselJone(G4double value)
00459 {
00460   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00461 
00462   modvalue = fabs(value);
00463 
00464   if ( modvalue < 8.0 ) 
00465   {
00466     value2 = value*value;
00467 
00468     fact1  = value*(72362614232.0 + value2*(-7895059235.0 
00469                                   + value2*( 242396853.1
00470                                   + value2*(-2972611.439 
00471                                   + value2*( 15704.48260 
00472                                   + value2*(-30.16036606  ) ) ) ) ) );
00473 
00474     fact2  = 144725228442.0 + value2*(2300535178.0 
00475                             + value2*(18583304.74
00476                             + value2*(99447.43394 
00477                             + value2*(376.9991397 
00478                             + value2*1.0             ) ) ) );
00479     bessel = fact1/fact2;
00480   } 
00481   else 
00482   {
00483     arg    = 8.0/modvalue;
00484 
00485     value2 = arg*arg;
00486 
00487     shift  = modvalue - 2.356194491;
00488 
00489     fact1  = 1.0 + value2*( 0.183105e-2 
00490                  + value2*(-0.3516396496e-4
00491                  + value2*(0.2457520174e-5 
00492                  + value2*(-0.240337019e-6          ) ) ) );
00493 
00494     fact2 = 0.04687499995 + value2*(-0.2002690873e-3
00495                           + value2*( 0.8449199096e-5
00496                           + value2*(-0.88228987e-6
00497                           + value2*0.105787412e-6       ) ) );
00498 
00499     bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
00500 
00501     if (value < 0.0) bessel = -bessel;
00502   }
00503   return bessel;
00504 }
00505 
00507 //
00508 // damp factor in diffraction x/sh(x), x was already *pi
00509 
00510 inline G4double G4NuclNuclDiffuseElastic::DampFactor(G4double x)
00511 {
00512   G4double df;
00513   G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
00514 
00515   // x *= pi;
00516 
00517   if( std::fabs(x) < 0.01 )
00518   { 
00519     df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
00520   }
00521   else
00522   {
00523     df = x/std::sinh(x); 
00524   }
00525   return df;
00526 }
00527 
00528 
00530 //
00531 // return J1(x)/x with special case for small x
00532 
00533 inline G4double G4NuclNuclDiffuseElastic::BesselOneByArg(G4double x)
00534 {
00535   G4double x2, result;
00536   
00537   if( std::fabs(x) < 0.01 )
00538   { 
00539    x     *= 0.5;
00540    x2     = x*x;
00541    result = 2. - x2 + x2*x2/6.;
00542   }
00543   else
00544   {
00545     result = BesselJone(x)/x; 
00546   }
00547   return result;
00548 }
00549 
00551 //
00552 // return particle beta
00553 
00554 inline  G4double G4NuclNuclDiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle, 
00555                                         G4double momentum    )
00556 {
00557   G4double mass = particle->GetPDGMass();
00558   G4double a    = momentum/mass;
00559   fBeta         = a/std::sqrt(1+a*a);
00560 
00561   return fBeta; 
00562 }
00563 
00565 //
00566 // return Zommerfeld parameter for Coulomb scattering
00567 
00568 inline  G4double G4NuclNuclDiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
00569 {
00570   fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
00571 
00572   return fZommerfeld; 
00573 }
00574 
00576 //
00577 // return Wentzel correction for Coulomb scattering
00578 
00579 inline  G4double G4NuclNuclDiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
00580 {
00581   G4double k   = momentum/CLHEP::hbarc;
00582   G4double ch  = 1.13 + 3.76*n*n;
00583   G4double zn  = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
00584   G4double zn2 = zn*zn;
00585   fAm          = ch/zn2;
00586 
00587   return fAm;
00588 }
00589 
00591 //
00592 // calculate nuclear radius for different atomic weights using different approximations
00593 
00594 inline  G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad( G4double A)
00595 {
00596   G4double r0 = 1.*CLHEP::fermi, radius;
00597   // r0 *= 1.12;
00598   // r0 *= 1.44;
00599   r0 *= fNuclearRadiusCof;
00600 
00601   /*
00602   if( A < 50. )
00603   {
00604     if( A > 10. ) r0  = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi;   // 1.08*fermi;
00605     else          r0  = 1.1*CLHEP::fermi;
00606 
00607     radius = r0*std::pow(A, 1./3.);
00608   }
00609   else
00610   {
00611     r0 = 1.7*CLHEP::fermi;   // 1.7*fermi;
00612 
00613     radius = r0*std::pow(A, 0.27); // 0.27);
00614   }
00615   */
00616   radius = r0*std::pow(A, 1./3.);
00617 
00618   return radius;
00619 }
00620 
00622 //
00623 // return Coulomb scattering differential xsc with Wentzel correction. Test function  
00624 
00625 inline  G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
00626                                  G4double theta, 
00627                                  G4double momentum, 
00628                                  G4double Z         )
00629 {
00630   G4double sinHalfTheta  = std::sin(0.5*theta);
00631   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00632   G4double beta          = CalculateParticleBeta( particle, momentum);
00633   G4double z             = particle->GetPDGCharge();
00634   G4double n             = CalculateZommerfeld( beta, z, Z );
00635   G4double am            = CalculateAm( momentum, n, Z);
00636   G4double k             = momentum/CLHEP::hbarc;
00637   G4double ch            = 0.5*n/k;
00638   G4double ch2           = ch*ch;
00639   G4double xsc           = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
00640 
00641   return xsc;
00642 }
00643 
00645 //
00646 // return Rutherford scattering differential xsc with Wentzel correction. For Sampling.  
00647 
00648 inline  G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc(   G4double theta  )
00649 {
00650   G4double sinHalfTheta  = std::sin(0.5*theta);
00651   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00652 
00653   G4double ch2           = fRutherfordRatio*fRutherfordRatio;
00654 
00655   G4double xsc           = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);
00656 
00657   return xsc;
00658 }
00659 
00660 
00662 //
00663 // return Coulomb scattering total xsc with Wentzel correction  
00664 
00665 inline  G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
00666                                                              G4double momentum, G4double Z  )
00667 {
00668   G4double beta          = CalculateParticleBeta( particle, momentum);
00669   G4cout<<"beta = "<<beta<<G4endl;
00670   G4double z             = particle->GetPDGCharge();
00671   G4double n             = CalculateZommerfeld( beta, z, Z );
00672   G4cout<<"fZomerfeld = "<<n<<G4endl;
00673   G4double am            = CalculateAm( momentum, n, Z);
00674   G4cout<<"cof Am = "<<am<<G4endl;
00675   G4double k             = momentum/CLHEP::hbarc;
00676   G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00677   G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
00678   G4double ch            = n/k;
00679   G4double ch2           = ch*ch;
00680   G4double xsc           = ch2*CLHEP::pi/(am +am*am);
00681 
00682   return xsc;
00683 }
00684 
00686 //
00687 // return Coulomb scattering xsc with Wentzel correction  integrated between
00688 // theta1 and < theta2
00689 
00690 inline  G4double G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
00691                                  G4double momentum, G4double Z, 
00692                                  G4double theta1, G4double theta2 )
00693 {
00694   G4double c1 = std::cos(theta1);
00695   G4cout<<"c1 = "<<c1<<G4endl;
00696   G4double c2 = std::cos(theta2);
00697   G4cout<<"c2 = "<<c2<<G4endl;
00698   G4double beta          = CalculateParticleBeta( particle, momentum);
00699   // G4cout<<"beta = "<<beta<<G4endl;
00700   G4double z             = particle->GetPDGCharge();
00701   G4double n             = CalculateZommerfeld( beta, z, Z );
00702   // G4cout<<"fZomerfeld = "<<n<<G4endl;
00703   G4double am            = CalculateAm( momentum, n, Z);
00704   // G4cout<<"cof Am = "<<am<<G4endl;
00705   G4double k             = momentum/CLHEP::hbarc;
00706   // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00707   // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
00708   G4double ch            = n/k;
00709   G4double ch2           = ch*ch;
00710   am *= 2.;
00711   G4double xsc           = ch2*CLHEP::twopi*(c1-c2);
00712            xsc          /= (1 - c1 + am)*(1 - c2 + am);
00713 
00714   return xsc;
00715 }
00716 
00718 //
00719 // For the calculation of arg Gamma(z) one needs complex extension 
00720 // of ln(Gamma(z))
00721 
00722 inline G4complex G4NuclNuclDiffuseElastic::GammaLogarithm(G4complex zz)
00723 {
00724   static G4double cof[6] = { 76.18009172947146,     -86.50532032941677,
00725                              24.01409824083091,      -1.231739572450155,
00726                               0.1208650973866179e-2, -0.5395239384953e-5  } ;
00727   register G4int j;
00728   G4complex z = zz - 1.0;
00729   G4complex tmp = z + 5.5;
00730   tmp -= (z + 0.5) * std::log(tmp);
00731   G4complex ser = G4complex(1.000000000190015,0.);
00732 
00733   for ( j = 0; j <= 5; j++ )
00734   {
00735     z += 1.0;
00736     ser += cof[j]/z;
00737   }
00738   return -tmp + std::log(2.5066282746310005*ser);
00739 }
00740 
00742 //
00743 // For the calculation of arg Gamma(z) one needs complex extension 
00744 // of ln(Gamma(z)) here is approximate algorithm
00745 
00746 inline G4complex G4NuclNuclDiffuseElastic::GammaLogB2n(G4complex z)
00747 {
00748   G4complex z1 = 12.*z;
00749   G4complex z2 = z*z;
00750   G4complex z3 = z2*z;
00751   G4complex z5 = z2*z3;
00752   G4complex z7 = z2*z5;
00753 
00754   z3 *= 360.;
00755   z5 *= 1260.;
00756   z7 *= 1680.;
00757 
00758   G4complex result  = (z-0.5)*std::log(z)-z+0.5*std::log(CLHEP::twopi);
00759             result += 1./z1 - 1./z3 +1./z5 -1./z7;
00760   return result;
00761 }
00762 
00764 //
00765 //
00766 
00767 inline G4double  G4NuclNuclDiffuseElastic::GetErf(G4double x)
00768 {
00769   G4double t, z, tmp, result;
00770 
00771   z   = std::fabs(x);
00772   t   = 1.0/(1.0+0.5*z);
00773 
00774   tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
00775                 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
00776                 t*(-0.82215223+t*0.17087277)))))))));
00777 
00778   if( x >= 0.) result = 1. - tmp;
00779   else         result = 1. + tmp; 
00780     
00781   return result;
00782 }
00783 
00785 //
00786 //
00787 
00788 inline G4complex G4NuclNuclDiffuseElastic::GetErfcComp(G4complex z, G4int nMax)
00789 {
00790   G4complex erfcz = 1. - GetErfComp( z, nMax);
00791   return erfcz;
00792 }
00793 
00795 //
00796 //
00797 
00798 inline G4complex G4NuclNuclDiffuseElastic::GetErfcSer(G4complex z, G4int nMax)
00799 {
00800   G4complex erfcz = 1. - GetErfSer( z, nMax);
00801   return erfcz;
00802 }
00803 
00805 //
00806 //
00807 
00808 inline G4complex G4NuclNuclDiffuseElastic::GetErfcInt(G4complex z) // , G4int nMax)
00809 {
00810   G4complex erfcz = 1. - GetErfInt( z); // , nMax);
00811   return erfcz;
00812 }
00813 
00814 inline  G4double G4NuclNuclDiffuseElastic::GetLegendrePol(G4int n, G4double theta)
00815 {
00816   G4double legPol, epsilon = 1.e-6;
00817   G4double x = std::cos(theta);
00818 
00819   if     ( n  < 0 ) legPol = 0.;
00820   else if( n == 0 ) legPol = 1.;
00821   else if( n == 1 ) legPol = x;
00822   else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
00823   else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
00824   else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
00825   else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
00826   else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
00827   else           
00828   {
00829     // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
00830 
00831     legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
00832   }
00833   return legPol; 
00834 }
00835 
00836 
00837 
00839 //
00840 //
00841 
00842 inline G4complex G4NuclNuclDiffuseElastic::TestErfcComp(G4complex z, G4int nMax)
00843 {
00844   G4complex miz = G4complex( z.imag(), -z.real() ); 
00845   G4complex erfcz = 1. - GetErfComp( miz, nMax);
00846   G4complex w = std::exp(-z*z)*erfcz;
00847   return w;
00848 }
00849 
00851 //
00852 //
00853 
00854 inline G4complex G4NuclNuclDiffuseElastic::TestErfcSer(G4complex z, G4int nMax)
00855 {
00856   G4complex miz = G4complex( z.imag(), -z.real() ); 
00857   G4complex erfcz = 1. - GetErfSer( miz, nMax);
00858   G4complex w = std::exp(-z*z)*erfcz;
00859   return w;
00860 }
00861 
00863 //
00864 //
00865 
00866 inline G4complex G4NuclNuclDiffuseElastic::TestErfcInt(G4complex z) // , G4int nMax)
00867 {
00868   G4complex miz = G4complex( z.imag(), -z.real() ); 
00869   G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
00870   G4complex w = std::exp(-z*z)*erfcz;
00871   return w;
00872 }
00873 
00875 //
00876 //
00877 
00878 inline G4complex G4NuclNuclDiffuseElastic::GetErfComp(G4complex z, G4int nMax)
00879 {
00880   G4int n;
00881   G4double n2, cofn, shny, chny, fn, gn;
00882 
00883   G4double x = z.real();
00884   G4double y = z.imag();
00885 
00886   G4double outRe = 0., outIm = 0.;
00887 
00888   G4double twox  = 2.*x;
00889   G4double twoxy = twox*y;
00890   G4double twox2 = twox*twox;
00891 
00892   G4double cof1 = std::exp(-x*x)/CLHEP::pi;
00893 
00894   G4double cos2xy = std::cos(twoxy);
00895   G4double sin2xy = std::sin(twoxy);
00896 
00897   G4double twoxcos2xy = twox*cos2xy;
00898   G4double twoxsin2xy = twox*sin2xy;
00899 
00900   for( n = 1; n <= nMax; n++)
00901   {
00902     n2   = n*n;
00903 
00904     cofn = std::exp(-0.5*n2)/(n2+twox2);  // /(n2+0.5*twox2);
00905 
00906     chny = std::cosh(n*y);
00907     shny = std::sinh(n*y);
00908 
00909     fn  = twox - twoxcos2xy*chny + n*sin2xy*shny;
00910     gn  =        twoxsin2xy*chny + n*cos2xy*shny;
00911 
00912     fn *= cofn;
00913     gn *= cofn;
00914 
00915     outRe += fn;
00916     outIm += gn;
00917   }
00918   outRe *= 2*cof1;
00919   outIm *= 2*cof1;
00920 
00921   if(std::abs(x) < 0.0001)
00922   {
00923     outRe += GetErf(x);
00924     outIm += cof1*y;
00925   }
00926   else
00927   {
00928     outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
00929     outIm += cof1*sin2xy/twox;
00930   }
00931   return G4complex(outRe, outIm);
00932 }
00933 
00935 //
00936 //
00937 
00938 inline G4complex G4NuclNuclDiffuseElastic::GetErfSer(G4complex z, G4int nMax)
00939 {
00940   G4int n;
00941   G4double a =1., b = 1., tmp;
00942   G4complex sum = z, d = z;
00943 
00944   for( n = 1; n <= nMax; n++)
00945   {
00946     a *= 2.;
00947     b *= 2.*n +1.;
00948     d *= z*z;
00949 
00950     tmp = a/b;
00951 
00952     sum += tmp*d;
00953   }
00954   sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);
00955 
00956   return sum;
00957 }
00958 
00960 
00961 inline  G4double  G4NuclNuclDiffuseElastic::GetExpCos(G4double x)
00962 {
00963   G4double result;
00964 
00965   result = std::exp(x*x-fReZ*fReZ);
00966   result *= std::cos(2.*x*fReZ);
00967   return result;
00968 }
00969 
00971 
00972 inline  G4double  G4NuclNuclDiffuseElastic::GetExpSin(G4double x)
00973 {
00974   G4double result;
00975 
00976   result = std::exp(x*x-fReZ*fReZ);
00977   result *= std::sin(2.*x*fReZ);
00978   return result;
00979 }
00980 
00981 
00982 
00984 //
00985 //
00986 
00987 inline G4complex G4NuclNuclDiffuseElastic::GetErfInt(G4complex z) // , G4int nMax)
00988 {
00989   G4double outRe, outIm;
00990 
00991   G4double x = z.real();
00992   G4double y = z.imag();
00993   fReZ       = x;
00994 
00995   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
00996 
00997   outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
00998   outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
00999 
01000   outRe *= 2./sqrt(CLHEP::pi);
01001   outIm *= 2./sqrt(CLHEP::pi);
01002 
01003   outRe += GetErf(x);
01004 
01005   return G4complex(outRe, outIm);
01006 }
01007 
01009 //
01010 //
01011 
01012 inline G4double G4NuclNuclDiffuseElastic::GetCint(G4double x)
01013 {
01014   G4double out;
01015 
01016   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
01017 
01018   out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetCosHaPit2, 0., x );
01019 
01020   return out;
01021 }
01022 
01024 //
01025 //
01026 
01027 inline G4double G4NuclNuclDiffuseElastic::GetSint(G4double x)
01028 {
01029   G4double out;
01030 
01031   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
01032 
01033   out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetSinHaPit2, 0., x );
01034 
01035   return out;
01036 }
01037 
01038 
01040 //
01041 //
01042 
01043 inline  G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude(G4double theta)
01044 {
01045   G4complex ca;
01046   
01047   G4double sinHalfTheta  = std::sin(0.5*theta);
01048   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 
01049   sinHalfTheta2         += fAm;
01050 
01051   G4double order         = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
01052   G4complex z            = G4complex(0., order);
01053   ca                     = std::exp(z);
01054 
01055   ca                    *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
01056 
01057   return ca; 
01058 }
01059 
01061 //
01062 //
01063 
01064 inline  G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2(G4double theta)
01065 {
01066   G4complex ca = CoulombAmplitude(theta);
01067   G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();
01068 
01069   return  out;
01070 }
01071 
01073 //
01074 //
01075 
01076 
01077 inline  void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero()
01078 {
01079   G4complex z        = G4complex(1,fZommerfeld); 
01080   // G4complex gammalog = GammaLogarithm(z);
01081   G4complex gammalog = GammaLogB2n(z);
01082   fCoulombPhase0     = gammalog.imag();
01083 }
01084 
01086 //
01087 //
01088 
01089 
01090 inline  G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase(G4int n)
01091 {
01092   G4complex z        = G4complex(1. + n, fZommerfeld); 
01093   // G4complex gammalog = GammaLogarithm(z);
01094   G4complex gammalog = GammaLogB2n(z);
01095   return gammalog.imag();
01096 }
01097 
01098 
01100 //
01101 //
01102 
01103 
01104 inline  void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
01105 {
01106   fHalfRutThetaTg   = fZommerfeld/fProfileLambda;  // (fWaveVector*fNuclearRadius);
01107   fRutherfordTheta  = 2.*std::atan(fHalfRutThetaTg);
01108   fHalfRutThetaTg2  = fHalfRutThetaTg*fHalfRutThetaTg;
01109   // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
01110 
01111 }
01112 
01114 //
01115 //
01116 
01117 inline   G4double G4NuclNuclDiffuseElastic::ProfileNear(G4double theta)
01118 {
01119   G4double dTheta = fRutherfordTheta - theta;
01120   G4double result = 0., argument = 0.;
01121 
01122   if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
01123   else
01124   {
01125     argument = fProfileDelta*dTheta;
01126     result   = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
01127     result  /= std::sinh(CLHEP::pi*argument);
01128     result  -= 1.;
01129     result  /= dTheta;
01130   }
01131   return result;
01132 }
01133 
01135 //
01136 //
01137 
01138 inline   G4double G4NuclNuclDiffuseElastic::ProfileFar(G4double theta)
01139 {
01140   G4double dTheta   = fRutherfordTheta + theta;
01141   G4double argument = fProfileDelta*dTheta;
01142 
01143   G4double result   = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
01144   result           /= std::sinh(CLHEP::pi*argument);
01145   result           /= dTheta;
01146 
01147   return result;
01148 }
01149 
01151 //
01152 //
01153 
01154 inline   G4double G4NuclNuclDiffuseElastic::Profile(G4double theta)
01155 {
01156   G4double dTheta = fRutherfordTheta - theta;
01157   G4double result = 0., argument = 0.;
01158 
01159   if(std::abs(dTheta) < 0.001) result = 1.;
01160   else
01161   {
01162     argument = fProfileDelta*dTheta;
01163     result   = CLHEP::pi*argument;
01164     result  /= std::sinh(CLHEP::pi*argument);
01165   }
01166   return result;
01167 }
01168 
01170 //
01171 //
01172 
01173 inline   G4complex G4NuclNuclDiffuseElastic::PhaseNear(G4double theta)
01174 {
01175   G4double twosigma = 2.*fCoulombPhase0; 
01176   twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
01177   twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
01178   twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;
01179 
01180   twosigma *= fCofPhase;
01181 
01182   G4complex z = G4complex(0., twosigma);
01183 
01184   return std::exp(z);
01185 }
01186 
01188 //
01189 //
01190 
01191 inline   G4complex G4NuclNuclDiffuseElastic::PhaseFar(G4double theta)
01192 {
01193   G4double twosigma = 2.*fCoulombPhase0; 
01194   twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
01195   twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
01196   twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;
01197 
01198   twosigma *= fCofPhase;
01199 
01200   G4complex z = G4complex(0., twosigma);
01201 
01202   return std::exp(z);
01203 }
01204 
01206 //
01207 //
01208 
01209 
01210 inline G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta)
01211 {
01212   G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01213   G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
01214 
01215   G4double u              = std::sqrt(0.5*fProfileLambda/sinThetaR);
01216   G4double kappa          = u/std::sqrt(CLHEP::pi);
01217   G4double dTheta         = theta - fRutherfordTheta;
01218   u                      *= dTheta;
01219   G4double u2             = u*u;
01220   G4double u2m2p3         = u2*2./3.;
01221 
01222   G4complex im            = G4complex(0.,1.);
01223   G4complex order         = G4complex(u,u);
01224   order                  /= std::sqrt(2.);
01225 
01226   G4complex gamma         = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
01227   G4complex a0            = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
01228   G4complex a1            = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
01229   G4complex out           = gamma*(1. - a1*dTheta) - a0;
01230 
01231   return out;
01232 }
01233 
01235 //
01236 //
01237 
01238 inline G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta)
01239 {
01240   G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01241   G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
01242 
01243   G4double u              = std::sqrt(0.5*fProfileLambda/sinThetaR);
01244   G4double kappa          = u/std::sqrt(CLHEP::pi);
01245   G4double dTheta         = theta - fRutherfordTheta;
01246   u                      *= dTheta;
01247   G4double u2             = u*u;
01248   G4double u2m2p3         = u2*2./3.;
01249 
01250   G4complex im            = G4complex(0.,1.);
01251   G4complex order         = G4complex(u,u);
01252   order                  /= std::sqrt(2.);
01253   G4complex gamma         = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
01254   G4complex a0            = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
01255   G4complex a1            = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
01256   G4complex out           = -gamma*(1. - a1*dTheta) - a0;
01257 
01258   return out;
01259 }
01260 
01262 //
01263 //
01264 
01265 inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeNear(G4double theta)
01266 {
01267   G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
01268   G4complex out = G4complex(kappa/fWaveVector,0.);
01269 
01270   out *= PhaseNear(theta);
01271 
01272   if( theta <= fRutherfordTheta )
01273   {
01274     out *= GammaLess(theta) + ProfileNear(theta);
01275     // out *= GammaMore(theta) + ProfileNear(theta);
01276     out += CoulombAmplitude(theta);
01277   }
01278   else
01279   {
01280     out *= GammaMore(theta) + ProfileNear(theta);
01281     // out *= GammaLess(theta) + ProfileNear(theta);
01282   }
01283   return out;
01284 }
01285 
01287 //
01288 //
01289 
01290 inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeFar(G4double theta)
01291 {
01292   G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
01293   G4complex out = G4complex(kappa/fWaveVector,0.);
01294   out *= ProfileFar(theta);
01295   out *= PhaseFar(theta);
01296   return out;
01297 }
01298 
01299 
01301 //
01302 //
01303 
01304 inline  G4complex G4NuclNuclDiffuseElastic::Amplitude(G4double theta)
01305 {
01306  
01307   G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
01308   // G4complex out = AmplitudeNear(theta);
01309   // G4complex out = AmplitudeFar(theta);
01310   return    out;
01311 }
01312 
01314 //
01315 //
01316 
01317 inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeMod2(G4double theta)
01318 {
01319   G4complex out = Amplitude(theta);
01320   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01321   return   mod2;
01322 }
01323 
01325 //
01326 //
01327 
01328 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeSim(G4double theta)
01329 {
01330   G4double sinThetaR  = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01331   G4double dTheta     = 0.5*(theta - fRutherfordTheta);
01332   G4double sindTheta  = std::sin(dTheta);
01333   G4double persqrt2   = std::sqrt(0.5);
01334 
01335   G4complex order     = G4complex(persqrt2,persqrt2);
01336   order              *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
01337   // order              *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
01338 
01339   G4complex out;
01340 
01341   if ( theta <= fRutherfordTheta )
01342   {
01343     out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
01344   }
01345   else
01346   {
01347     out = 0.5*GetErfcInt(order)*ProfileNear(theta);
01348   }
01349 
01350   out *= CoulombAmplitude(theta);
01351 
01352   return out;
01353 }
01354 
01356 //
01357 //
01358 
01359 inline G4double G4NuclNuclDiffuseElastic::GetRatioSim(G4double theta)
01360 {
01361   G4double sinThetaR  = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01362   G4double dTheta     = 0.5*(theta - fRutherfordTheta);
01363   G4double sindTheta  = std::sin(dTheta);
01364 
01365   G4double order      = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
01366   // G4cout<<"order = "<<order<<G4endl;  
01367   G4double cosFresnel = 0.5 - GetCint(order);  
01368   G4double sinFresnel = 0.5 - GetSint(order);  
01369 
01370   G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
01371 
01372   return out;
01373 }
01374 
01376 //
01377 // The ratio el/ruth for Fresnel smooth nucleus profile
01378 
01379 inline G4double G4NuclNuclDiffuseElastic::GetRatioGen(G4double theta)
01380 {
01381   G4double sinThetaR  = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01382   G4double dTheta     = 0.5*(theta - fRutherfordTheta);
01383   G4double sindTheta  = std::sin(dTheta);
01384 
01385   G4double prof       = Profile(theta);
01386   G4double prof2      = prof*prof;
01387   // G4double profmod    = std::abs(prof);
01388   G4double order      = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
01389 
01390   order = std::abs(order); // since sin changes sign!
01391   // G4cout<<"order = "<<order<<G4endl; 
01392  
01393   G4double cosFresnel = GetCint(order);  
01394   G4double sinFresnel = GetSint(order);  
01395 
01396   G4double out;
01397 
01398   if ( theta <= fRutherfordTheta )
01399   {
01400     out  = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2; 
01401     out += ( cosFresnel + sinFresnel - 1. )*prof;
01402   }
01403   else
01404   {
01405     out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
01406   }
01407 
01408   return out;
01409 }
01410 
01412 //
01413 // The xsc for Fresnel smooth nucleus profile
01414 
01415 inline G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc(G4double theta)
01416 {
01417   G4double ratio   = GetRatioGen(theta);
01418   G4double ruthXsc = GetRutherfordXsc(theta);
01419   G4double xsc     = ratio*ruthXsc;
01420   return xsc;
01421 }
01422 
01424 //
01425 // The xsc for Fresnel smooth nucleus profile for integration
01426 
01427 inline G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc(G4double alpha)
01428 {
01429   G4double theta = std::sqrt(alpha);
01430   G4double xsc     = GetFresnelDiffuseXsc(theta);
01431   return xsc;
01432 }
01433 
01434 
01435 
01437 //
01438 //
01439 
01440 inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeSimMod2(G4double theta)
01441 {
01442   G4complex out = AmplitudeSim(theta);
01443   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01444   return   mod2;
01445 }
01446 
01448 //
01449 //
01450 
01451 inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeGla(G4double theta)
01452 {
01453   G4int n;
01454   G4double T12b, b, b2; // cosTheta = std::cos(theta);
01455   G4complex out = G4complex(0.,0.), shiftC, shiftN; 
01456   G4complex im  = G4complex(0.,1.);
01457 
01458   for( n = 0; n < fMaxL; n++)
01459   {
01460     shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
01461     // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
01462     b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
01463     b2 = b*b;
01464     T12b = fSumSigma*std::exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;         
01465     shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
01466     out +=  (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);   
01467   }
01468   out /= 2.*im*fWaveVector;
01469   out += CoulombAmplitude(theta);
01470   return out;
01471 }
01472 
01474 //
01475 //
01476 
01477 inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeGlaMod2(G4double theta)
01478 {
01479   G4complex out = AmplitudeGla(theta);
01480   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01481   return   mod2;
01482 }
01483 
01484 
01486 //
01487 //
01488 
01489 inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeGG(G4double theta)
01490 {
01491   G4int n;
01492   G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
01493   G4double sinThetaH2 = sinThetaH*sinThetaH;
01494   G4complex out = G4complex(0.,0.); 
01495   G4complex im  = G4complex(0.,1.);
01496 
01497   a  = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
01498   b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
01499 
01500   aTemp = a;
01501 
01502   for( n = 1; n < fMaxL; n++)
01503   {    
01504     T12b   = aTemp*std::exp(-b2/n)/n;         
01505     aTemp *= a;  
01506     out   += T12b; 
01507     G4cout<<"out = "<<out<<G4endl;  
01508   }
01509   out *= -4.*im*fWaveVector/CLHEP::pi;
01510   out += CoulombAmplitude(theta);
01511   return out;
01512 }
01513 
01515 //
01516 //
01517 
01518 inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeGGMod2(G4double theta)
01519 {
01520   G4complex out = AmplitudeGG(theta);
01521   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01522   return   mod2;
01523 }
01524 
01525 
01527 //
01528 // Test for given particle and element table of momentum, angle probability.
01529 // For the partMom in CMS. 
01530 
01531 inline void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle,  
01532                                           G4double partMom, G4double Z, G4double A) 
01533 {
01534   fAtomicNumber  = Z;     // atomic number
01535   fAtomicWeight  = A;     // number of nucleons
01536 
01537   fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
01538   G4double A1     = G4double( theParticle->GetBaryonNumber() );   
01539   fNuclearRadius1 = CalculateNuclearRad(A1);
01540   // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
01541   fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
01542 
01543   G4double a = 0.;
01544   G4double z = theParticle->GetPDGCharge();
01545   G4double m1 = theParticle->GetPDGMass();
01546 
01547   fWaveVector = partMom/CLHEP::hbarc;
01548 
01549   G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
01550   G4cout<<"kR = "<<lambda<<G4endl;
01551 
01552   if( z )
01553   {
01554     a           = partMom/m1; // beta*gamma for m1
01555     fBeta       = a/std::sqrt(1+a*a);
01556     fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01557     fRutherfordRatio = fZommerfeld/fWaveVector; 
01558     fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01559   }
01560   G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
01561   fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
01562   G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
01563   fProfileDelta  = fCofDelta*fProfileLambda;
01564   fProfileAlpha  = fCofAlpha*fProfileLambda;
01565 
01566   CalculateCoulombPhaseZero();
01567   CalculateRutherfordAnglePar();
01568 
01569   return;
01570 }
01572 //
01573 // Test for given particle and element table of momentum, angle probability.
01574 // For the partMom in CMS. 
01575 
01576 inline void G4NuclNuclDiffuseElastic::InitDynParameters(const G4ParticleDefinition* theParticle,  
01577                                           G4double partMom) 
01578 {
01579   G4double a = 0.;
01580   G4double z = theParticle->GetPDGCharge();
01581   G4double m1 = theParticle->GetPDGMass();
01582 
01583   fWaveVector = partMom/CLHEP::hbarc;
01584 
01585   G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
01586 
01587   if( z )
01588   {
01589     a           = partMom/m1; // beta*gamma for m1
01590     fBeta       = a/std::sqrt(1+a*a);
01591     fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01592     fRutherfordRatio = fZommerfeld/fWaveVector; 
01593     fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01594   }
01595   fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
01596   fProfileDelta  = fCofDelta*fProfileLambda;
01597   fProfileAlpha  = fCofAlpha*fProfileLambda;
01598 
01599   CalculateCoulombPhaseZero();
01600   CalculateRutherfordAnglePar();
01601 
01602   return;
01603 }
01604 
01605 
01607 //
01608 // Test for given particle and element table of momentum, angle probability.
01609 // For the partMom in CMS. 
01610 
01611 inline void G4NuclNuclDiffuseElastic::InitParametersGla(const G4DynamicParticle* aParticle,  
01612                                           G4double partMom, G4double Z, G4double A) 
01613 {
01614   fAtomicNumber  = Z;     // target atomic number
01615   fAtomicWeight  = A;     // target number of nucleons
01616 
01617   fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
01618   G4double A1     = G4double( aParticle->GetDefinition()->GetBaryonNumber() );   
01619   fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
01620   fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
01621  
01622 
01623   G4double a  = 0., kR12;
01624   G4double z  = aParticle->GetDefinition()->GetPDGCharge();
01625   G4double m1 = aParticle->GetDefinition()->GetPDGMass();
01626 
01627   fWaveVector = partMom/CLHEP::hbarc;
01628 
01629   G4double pN = A1 - z;
01630   if( pN < 0. ) pN = 0.;
01631 
01632   G4double tN = A - Z;
01633   if( tN < 0. ) tN = 0.;
01634 
01635   G4double pTkin = aParticle->GetKineticEnergy();  
01636   pTkin /= A1;
01637 
01638 
01639   fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
01640               (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
01641 
01642   G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
01643   G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl;
01644   kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
01645   G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
01646   fMaxL = (G4int(kR12)+1)*4;
01647   G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
01648 
01649   if( z )
01650   {
01651       a           = partMom/m1; // beta*gamma for m1
01652       fBeta       = a/std::sqrt(1+a*a);
01653       fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01654       fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01655   }
01656 
01657   CalculateCoulombPhaseZero();
01658  
01659 
01660   return;
01661 }
01662 
01663 
01665 //
01666 // Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
01667 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
01668 // projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
01669 
01670 inline G4double 
01671 G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle, 
01672                                                  G4double pTkin, 
01673                                                  G4ParticleDefinition* tParticle)
01674 {
01675   G4double xsection(0), /*Delta,*/ A0, B0;
01676   G4double hpXsc(0);
01677   G4double hnXsc(0);
01678 
01679 
01680   G4double targ_mass     = tParticle->GetPDGMass();
01681   G4double proj_mass     = pParticle->GetPDGMass(); 
01682 
01683   G4double proj_energy   = proj_mass + pTkin; 
01684   G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
01685 
01686   G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
01687 
01688   sMand         /= CLHEP::GeV*CLHEP::GeV;  // in GeV for parametrisation
01689   proj_momentum /= CLHEP::GeV;
01690   proj_energy   /= CLHEP::GeV;
01691   proj_mass     /= CLHEP::GeV;
01692   G4double logS = std::log(sMand);
01693 
01694   // General PDG fit constants
01695 
01696 
01697   // fEtaRatio=Re[f(0)]/Im[f(0)]
01698 
01699   if( proj_momentum >= 1.2 )
01700   {
01701     fEtaRatio  = 0.13*(logS - 5.8579332)*std::pow(sMand,-0.18);
01702   }       
01703   else if( proj_momentum >= 0.6 )
01704   { 
01705     fEtaRatio = -75.5*(std::pow(proj_momentum,0.25)-0.95)/
01706           (std::pow(3*proj_momentum,2.2)+1);     
01707   }
01708   else 
01709   {
01710     fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
01711   }
01712   G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
01713 
01714   // xsc
01715   
01716   if( proj_momentum >= 10. ) // high energy: pp = nn = np
01717     // if( proj_momentum >= 2.)
01718   {
01719     //Delta = 1.;
01720 
01721     //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
01722 
01723     if( proj_momentum >= 10.)
01724     {
01725         B0 = 7.5;
01726         A0 = 100. - B0*std::log(3.0e7);
01727 
01728         xsection = A0 + B0*std::log(proj_energy) - 11
01729                   + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
01730                      0.93827*0.93827,-0.165);        //  mb
01731     }
01732   }
01733   else // low energy pp = nn != np
01734   {
01735       if(pParticle == tParticle) // pp or nn      // nn to be pp
01736       {
01737         if( proj_momentum < 0.73 )
01738         {
01739           hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
01740         }
01741         else if( proj_momentum < 1.05  )
01742         {
01743           hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
01744                          (std::log(proj_momentum/0.73));
01745         }
01746         else  // if( proj_momentum < 10.  )
01747         {
01748           hnXsc = 39.0 + 
01749               75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
01750         }
01751         xsection = hnXsc;
01752       }
01753       else  // pn to be np
01754       {
01755         if( proj_momentum < 0.8 )
01756         {
01757           hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
01758         }      
01759         else if( proj_momentum < 1.4 )
01760         {
01761           hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
01762         }
01763         else    // if( proj_momentum < 10.  )
01764         {
01765           hpXsc = 33.3+
01766               20.8*(std::pow(proj_momentum,2.0)-1.35)/
01767                  (std::pow(proj_momentum,2.50)+0.95);
01768         }
01769         xsection = hpXsc;
01770       }
01771   }
01772   xsection *= CLHEP::millibarn; // parametrised in mb
01773   G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
01774   return xsection;
01775 }
01776 
01778 //
01779 //
01780 
01781 inline G4double G4NuclNuclDiffuseElastic::CalcMandelstamS( const G4double mp , 
01782                                                        const G4double mt , 
01783                                                        const G4double Plab )
01784 {
01785   G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01786   G4double sMand  = mp*mp + mt*mt + 2*Elab*mt ;
01787 
01788   return sMand;
01789 }
01790 
01791 
01792 
01793 #endif

Generated on Mon May 27 17:49:07 2013 for Geant4 by  doxygen 1.4.7