G4GlauberGribovCrossSection Class Reference

#include <G4GlauberGribovCrossSection.hh>

Inheritance diagram for G4GlauberGribovCrossSection:

G4VCrossSectionDataSet

Public Member Functions

 G4GlauberGribovCrossSection ()
virtual ~G4GlauberGribovCrossSection ()
virtual G4bool IsIsoApplicable (const G4DynamicParticle *aDP, G4int Z, G4int A, const G4Element *elm=0, const G4Material *mat=0)
virtual G4double GetIsoCrossSection (const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
G4double GetRatioSD (const G4DynamicParticle *, G4int At, G4int Zt)
G4double GetRatioQE (const G4DynamicParticle *, G4int At, G4int Zt)
G4double GetHadronNucleonXsc (const G4DynamicParticle *, const G4Element *)
G4double GetHadronNucleonXsc (const G4DynamicParticle *, G4int At, G4int Zt)
G4double GetHadronNucleonXscPDG (const G4DynamicParticle *, const G4Element *)
G4double GetHadronNucleonXscPDG (const G4DynamicParticle *, G4int At, G4int Zt)
G4double GetHadronNucleonXscNS (const G4DynamicParticle *, const G4Element *)
G4double GetHadronNucleonXscNS (const G4DynamicParticle *, G4int At, G4int Zt)
G4double GetKaonNucleonXscVector (const G4DynamicParticle *, G4int At, G4int Zt)
G4double GetHNinelasticXsc (const G4DynamicParticle *, const G4Element *)
G4double GetHNinelasticXsc (const G4DynamicParticle *, G4int At, G4int Zt)
G4double GetHNinelasticXscVU (const G4DynamicParticle *, G4int At, G4int Zt)
G4double CalculateEcmValue (const G4double, const G4double, const G4double)
G4double CalcMandelstamS (const G4double, const G4double, const G4double)
G4double GetNucleusRadius (const G4DynamicParticle *, const G4Element *)
G4double GetNucleusRadius (G4int At)
virtual void CrossSectionDescription (std::ostream &) const
G4double GetElasticGlauberGribov (const G4DynamicParticle *, G4int Z, G4int A)
G4double GetInelasticGlauberGribov (const G4DynamicParticle *, G4int Z, G4int A)
G4double GetTotalGlauberGribovXsc ()
G4double GetElasticGlauberGribovXsc ()
G4double GetInelasticGlauberGribovXsc ()
G4double GetProductionGlauberGribovXsc ()
G4double GetDiffractionGlauberGribovXsc ()
G4double GetRadiusConst ()
G4double GetParticleBarCorTot (const G4ParticleDefinition *theParticle, G4int Z)
G4double GetParticleBarCorIn (const G4ParticleDefinition *theParticle, G4int Z)
void SetEnergyLowerLimit (G4double E)

Static Public Member Functions

static const char * Default_Name ()

Detailed Description

Definition at line 55 of file G4GlauberGribovCrossSection.hh.


Constructor & Destructor Documentation

G4GlauberGribovCrossSection::G4GlauberGribovCrossSection (  ) 

Definition at line 230 of file G4GlauberGribovCrossSection.cc.

References G4Alpha::Alpha(), G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), G4Deuteron::Deuteron(), G4Gamma::Gamma(), G4He3::He3(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4Neutron::Neutron(), G4OmegaMinus::OmegaMinus(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4Proton::Proton(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), G4SigmaZero::SigmaZero(), G4Triton::Triton(), G4XiMinus::XiMinus(), and G4XiZero::XiZero().

00231  : G4VCrossSectionDataSet(Default_Name()),
00232    fUpperLimit(100000*GeV), fLowerLimit(10.*MeV),// fLowerLimit(3*GeV),
00233    fRadiusConst(1.08*fermi),  // 1.1, 1.3 ?
00234    fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
00235    fDiffractionXsc(0.0), fHadronNucleonXsc(0.0)
00236 {
00237   theGamma    = G4Gamma::Gamma();
00238   theProton   = G4Proton::Proton();
00239   theNeutron  = G4Neutron::Neutron();
00240   theAProton  = G4AntiProton::AntiProton();
00241   theANeutron = G4AntiNeutron::AntiNeutron();
00242   thePiPlus   = G4PionPlus::PionPlus();
00243   thePiMinus  = G4PionMinus::PionMinus();
00244   thePiZero   = G4PionZero::PionZero();
00245   theKPlus    = G4KaonPlus::KaonPlus();
00246   theKMinus   = G4KaonMinus::KaonMinus();
00247   theK0S      = G4KaonZeroShort::KaonZeroShort();
00248   theK0L      = G4KaonZeroLong::KaonZeroLong();
00249   theL        = G4Lambda::Lambda();
00250   theAntiL    = G4AntiLambda::AntiLambda();
00251   theSPlus    = G4SigmaPlus::SigmaPlus();
00252   theASPlus   = G4AntiSigmaPlus::AntiSigmaPlus();
00253   theSMinus   = G4SigmaMinus::SigmaMinus();
00254   theASMinus  = G4AntiSigmaMinus::AntiSigmaMinus();
00255   theS0       = G4SigmaZero::SigmaZero();
00256   theAS0      = G4AntiSigmaZero::AntiSigmaZero();
00257   theXiMinus  = G4XiMinus::XiMinus();
00258   theXi0      = G4XiZero::XiZero();
00259   theAXiMinus = G4AntiXiMinus::AntiXiMinus();
00260   theAXi0     = G4AntiXiZero::AntiXiZero();
00261   theOmega    = G4OmegaMinus::OmegaMinus();
00262   theAOmega   = G4AntiOmegaMinus::AntiOmegaMinus();
00263   theD        = G4Deuteron::Deuteron();
00264   theT        = G4Triton::Triton();
00265   theA        = G4Alpha::Alpha();
00266   theHe3      = G4He3::He3();
00267 
00268   hnXsc = new G4HadronNucleonXsc();
00269 }

G4GlauberGribovCrossSection::~G4GlauberGribovCrossSection (  )  [virtual]

Definition at line 275 of file G4GlauberGribovCrossSection.cc.

00276 {
00277   if (hnXsc) delete hnXsc;
00278 }


Member Function Documentation

G4double G4GlauberGribovCrossSection::CalcMandelstamS ( const   G4double,
const   G4double,
const   G4double 
)

Definition at line 1507 of file G4GlauberGribovCrossSection.cc.

Referenced by GetHadronNucleonXsc(), GetHadronNucleonXscNS(), and GetHadronNucleonXscPDG().

01510 {
01511   G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01512   G4double sMand  = mp*mp + mt*mt + 2*Elab*mt ;
01513 
01514   return sMand;
01515 }

G4double G4GlauberGribovCrossSection::CalculateEcmValue ( const   G4double,
const   G4double,
const   G4double 
)

Definition at line 1491 of file G4GlauberGribovCrossSection.cc.

01494 {
01495   G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01496   G4double Ecm  = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
01497   // G4double Pcm  = Plab * mt / Ecm;
01498   // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
01499 
01500   return Ecm ; // KEcm;
01501 }

void G4GlauberGribovCrossSection::CrossSectionDescription ( std::ostream &   )  const [virtual]

Reimplemented from G4VCrossSectionDataSet.

Definition at line 1521 of file G4GlauberGribovCrossSection.cc.

01522 {
01523   outFile << "G4GlauberGribovCrossSection calculates total, inelastic and\n"
01524           << "elastic cross sections for hadron-nucleus cross sections using\n"
01525           << "the Glauber model with Gribov corrections.  It is valid for all\n"
01526           << "targets except hydrogen, and for incident p, pbar, n, sigma-,\n"
01527           << "pi+, pi-, K+, K- and gammas with energies above 3 GeV.  This is\n"
01528           << "a cross section component which is to be used to build a cross\n"
01529           << "data set.\n";
01530 }

static const char* G4GlauberGribovCrossSection::Default_Name (  )  [inline, static]

Definition at line 62 of file G4GlauberGribovCrossSection.hh.

Referenced by G4BGGNucleonInelasticXS::BuildPhysicsTable(), and G4BGGNucleonElasticXS::BuildPhysicsTable().

00062 {return "Glauber-Gribov";}

G4double G4GlauberGribovCrossSection::GetDiffractionGlauberGribovXsc (  )  [inline]

Definition at line 107 of file G4GlauberGribovCrossSection.hh.

00107 { return fDiffractionXsc; }; 

G4double G4GlauberGribovCrossSection::GetElasticGlauberGribov ( const G4DynamicParticle ,
G4int  Z,
G4int  A 
) [inline]

Definition at line 177 of file G4GlauberGribovCrossSection.hh.

References GetIsoCrossSection().

Referenced by G4BGGPionElasticXS::BuildPhysicsTable(), G4BGGNucleonElasticXS::BuildPhysicsTable(), G4BGGPionElasticXS::GetElementCrossSection(), and G4BGGNucleonElasticXS::GetElementCrossSection().

00179 {
00180   GetIsoCrossSection(dp, Z, A);
00181   return fElasticXsc;
00182 }

G4double G4GlauberGribovCrossSection::GetElasticGlauberGribovXsc (  )  [inline]

Definition at line 104 of file G4GlauberGribovCrossSection.hh.

Referenced by G4NeutronElasticXS::GetElementCrossSection().

00104 { return fElasticXsc;   }; 

G4double G4GlauberGribovCrossSection::GetHadronNucleonXsc ( const G4DynamicParticle ,
G4int  At,
G4int  Zt 
)

Definition at line 539 of file G4GlauberGribovCrossSection.cc.

References CalcMandelstamS(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMass(), and G4DynamicParticle::GetMomentum().

00541 {
00542   G4double xsection;
00543 
00544   //G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);
00545 
00546   G4double targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
00547 
00548   G4double proj_mass     = aParticle->GetMass();
00549   G4double proj_momentum = aParticle->GetMomentum().mag();
00550   G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00551 
00552   sMand /= GeV*GeV;  // in GeV for parametrisation
00553   proj_momentum /= GeV;
00554 
00555   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00556   
00557   G4double aa = At;
00558 
00559   if(theParticle == theGamma) 
00560   {
00561     xsection = aa*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
00562   } 
00563   else if(theParticle == theNeutron) // as proton ??? 
00564   {
00565     xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00566   } 
00567   else if(theParticle == theProton) 
00568   {
00569     xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00570     // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
00571     // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
00572   } 
00573   else if(theParticle == theAProton) 
00574   {
00575     xsection = aa*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
00576   } 
00577   else if(theParticle == thePiPlus) 
00578   {
00579     xsection = aa*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
00580   } 
00581   else if(theParticle == thePiMinus) 
00582   {
00583     // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
00584     xsection = aa*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
00585   } 
00586   else if(theParticle == theKPlus) 
00587   {
00588     xsection = aa*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
00589   } 
00590   else if(theParticle == theKMinus) 
00591   {
00592     xsection = aa*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
00593   }
00594   else  // as proton ??? 
00595   {
00596     xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00597   } 
00598   xsection *= millibarn;
00599   return xsection;
00600 }

G4double G4GlauberGribovCrossSection::GetHadronNucleonXsc ( const G4DynamicParticle ,
const G4Element  
)

Definition at line 522 of file G4GlauberGribovCrossSection.cc.

References G4lrint(), G4Element::GetN(), and G4Element::GetZ().

00524 {
00525   G4int At = G4lrint(anElement->GetN());  // number of nucleons 
00526   G4int Zt = G4lrint(anElement->GetZ());  // number of protons
00527 
00528   return GetHadronNucleonXsc(aParticle, At, Zt);
00529 }

G4double G4GlauberGribovCrossSection::GetHadronNucleonXscNS ( const G4DynamicParticle ,
G4int  At,
G4int  Zt 
)

Definition at line 761 of file G4GlauberGribovCrossSection.cc.

References CalcMandelstamS(), G4DynamicParticle::GetDefinition(), GetHadronNucleonXscPDG(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), G4ParticleTable::GetParticleTable(), G4DynamicParticle::GetTotalEnergy(), G4InuclParticleNames::nn, and G4InuclParticleNames::s0.

00763 {
00764   G4double xsection(0);
00765   // G4double Delta;   DHW 19 May 2011: variable set but not used
00766   G4double A0, B0;
00767   G4double hpXscv(0);
00768   G4double hnXscv(0);
00769 
00770   G4int Nt = At-Zt;              // number of neutrons
00771   if (Nt < 0) Nt = 0;  
00772 
00773   G4double aa = At;
00774   G4double zz = Zt;
00775   G4double nn = Nt;
00776 
00777   G4double targ_mass = G4ParticleTable::GetParticleTable()->
00778   GetIonTable()->GetIonMass(Zt, At);
00779 
00780   targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
00781 
00782   G4double proj_mass     = aParticle->GetMass();
00783   G4double proj_energy   = aParticle->GetTotalEnergy(); 
00784   G4double proj_momentum = aParticle->GetMomentum().mag();
00785 
00786   G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00787 
00788   sMand         /= GeV*GeV;  // in GeV for parametrisation
00789   proj_momentum /= GeV;
00790   proj_energy   /= GeV;
00791   proj_mass     /= GeV;
00792 
00793   // General PDG fit constants
00794 
00795   G4double s0   = 5.38*5.38; // in Gev^2
00796   G4double eta1 = 0.458;
00797   G4double eta2 = 0.458;
00798   G4double B    = 0.308;
00799 
00800 
00801   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00802   
00803 
00804   if(theParticle == theNeutron) 
00805   {
00806     if( proj_momentum >= 373.)
00807     {
00808       return GetHadronNucleonXscPDG(aParticle,At,Zt);
00809     }
00810     else if( proj_momentum >= 10.)
00811     // if( proj_momentum >= 2.)
00812     {
00813       //  Delta = 1.;  // DHW 19 May 2011: variable set but not used
00814       // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
00815 
00816       if(proj_momentum >= 10.)
00817       {
00818         B0 = 7.5;
00819         A0 = 100. - B0*std::log(3.0e7);
00820 
00821         xsection = A0 + B0*std::log(proj_energy) - 11
00822                   + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
00823                      0.93827*0.93827,-0.165);        //  mb
00824       }
00825       xsection *= zz + nn;
00826     }
00827     else
00828     {
00829       // nn to be pp
00830 
00831       if( proj_momentum < 0.73 )
00832       {
00833         hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
00834       }
00835       else if( proj_momentum < 1.05  )
00836       {
00837        hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
00838                          (std::log(proj_momentum/0.73));
00839       }
00840       else  // if( proj_momentum < 10.  )
00841       {
00842          hnXscv = 39.0+
00843               75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
00844       }
00845       // pn to be np
00846 
00847       if( proj_momentum < 0.8 )
00848       {
00849         hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
00850       }      
00851       else if( proj_momentum < 1.4 )
00852       {
00853         hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
00854       }
00855       else    // if( proj_momentum < 10.  )
00856       {
00857         hpXscv = 33.3+
00858               20.8*(std::pow(proj_momentum,2.0)-1.35)/
00859                  (std::pow(proj_momentum,2.50)+0.95);
00860       }
00861       xsection = hpXscv*zz + hnXscv*nn;
00862     }
00863   } 
00864   else if(theParticle == theProton) 
00865   {
00866     if( proj_momentum >= 373.)
00867     {
00868       return GetHadronNucleonXscPDG(aParticle,At,Zt);
00869     }
00870     else if( proj_momentum >= 10.)
00871     // if( proj_momentum >= 2.)
00872     {
00873       // Delta = 1.;  DHW 19 May 2011: variable set but not used
00874       // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
00875 
00876       if(proj_momentum >= 10.)
00877       {
00878         B0 = 7.5;
00879         A0 = 100. - B0*std::log(3.0e7);
00880 
00881         xsection = A0 + B0*std::log(proj_energy) - 11
00882                   + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
00883                      0.93827*0.93827,-0.165);        //  mb
00884       }
00885       xsection *= zz + nn;
00886     }
00887     else
00888     {
00889       // pp
00890 
00891       if( proj_momentum < 0.73 )
00892       {
00893         hpXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
00894       }
00895       else if( proj_momentum < 1.05  )
00896       {
00897        hpXscv = 23 + 40*(std::log(proj_momentum/0.73))*
00898                          (std::log(proj_momentum/0.73));
00899       }
00900       else    // if( proj_momentum < 10.  )
00901       {
00902          hpXscv = 39.0+
00903               75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
00904       }
00905       // pn to be np
00906 
00907       if( proj_momentum < 0.8 )
00908       {
00909         hnXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
00910       }      
00911       else if( proj_momentum < 1.4 )
00912       {
00913         hnXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
00914       }
00915       else   // if( proj_momentum < 10.  )
00916       {
00917         hnXscv = 33.3+
00918               20.8*(std::pow(proj_momentum,2.0)-1.35)/
00919                  (std::pow(proj_momentum,2.50)+0.95);
00920       }
00921       xsection = hpXscv*zz + hnXscv*nn;
00922       // xsection = hpXscv*(Zt + Nt);
00923       // xsection = hnXscv*(Zt + Nt);
00924     }    
00925     // xsection *= 0.95;
00926   } 
00927   else if( theParticle == theAProton ) 
00928   {
00929     // xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00930     //                       + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
00931 
00932     // xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00933     //                    + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
00934 
00935     G4double logP = std::log(proj_momentum);
00936 
00937     if( proj_momentum <= 1.0 )
00938     {
00939       xsection  = zz*(65.55 + 53.84/(proj_momentum+1.e-6)  );
00940     }
00941     else
00942     {
00943       xsection  = zz*( 41.1 + 77.2*std::pow( proj_momentum, -0.68) 
00944                        + 0.293*logP*logP - 1.82*logP );
00945     }
00946     if ( nn > 0.)  
00947     {
00948       xsection += nn*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP);
00949     }
00950     else // H
00951     {
00952       fInelasticXsc =   38.0 + 38.0*std::pow( proj_momentum, -0.96) 
00953                         - 0.169*logP*logP;
00954       fInelasticXsc *=  millibarn;
00955     }    
00956   } 
00957   else if( theParticle == thePiPlus ) 
00958   {
00959     if(proj_momentum < 0.4)
00960     {
00961       G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
00962       hpXscv      = Ex3+20.0;
00963     }
00964     else if( proj_momentum < 1.15 )
00965     {
00966       G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
00967       hpXscv = Ex4+14.0;
00968     }
00969     else if(proj_momentum < 3.5)
00970     {
00971       G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
00972       G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
00973       hpXscv = Ex1+Ex2+27.5;
00974     }
00975     else //  if(proj_momentum > 3.5) // mb
00976     {
00977       hpXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
00978     }
00979     // pi+n = pi-p??
00980 
00981     if(proj_momentum < 0.37)
00982     {
00983       hnXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
00984     }
00985     else if(proj_momentum<0.65)
00986     {
00987        hnXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
00988     }
00989     else if(proj_momentum<1.3)
00990     {
00991       hnXscv = 36.1+
00992                 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
00993                 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
00994     }
00995     else if(proj_momentum<3.0)
00996     {
00997       hnXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
00998                 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
00999                 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
01000     }
01001     else   // mb
01002     {
01003       hnXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); 
01004     }
01005     xsection = hpXscv*zz + hnXscv*nn;
01006   } 
01007   else if(theParticle == thePiMinus) 
01008   {
01009     // pi-n = pi+p??
01010 
01011     if(proj_momentum < 0.4)
01012     {
01013       G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
01014       hnXscv      = Ex3+20.0;
01015     }
01016     else if(proj_momentum < 1.15)
01017     {
01018       G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
01019       hnXscv = Ex4+14.0;
01020     }
01021     else if(proj_momentum < 3.5)
01022     {
01023       G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
01024       G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
01025       hnXscv = Ex1+Ex2+27.5;
01026     }
01027     else //  if(proj_momentum > 3.5) // mb
01028     {
01029       hnXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
01030     }
01031     // pi-p
01032 
01033     if(proj_momentum < 0.37)
01034     {
01035       hpXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
01036     }
01037     else if(proj_momentum<0.65)
01038     {
01039        hpXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
01040     }
01041     else if(proj_momentum<1.3)
01042     {
01043       hpXscv = 36.1+
01044                 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
01045                 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
01046     }
01047     else if(proj_momentum<3.0)
01048     {
01049       hpXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
01050                 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
01051                 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
01052     }
01053     else   // mb
01054     {
01055       hpXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); 
01056     }
01057     xsection = hpXscv*zz + hnXscv*nn;
01058   } 
01059   else if(theParticle == theKPlus) 
01060   {
01061     xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
01062                           + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
01063 
01064     xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
01065                           + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
01066   } 
01067   else if(theParticle == theKMinus) 
01068   {
01069     xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
01070                           + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
01071 
01072     xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
01073                           + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
01074   }
01075   else if(theParticle == theSMinus) 
01076   {
01077     xsection  = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.) 
01078                           - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
01079   } 
01080   else if(theParticle == theGamma) // modify later on
01081   {
01082     xsection  = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.) 
01083                           + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
01084    
01085   } 
01086   else  // as proton ??? 
01087   {
01088     xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
01089                           + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
01090 
01091     xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
01092                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
01093   } 
01094   xsection *= millibarn; // parametrised in mb
01095   return xsection;
01096 }

G4double G4GlauberGribovCrossSection::GetHadronNucleonXscNS ( const G4DynamicParticle ,
const G4Element  
)

Definition at line 743 of file G4GlauberGribovCrossSection.cc.

References G4lrint(), G4Element::GetN(), and G4Element::GetZ().

Referenced by GetHNinelasticXsc(), GetIsoCrossSection(), GetRatioQE(), and GetRatioSD().

00745 {
00746   G4int At = G4lrint(anElement->GetN());  // number of nucleons 
00747   G4int Zt = G4lrint(anElement->GetZ());  // number of protons
00748 
00749   return GetHadronNucleonXscNS(aParticle, At, Zt);
00750 }

G4double G4GlauberGribovCrossSection::GetHadronNucleonXscPDG ( const G4DynamicParticle ,
G4int  At,
G4int  Zt 
)

Definition at line 628 of file G4GlauberGribovCrossSection.cc.

References CalcMandelstamS(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), G4ParticleTable::GetParticleTable(), G4InuclParticleNames::nn, and G4InuclParticleNames::s0.

00630 {
00631   G4double xsection;
00632 
00633   G4int Nt = At-Zt;              // number of neutrons
00634   if (Nt < 0) Nt = 0;
00635   
00636   G4double zz = Zt;
00637   G4double aa = At;
00638   G4double nn = Nt;
00639 
00640   G4double targ_mass = G4ParticleTable::GetParticleTable()->
00641     GetIonTable()->GetIonMass(Zt, At);
00642 
00643   targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
00644 
00645   G4double proj_mass     = aParticle->GetMass(); 
00646   G4double proj_momentum = aParticle->GetMomentum().mag();
00647 
00648   G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00649 
00650   sMand         /= GeV*GeV;  // in GeV for parametrisation
00651 
00652   // General PDG fit constants
00653 
00654   G4double s0   = 5.38*5.38; // in Gev^2
00655   G4double eta1 = 0.458;
00656   G4double eta2 = 0.458;
00657   G4double B    = 0.308;
00658 
00659 
00660   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00661   
00662 
00663   if(theParticle == theNeutron) // proton-neutron fit 
00664   {
00665     xsection = zz*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00666                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00667     xsection  += nn*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00668                       + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
00669   } 
00670   else if(theParticle == theProton) 
00671   {
00672       
00673       xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00674                           + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00675 
00676       xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00677                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00678   } 
00679   else if(theParticle == theAProton) 
00680   {
00681     xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00682                           + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
00683 
00684     xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00685                           + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
00686   } 
00687   else if(theParticle == thePiPlus) 
00688   {
00689     xsection  = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.) 
00690                           + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
00691   } 
00692   else if(theParticle == thePiMinus) 
00693   {
00694     xsection  = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.) 
00695                           + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
00696   } 
00697   else if(theParticle == theKPlus || theParticle == theK0L ) 
00698   {
00699     xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
00700                           + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
00701 
00702     xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
00703                           + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
00704   } 
00705   else if(theParticle == theKMinus || theParticle == theK0S ) 
00706   {
00707     xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
00708                           + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
00709 
00710     xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
00711                           + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
00712   }
00713   else if(theParticle == theSMinus) 
00714   {
00715     xsection  = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.) 
00716                           - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
00717   } 
00718   else if(theParticle == theGamma) // modify later on
00719   {
00720     xsection  = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.) 
00721                           + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
00722    
00723   } 
00724   else  // as proton ??? 
00725   {
00726     xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
00727                           + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00728 
00729     xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
00730                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00731   } 
00732   xsection *= millibarn; // parametrised in mb
00733   return xsection;
00734 }

G4double G4GlauberGribovCrossSection::GetHadronNucleonXscPDG ( const G4DynamicParticle ,
const G4Element  
)

Definition at line 609 of file G4GlauberGribovCrossSection.cc.

References G4lrint(), G4Element::GetN(), and G4Element::GetZ().

Referenced by GetHadronNucleonXscNS(), and GetKaonNucleonXscVector().

00611 {
00612   G4int At = G4lrint(anElement->GetN());  // number of nucleons 
00613   G4int Zt = G4lrint(anElement->GetZ());  // number of protons
00614 
00615   return GetHadronNucleonXscPDG(aParticle, At, Zt);
00616 }

G4double G4GlauberGribovCrossSection::GetHNinelasticXsc ( const G4DynamicParticle ,
G4int  At,
G4int  Zt 
)

Definition at line 1151 of file G4GlauberGribovCrossSection.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), and GetHNinelasticXscVU().

01153 {
01154   G4ParticleDefinition* hadron = aParticle->GetDefinition();
01155   G4double sumInelastic;
01156   G4int Nt = At - Zt;
01157   if(Nt < 0) Nt = 0;
01158   
01159   if( hadron == theKPlus )
01160   {
01161     sumInelastic =  GetHNinelasticXscVU(aParticle, At, Zt);
01162   }
01163   else
01164   {
01165     //sumInelastic  = Zt*GetHadronNucleonXscMK(aParticle, theProton);
01166     // sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);    
01167     sumInelastic  = G4double(Zt)*GetHadronNucleonXscNS(aParticle, 1, 1);
01168     sumInelastic += G4double(Nt)*GetHadronNucleonXscNS(aParticle, 1, 0);    
01169   } 
01170   return sumInelastic;
01171 }

G4double G4GlauberGribovCrossSection::GetHNinelasticXsc ( const G4DynamicParticle ,
const G4Element  
)

Definition at line 1137 of file G4GlauberGribovCrossSection.cc.

References G4lrint(), G4Element::GetN(), and G4Element::GetZ().

Referenced by GetIsoCrossSection(), and GetRatioQE().

01139 {
01140   G4int At = G4lrint(anElement->GetN());  // number of nucleons 
01141   G4int Zt = G4lrint(anElement->GetZ());  // number of protons
01142 
01143   return GetHNinelasticXsc(aParticle, At, Zt);
01144 }

G4double G4GlauberGribovCrossSection::GetHNinelasticXscVU ( const G4DynamicParticle ,
G4int  At,
G4int  Zt 
)

Definition at line 1179 of file G4GlauberGribovCrossSection.cc.

References G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMomentum(), G4ParticleDefinition::GetPDGEncoding(), and G4DynamicParticle::GetTotalEnergy().

Referenced by GetHNinelasticXsc().

01181 {
01182   G4int PDGcode    = aParticle->GetDefinition()->GetPDGEncoding();
01183   G4int absPDGcode = std::abs(PDGcode);
01184 
01185   G4double Elab = aParticle->GetTotalEnergy();              
01186                           // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
01187   G4double Plab = aParticle->GetMomentum().mag();            
01188                           // std::sqrt(Elab * Elab - 0.88);
01189 
01190   Elab /= GeV;
01191   Plab /= GeV;
01192 
01193   G4double LogPlab    = std::log( Plab );
01194   G4double sqrLogPlab = LogPlab * LogPlab;
01195 
01196   //G4cout<<"Plab = "<<Plab<<G4endl;
01197 
01198   G4double NumberOfTargetProtons = G4double(Zt); 
01199   G4double NumberOfTargetNucleons = G4double(At);
01200   G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
01201 
01202   if(NumberOfTargetNeutrons < 0.0) NumberOfTargetNeutrons = 0.0;
01203 
01204   G4double Xtotal, Xelastic, Xinelastic;
01205 
01206   if( absPDGcode > 1000 )  //------Projectile is baryon --------
01207   {
01208        G4double XtotPP = 48.0 +  0. *std::pow(Plab, 0.  ) +
01209                          0.522*sqrLogPlab - 4.51*LogPlab;
01210 
01211        G4double XtotPN = 47.3 +  0. *std::pow(Plab, 0.  ) +
01212                          0.513*sqrLogPlab - 4.27*LogPlab;
01213 
01214        G4double XelPP  = 11.9 + 26.9*std::pow(Plab,-1.21) +
01215                          0.169*sqrLogPlab - 1.85*LogPlab;
01216 
01217        G4double XelPN  = 11.9 + 26.9*std::pow(Plab,-1.21) +
01218                          0.169*sqrLogPlab - 1.85*LogPlab;
01219 
01220        Xtotal          = (NumberOfTargetProtons * XtotPP +
01221                           NumberOfTargetNeutrons * XtotPN);
01222 
01223        Xelastic        = (NumberOfTargetProtons * XelPP +
01224                           NumberOfTargetNeutrons * XelPN);
01225   }
01226   else if( PDGcode ==  211 ) //------Projectile is PionPlus -------
01227   {
01228        G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) +
01229                           0.19 *sqrLogPlab - 0.0 *LogPlab;
01230 
01231        G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) +
01232                           0.456*sqrLogPlab - 4.03*LogPlab;
01233 
01234        G4double XelPiP  =  0.0 + 11.4*std::pow(Plab,-0.40) +
01235                            0.079*sqrLogPlab - 0.0 *LogPlab;
01236 
01237        G4double XelPiN  = 1.76 + 11.2*std::pow(Plab,-0.64) +
01238                           0.043*sqrLogPlab - 0.0 *LogPlab;
01239 
01240        Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
01241                             NumberOfTargetNeutrons * XtotPiN  );
01242 
01243        Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
01244                             NumberOfTargetNeutrons * XelPiN   );
01245   }
01246   else if( PDGcode == -211 ) //------Projectile is PionMinus -------
01247   {
01248        G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) +
01249                           0.456*sqrLogPlab - 4.03*LogPlab;
01250 
01251        G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) +
01252                           0.19 *sqrLogPlab - 0.0 *LogPlab;
01253 
01254        G4double XelPiP  = 1.76 + 11.2*std::pow(Plab,-0.64) +
01255                           0.043*sqrLogPlab - 0.0 *LogPlab;
01256 
01257        G4double XelPiN  =  0.0 + 11.4*std::pow(Plab,-0.40) +
01258                            0.079*sqrLogPlab - 0.0 *LogPlab;
01259 
01260        Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
01261                             NumberOfTargetNeutrons * XtotPiN  );
01262 
01263        Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
01264                             NumberOfTargetNeutrons * XelPiN   );
01265   }
01266   else if( PDGcode ==  111 )  //------Projectile is PionZero  -------
01267   {
01268        G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) +
01269                           0.19 *sqrLogPlab - 0.0 *LogPlab +   //Pi+
01270                           33.0 + 14.0 *std::pow(Plab,-1.36) +
01271                           0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
01272 
01273        G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) +
01274                           0.456*sqrLogPlab - 4.03*LogPlab +   //Pi+
01275                           16.4 + 19.3 *std::pow(Plab,-0.42) +
01276                           0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
01277 
01278        G4double XelPiP  =( 0.0 + 11.4*std::pow(Plab,-0.40) +
01279                            0.079*sqrLogPlab - 0.0 *LogPlab +    //Pi+
01280                            1.76 + 11.2*std::pow(Plab,-0.64) +
01281                            0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
01282 
01283        G4double XelPiN  =( 1.76 + 11.2*std::pow(Plab,-0.64) +
01284                            0.043*sqrLogPlab - 0.0 *LogPlab +   //Pi+
01285                            0.0  + 11.4*std::pow(Plab,-0.40) +
01286                            0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
01287 
01288        Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
01289                             NumberOfTargetNeutrons * XtotPiN  );
01290 
01291        Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
01292                             NumberOfTargetNeutrons * XelPiN   );
01293   }
01294   else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
01295   {
01296        G4double XtotKP = 18.1 +  0. *std::pow(Plab, 0.  ) +
01297                          0.26 *sqrLogPlab - 1.0 *LogPlab;
01298        G4double XtotKN = 18.7 +  0. *std::pow(Plab, 0.  ) +
01299                          0.21 *sqrLogPlab - 0.89*LogPlab;
01300 
01301        G4double XelKP  =  5.0 +  8.1*std::pow(Plab,-1.8 ) +
01302                           0.16 *sqrLogPlab - 1.3 *LogPlab;
01303 
01304        G4double XelKN  =  7.3 +  0. *std::pow(Plab,-0.  ) +
01305                           0.29 *sqrLogPlab - 2.4 *LogPlab;
01306 
01307        Xtotal          = ( NumberOfTargetProtons  * XtotKP +
01308                            NumberOfTargetNeutrons * XtotKN  );
01309 
01310        Xelastic        = ( NumberOfTargetProtons  * XelKP  +
01311                            NumberOfTargetNeutrons * XelKN   );
01312   }
01313   else if( PDGcode ==-321 )  //------Projectile is KaonMinus ------
01314   {
01315        G4double XtotKP = 32.1 +  0. *std::pow(Plab, 0.  ) +
01316                          0.66 *sqrLogPlab - 5.6 *LogPlab;
01317        G4double XtotKN = 25.2 +  0. *std::pow(Plab, 0.  ) +
01318                          0.38 *sqrLogPlab - 2.9 *LogPlab;
01319 
01320        G4double XelKP  =  7.3 +  0. *std::pow(Plab,-0.  ) +
01321                           0.29 *sqrLogPlab - 2.4 *LogPlab;
01322 
01323        G4double XelKN  =  5.0 +  8.1*std::pow(Plab,-1.8 ) +
01324                           0.16 *sqrLogPlab - 1.3 *LogPlab;
01325 
01326        Xtotal          = ( NumberOfTargetProtons  * XtotKP +
01327                            NumberOfTargetNeutrons * XtotKN  );
01328 
01329        Xelastic        = ( NumberOfTargetProtons  * XelKP  +
01330                            NumberOfTargetNeutrons * XelKN   );
01331   }
01332   else if( PDGcode == 311 ) //------Projectile is KaonZero ------
01333   {
01334        G4double XtotKP = ( 18.1 +  0. *std::pow(Plab, 0.  ) +
01335                           0.26 *sqrLogPlab - 1.0 *LogPlab +   //K+
01336                           32.1 +  0. *std::pow(Plab, 0.  ) +
01337                           0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
01338 
01339        G4double XtotKN = ( 18.7 +  0. *std::pow(Plab, 0.  ) +
01340                           0.21 *sqrLogPlab - 0.89*LogPlab +   //K+
01341                           25.2 +  0. *std::pow(Plab, 0.  ) +
01342                           0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
01343 
01344        G4double XelKP  = (  5.0 +  8.1*std::pow(Plab,-1.8 )
01345                            + 0.16 *sqrLogPlab - 1.3 *LogPlab +   //K+
01346                            7.3 +  0. *std::pow(Plab,-0.  ) +
01347                            0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
01348 
01349        G4double XelKN  = (  7.3 +  0. *std::pow(Plab,-0.  ) +
01350                            0.29 *sqrLogPlab - 2.4 *LogPlab +   //K+
01351                            5.0 +  8.1*std::pow(Plab,-1.8 ) +
01352                            0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
01353 
01354        Xtotal          = ( NumberOfTargetProtons  * XtotKP +
01355                            NumberOfTargetNeutrons * XtotKN  );
01356 
01357        Xelastic        = ( NumberOfTargetProtons  * XelKP  +
01358                            NumberOfTargetNeutrons * XelKN   );
01359   }
01360   else  //------Projectile is undefined, Nucleon assumed
01361   {
01362        G4double XtotPP = 48.0 +  0. *std::pow(Plab, 0.  ) +
01363                          0.522*sqrLogPlab - 4.51*LogPlab;
01364 
01365        G4double XtotPN = 47.3 +  0. *std::pow(Plab, 0.  ) +
01366                          0.513*sqrLogPlab - 4.27*LogPlab;
01367 
01368        G4double XelPP  = 11.9 + 26.9*std::pow(Plab,-1.21) +
01369                          0.169*sqrLogPlab - 1.85*LogPlab;
01370        G4double XelPN  = 11.9 + 26.9*std::pow(Plab,-1.21) +
01371                          0.169*sqrLogPlab - 1.85*LogPlab;
01372 
01373        Xtotal          = ( NumberOfTargetProtons  * XtotPP +
01374                            NumberOfTargetNeutrons * XtotPN  );
01375 
01376        Xelastic        = ( NumberOfTargetProtons  * XelPP  +
01377                            NumberOfTargetNeutrons * XelPN   );
01378   }
01379   Xinelastic = Xtotal - Xelastic;
01380 
01381   if( Xinelastic < 0.) Xinelastic = 0.;
01382 
01383   return Xinelastic*= millibarn;
01384 }

G4double G4GlauberGribovCrossSection::GetInelasticGlauberGribov ( const G4DynamicParticle ,
G4int  Z,
G4int  A 
) [inline]

Definition at line 188 of file G4GlauberGribovCrossSection.hh.

References GetIsoCrossSection().

Referenced by G4CrossSectionPairGG::BuildPhysicsTable(), G4BGGPionInelasticXS::BuildPhysicsTable(), G4BGGNucleonInelasticXS::BuildPhysicsTable(), G4CrossSectionPairGG::GetElementCrossSection(), G4BGGPionInelasticXS::GetElementCrossSection(), and G4BGGNucleonInelasticXS::GetElementCrossSection().

00190 {
00191   GetIsoCrossSection(dp, Z, A);
00192   return fInelasticXsc;
00193 }

G4double G4GlauberGribovCrossSection::GetInelasticGlauberGribovXsc (  )  [inline]

Definition at line 105 of file G4GlauberGribovCrossSection.hh.

Referenced by G4NeutronInelasticXS::GetElementCrossSection().

00105 { return fInelasticXsc; }; 

G4double G4GlauberGribovCrossSection::GetIsoCrossSection ( const G4DynamicParticle ,
G4int  Z,
G4int  A,
const G4Isotope iso = 0,
const G4Element elm = 0,
const G4Material mat = 0 
) [virtual]

Reimplemented from G4VCrossSectionDataSet.

Definition at line 319 of file G4GlauberGribovCrossSection.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), G4HadronNucleonXsc::GetHadronNucleonXscNS(), GetHNinelasticXsc(), G4HadronNucleonXsc::GetInelasticHadronNucleonXsc(), GetKaonNucleonXscVector(), GetNucleusRadius(), GetParticleBarCorIn(), GetParticleBarCorTot(), and G4INCL::Math::pi.

Referenced by GetElasticGlauberGribov(), G4NeutronInelasticXS::GetElementCrossSection(), G4NeutronElasticXS::GetElementCrossSection(), and GetInelasticGlauberGribov().

00324 {
00325   G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00326   G4double hpInXsc(0.), hnInXsc(0.);
00327   G4double R             = GetNucleusRadius(A); 
00328 
00329   G4int N = A - Z;              // number of neutrons
00330   if (N < 0) N = 0;
00331 
00332   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00333 
00334   if( theParticle == theProton  || 
00335       theParticle == theNeutron ||
00336       theParticle == thePiPlus  || 
00337       theParticle == thePiMinus      )
00338   {
00339     // sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00340 
00341     sigma = Z*hnXsc->GetHadronNucleonXscNS(aParticle, theProton);
00342 
00343     hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
00344 
00345     sigma += N*hnXsc->GetHadronNucleonXscNS(aParticle, theNeutron);
00346 
00347     hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
00348 
00349     cofInelastic = 2.4;
00350     cofTotal     = 2.0;
00351   }
00352   else if( theParticle == theKPlus   || 
00353            theParticle == theKMinus  || 
00354            theParticle == theK0S     || 
00355            theParticle == theK0L        ) 
00356   {
00357     sigma        = GetKaonNucleonXscVector(aParticle, A, Z);
00358     cofInelastic = 2.2;
00359     cofTotal     = 2.0;
00360     R = 1.3*fermi;
00361     R *= std::pow(G4double(A), 0.3333);
00362   }
00363   else
00364   {
00365     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00366     cofInelastic = 2.2;
00367     cofTotal     = 2.0;
00368   }
00369   // cofInelastic = 2.0;
00370 
00371   if( A > 1 )
00372   { 
00373     nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
00374     ratio = sigma/nucleusSquare;
00375 
00376     xsection =  nucleusSquare*std::log( 1. + ratio );
00377 
00378     xsection *= GetParticleBarCorTot(theParticle, Z);
00379 
00380     fTotalXsc = xsection;
00381 
00382   
00383 
00384     fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00385 
00386     fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);
00387 
00388     fElasticXsc   = fTotalXsc - fInelasticXsc;
00389 
00390     if(fElasticXsc < 0.) fElasticXsc = 0.;
00391     
00392     G4double difratio = ratio/(1.+ratio);
00393 
00394     fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
00395 
00396 
00397     // sigma = GetHNinelasticXsc(aParticle, A, Z);
00398 
00399     sigma = Z*hpInXsc + N*hnInXsc;
00400 
00401     ratio = sigma/nucleusSquare;
00402 
00403     fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00404 
00405     if (fElasticXsc < 0.) fElasticXsc = 0.;
00406   }
00407   else // H
00408   {
00409     fTotalXsc = sigma;
00410     xsection  = sigma;
00411     
00412     if ( theParticle != theAProton ) 
00413     {
00414       sigma         = GetHNinelasticXsc(aParticle, A, Z);
00415       fInelasticXsc = sigma;
00416       fElasticXsc   = fTotalXsc - fInelasticXsc;      
00417     }
00418     else
00419     {
00420       fElasticXsc   = fTotalXsc - fInelasticXsc;
00421     }
00422     if (fElasticXsc < 0.) fElasticXsc = 0.;
00423       
00424   }
00425   return xsection; 
00426 }

G4double G4GlauberGribovCrossSection::GetKaonNucleonXscVector ( const G4DynamicParticle ,
G4int  At,
G4int  Zt 
)

Definition at line 1099 of file G4GlauberGribovCrossSection.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscPDG(), G4DynamicParticle::GetKineticEnergy(), G4HadronNucleonXsc::GetKmNeutronTotXscVector(), G4HadronNucleonXsc::GetKmProtonTotXscVector(), G4HadronNucleonXsc::GetKpNeutronTotXscVector(), and G4HadronNucleonXsc::GetKpProtonTotXscVector().

Referenced by GetIsoCrossSection().

01101 {
01102   G4double Tkin, logTkin, xsc, xscP, xscN;
01103   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01104 
01105   G4int Nt = At-Zt;              // number of neutrons
01106   if (Nt < 0) Nt = 0;  
01107 
01108   Tkin = aParticle->GetKineticEnergy(); // Tkin in MeV
01109 
01110   if( Tkin > 70*GeV ) return GetHadronNucleonXscPDG(aParticle,At,Zt);
01111 
01112   logTkin = std::log(Tkin); // Tkin in MeV!!!
01113 
01114  if( theParticle == theKPlus )
01115  {
01116    xscP = hnXsc->GetKpProtonTotXscVector(logTkin);
01117    xscN = hnXsc->GetKpNeutronTotXscVector(logTkin);
01118  }
01119  else if( theParticle == theKMinus )
01120  {
01121    xscP = hnXsc->GetKmProtonTotXscVector(logTkin);
01122    xscN = hnXsc->GetKmNeutronTotXscVector(logTkin);
01123  }
01124  else // K-zero as half of K+ and K-
01125  {
01126    xscP = (hnXsc->GetKpProtonTotXscVector(logTkin)+hnXsc->GetKmProtonTotXscVector(logTkin))*0.5;
01127    xscN = (hnXsc->GetKpNeutronTotXscVector(logTkin)+hnXsc->GetKmNeutronTotXscVector(logTkin))*0.5;
01128  }
01129  xsc = xscP*Zt + xscN*Nt;
01130   return xsc;
01131 }

G4double G4GlauberGribovCrossSection::GetNucleusRadius ( G4int  At  ) 

Definition at line 1448 of file G4GlauberGribovCrossSection.cc.

References G4INCL::Math::oneThird.

01449 {
01450   G4double oneThird = 1.0/3.0;
01451   G4double cubicrAt = std::pow(G4double(At), oneThird); 
01452 
01453   G4double R;  // = fRadiusConst*cubicrAt;
01454 
01455   /*
01456   G4double tmp = std::pow( cubicrAt-1., 3.);
01457   tmp         += At;
01458   tmp         *= 0.5;
01459 
01460   if (At > 20.)
01461   {
01462     R = fRadiusConst*std::pow (tmp, oneThird); 
01463   }
01464   else
01465   {
01466     R = fRadiusConst*cubicrAt; 
01467   }
01468   */
01469 
01470   R = fRadiusConst*cubicrAt;
01471 
01472   G4double meanA = 20.;
01473   G4double tauA  = 20.; 
01474 
01475   if (At > 20)   // 20.
01476   {
01477     R *= ( 0.8 + 0.2*std::exp( -(G4double(At) - meanA)/tauA) ); 
01478   }
01479   else
01480   {
01481     R *= ( 1.0 + 0.1*( 1. - std::exp( (G4double(At) - meanA)/tauA) ) ); 
01482   }
01483 
01484   return R;
01485 }

G4double G4GlauberGribovCrossSection::GetNucleusRadius ( const G4DynamicParticle ,
const G4Element  
)

Definition at line 1391 of file G4GlauberGribovCrossSection.cc.

References G4lrint(), G4Element::GetN(), and G4INCL::Math::oneThird.

Referenced by GetIsoCrossSection(), GetRatioQE(), and GetRatioSD().

01393 {
01394   G4int At = G4lrint(anElement->GetN());
01395   G4double oneThird = 1.0/3.0;
01396   G4double cubicrAt = std::pow(G4double(At), oneThird); 
01397 
01398   G4double R;  // = fRadiusConst*cubicrAt;
01399   /*  
01400   G4double tmp = std::pow( cubicrAt-1., 3.);
01401   tmp         += At;
01402   tmp         *= 0.5;
01403 
01404   if (At > 20.)   // 20.
01405   {
01406     R = fRadiusConst*std::pow (tmp, oneThird); 
01407   }
01408   else
01409   {
01410     R = fRadiusConst*cubicrAt; 
01411   }
01412   */
01413   
01414   R = fRadiusConst*cubicrAt;
01415 
01416   G4double meanA  = 21.;
01417 
01418   G4double tauA1  = 40.; 
01419   G4double tauA2  = 10.; 
01420   G4double tauA3  = 5.; 
01421 
01422   G4double a1 = 0.85;
01423   G4double b1 = 1. - a1;
01424 
01425   G4double b2 = 0.3;
01426   G4double b3 = 4.;
01427 
01428   if (At > 20)   // 20.
01429   {
01430     R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) ); 
01431   }
01432   else if (At > 3)
01433   {
01434     R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) ); 
01435   }
01436   else 
01437   {
01438     R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) ); 
01439   }  
01440   return R;
01441  
01442 }

G4double G4GlauberGribovCrossSection::GetParticleBarCorIn ( const G4ParticleDefinition theParticle,
G4int  Z 
) [inline]

Definition at line 221 of file G4GlauberGribovCrossSection.hh.

Referenced by GetIsoCrossSection().

00223 {
00224   if(Z >= 2 && Z <= 92)
00225   {
00226     if(      theParticle == theProton ) return fProtonBarCorrectionIn[Z]; 
00227     else if( theParticle == theNeutron) return fNeutronBarCorrectionIn[Z]; 
00228     else if( theParticle == thePiPlus ) return fPionPlusBarCorrectionIn[Z];
00229     else if( theParticle == thePiMinus) return fPionMinusBarCorrectionIn[Z];
00230     else return 1.0;
00231   }
00232   else return 1.0;
00233 }

G4double G4GlauberGribovCrossSection::GetParticleBarCorTot ( const G4ParticleDefinition theParticle,
G4int  Z 
) [inline]

Definition at line 201 of file G4GlauberGribovCrossSection.hh.

Referenced by GetIsoCrossSection().

00203 {
00204   if(Z >= 2 && Z <= 92)
00205   {
00206     if(      theParticle == theProton ) return fProtonBarCorrectionTot[Z]; 
00207     else if( theParticle == theNeutron) return fNeutronBarCorrectionTot[Z]; 
00208     else if( theParticle == thePiPlus ) return fPionPlusBarCorrectionTot[Z];
00209     else if( theParticle == thePiMinus) return fPionMinusBarCorrectionTot[Z];
00210     else return 1.0;
00211   }
00212   else return 1.0;
00213 }

G4double G4GlauberGribovCrossSection::GetProductionGlauberGribovXsc (  )  [inline]

Definition at line 106 of file G4GlauberGribovCrossSection.hh.

00106 { return fProductionXsc; }; 

G4double G4GlauberGribovCrossSection::GetRadiusConst (  )  [inline]

Definition at line 108 of file G4GlauberGribovCrossSection.hh.

00108 { return fRadiusConst;  }; 

G4double G4GlauberGribovCrossSection::GetRatioQE ( const G4DynamicParticle ,
G4int  At,
G4int  Zt 
)

Definition at line 475 of file G4GlauberGribovCrossSection.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), GetHNinelasticXsc(), GetNucleusRadius(), and G4INCL::Math::pi.

00476 {
00477   G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00478   G4double R             = GetNucleusRadius(A); 
00479 
00480   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00481 
00482   if( theParticle == theProton  || 
00483       theParticle == theNeutron ||
00484       theParticle == thePiPlus  || 
00485       theParticle == thePiMinus      )
00486   {
00487     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00488     cofInelastic = 2.4;
00489     cofTotal     = 2.0;
00490   }
00491   else
00492   {
00493     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00494     cofInelastic = 2.2;
00495     cofTotal     = 2.0;
00496   }
00497   nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
00498   ratio = sigma/nucleusSquare;
00499 
00500   fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00501 
00502   sigma = GetHNinelasticXsc(aParticle, A, Z);
00503   ratio = sigma/nucleusSquare;
00504 
00505   fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00506 
00507   if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
00508   else                                ratio = 0.;
00509   if ( ratio < 0. )                   ratio = 0.;
00510 
00511   return ratio; 
00512 }

G4double G4GlauberGribovCrossSection::GetRatioSD ( const G4DynamicParticle ,
G4int  At,
G4int  Zt 
)

Definition at line 433 of file G4GlauberGribovCrossSection.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), GetNucleusRadius(), and G4INCL::Math::pi.

00434 {
00435   G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00436   G4double R             = GetNucleusRadius(A); 
00437 
00438   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00439 
00440   if( theParticle == theProton  || 
00441       theParticle == theNeutron ||
00442       theParticle == thePiPlus  || 
00443       theParticle == thePiMinus      )
00444   {
00445     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00446     cofInelastic = 2.4;
00447     cofTotal     = 2.0;
00448   }
00449   else
00450   {
00451     sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
00452     cofInelastic = 2.2;
00453     cofTotal     = 2.0;
00454   }
00455   nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
00456   ratio = sigma/nucleusSquare;
00457 
00458   fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00459    
00460   G4double difratio = ratio/(1.+ratio);
00461 
00462   fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
00463 
00464   if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
00465   else                    ratio = 0.;
00466 
00467   return ratio; 
00468 }

G4double G4GlauberGribovCrossSection::GetTotalGlauberGribovXsc (  )  [inline]

Definition at line 103 of file G4GlauberGribovCrossSection.hh.

00103 { return fTotalXsc;     }; 

G4bool G4GlauberGribovCrossSection::IsIsoApplicable ( const G4DynamicParticle aDP,
G4int  Z,
G4int  A,
const G4Element elm = 0,
const G4Material mat = 0 
) [virtual]

Reimplemented from G4VCrossSectionDataSet.

Definition at line 283 of file G4GlauberGribovCrossSection.cc.

References G4DynamicParticle::GetDefinition(), and G4DynamicParticle::GetKineticEnergy().

00287 {
00288   G4bool applicable      = false;
00289   // G4int baryonNumber     = aDP->GetDefinition()->GetBaryonNumber();
00290   G4double kineticEnergy = aDP->GetKineticEnergy();
00291 
00292   const G4ParticleDefinition* theParticle = aDP->GetDefinition();
00293  
00294   if ( ( kineticEnergy  >= fLowerLimit &&
00295          Z > 1 &&      // >=  He
00296        ( theParticle == theAProton   ||
00297          theParticle == theGamma     ||
00298          theParticle == theKPlus     ||
00299          theParticle == theKMinus    || 
00300          theParticle == theK0L     ||
00301          theParticle == theK0S    || 
00302          theParticle == theSMinus    ||  
00303          theParticle == theProton    ||
00304          theParticle == theNeutron   ||   
00305          theParticle == thePiPlus    ||
00306          theParticle == thePiMinus       ) )    ) applicable = true;
00307 
00308   return applicable;
00309 }

void G4GlauberGribovCrossSection::SetEnergyLowerLimit ( G4double  E  )  [inline]

Definition at line 113 of file G4GlauberGribovCrossSection.hh.

00113 {fLowerLimit=E;};


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:05 2013 for Geant4 by  doxygen 1.4.7