Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4HadronNucleonXsc Class Reference

#include <G4HadronNucleonXsc.hh>

Public Member Functions

 G4HadronNucleonXsc ()
 
virtual ~G4HadronNucleonXsc ()
 
virtual G4bool IsApplicable (const G4DynamicParticle *aDP, const G4Element *)
 
virtual G4bool IsIsoApplicable (const G4DynamicParticle *aDP, G4int Z, G4int A)
 
virtual void DumpPhysicsTable (const G4ParticleDefinition &)
 
void CrossSectionDescription (std::ostream &) const
 
G4double GetHadronNucleonXscEL (const G4DynamicParticle *, const G4ParticleDefinition *)
 
G4double GetHadronNucleonXscPDG (const G4DynamicParticle *, const G4ParticleDefinition *)
 
G4double GetHadronNucleonXscNS (const G4DynamicParticle *, const G4ParticleDefinition *)
 
G4double GetKaonNucleonXscGG (const G4DynamicParticle *, const G4ParticleDefinition *)
 
G4double GetHadronNucleonXscVU (const G4DynamicParticle *, const G4ParticleDefinition *)
 
G4double CalculateEcmValue (const G4double, const G4double, const G4double)
 
G4double CalcMandelstamS (const G4double, const G4double, const G4double)
 
G4double GetCoulombBarrier (const G4DynamicParticle *aParticle, const G4ParticleDefinition *nucleon)
 
G4double GetTotalHadronNucleonXsc ()
 
G4double GetElasticHadronNucleonXsc ()
 
G4double GetInelasticHadronNucleonXsc ()
 
void InitialiseKaonNucleonTotXsc ()
 
G4double GetKpProtonTotXscVector (G4double logEnergy)
 
G4double GetKpNeutronTotXscVector (G4double logEnergy)
 
G4double GetKmProtonTotXscVector (G4double logEnergy)
 
G4double GetKmNeutronTotXscVector (G4double logEnergy)
 

Detailed Description

Definition at line 51 of file G4HadronNucleonXsc.hh.

Constructor & Destructor Documentation

G4HadronNucleonXsc::G4HadronNucleonXsc ( )

Definition at line 39 of file G4HadronNucleonXsc.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(), InitialiseKaonNucleonTotXsc(), 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().

40 :
41 // fUpperLimit( 10000 * GeV ),
42  fLowerLimit( 0.03 * MeV ),
43  fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0)
44 // , fHadronNucleonXsc(0.0)
45 {
46  theGamma = G4Gamma::Gamma();
47  theProton = G4Proton::Proton();
48  theNeutron = G4Neutron::Neutron();
49  theAProton = G4AntiProton::AntiProton();
50  theANeutron = G4AntiNeutron::AntiNeutron();
51  thePiPlus = G4PionPlus::PionPlus();
52  thePiMinus = G4PionMinus::PionMinus();
53  thePiZero = G4PionZero::PionZero();
54  theKPlus = G4KaonPlus::KaonPlus();
55  theKMinus = G4KaonMinus::KaonMinus();
56  theK0S = G4KaonZeroShort::KaonZeroShort();
57  theK0L = G4KaonZeroLong::KaonZeroLong();
58  theL = G4Lambda::Lambda();
59  theAntiL = G4AntiLambda::AntiLambda();
60  theSPlus = G4SigmaPlus::SigmaPlus();
61  theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
62  theSMinus = G4SigmaMinus::SigmaMinus();
63  theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
64  theS0 = G4SigmaZero::SigmaZero();
65  theAS0 = G4AntiSigmaZero::AntiSigmaZero();
66  theXiMinus = G4XiMinus::XiMinus();
67  theXi0 = G4XiZero::XiZero();
68  theAXiMinus = G4AntiXiMinus::AntiXiMinus();
69  theAXi0 = G4AntiXiZero::AntiXiZero();
70  theOmega = G4OmegaMinus::OmegaMinus();
71  theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
72  theD = G4Deuteron::Deuteron();
73  theT = G4Triton::Triton();
74  theA = G4Alpha::Alpha();
75  theHe3 = G4He3::He3();
76 
77  InitialiseKaonNucleonTotXsc();
78 }
G4AntiOmegaMinus::AntiOmegaMinus
static G4AntiOmegaMinus * AntiOmegaMinus()
Definition: G4AntiOmegaMinus.cc:110
G4OmegaMinus::OmegaMinus
static G4OmegaMinus * OmegaMinus()
Definition: G4OmegaMinus.cc:110
G4KaonZeroLong::KaonZeroLong
static G4KaonZeroLong * KaonZeroLong()
Definition: G4KaonZeroLong.cc:115
G4AntiSigmaPlus::AntiSigmaPlus
static G4AntiSigmaPlus * AntiSigmaPlus()
Definition: G4AntiSigmaPlus.cc:108
G4SigmaZero::SigmaZero
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:99
G4KaonMinus::KaonMinus
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
G4AntiSigmaMinus::AntiSigmaMinus
static G4AntiSigmaMinus * AntiSigmaMinus()
Definition: G4AntiSigmaMinus.cc:106
G4XiZero::XiZero
static G4XiZero * XiZero()
Definition: G4XiZero.cc:106
G4KaonZeroShort::KaonZeroShort
static G4KaonZeroShort * KaonZeroShort()
Definition: G4KaonZeroShort.cc:104
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4XiMinus::XiMinus
static G4XiMinus * XiMinus()
Definition: G4XiMinus.cc:106
G4AntiXiMinus::AntiXiMinus
static G4AntiXiMinus * AntiXiMinus()
Definition: G4AntiXiMinus.cc:105
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4SigmaMinus::SigmaMinus
static G4SigmaMinus * SigmaMinus()
Definition: G4SigmaMinus.cc:106
G4AntiLambda::AntiLambda
static G4AntiLambda * AntiLambda()
Definition: G4AntiLambda.cc:108
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4AntiSigmaZero::AntiSigmaZero
static G4AntiSigmaZero * AntiSigmaZero()
Definition: G4AntiSigmaZero.cc:99
G4AntiXiZero::AntiXiZero
static G4AntiXiZero * AntiXiZero()
Definition: G4AntiXiZero.cc:106
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4SigmaPlus::SigmaPlus
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
G4Lambda::Lambda
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
G4KaonPlus::KaonPlus
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
G4HadronNucleonXsc::~G4HadronNucleonXsc ( )
virtual

Definition at line 81 of file G4HadronNucleonXsc.cc.

82 {}

Member Function Documentation

G4double G4HadronNucleonXsc::CalcMandelstamS ( const G4double  mp,
const G4double  mt,
const G4double  Plab 
)

Definition at line 1398 of file G4HadronNucleonXsc.cc.

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

1401 {
1402  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1403  G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1404 
1405  return sMand;
1406 }
G4double
double G4double
Definition: G4Types.hh:76
G4double G4HadronNucleonXsc::CalculateEcmValue ( const G4double  mp,
const G4double  mt,
const G4double  Plab 
)

Definition at line 1381 of file G4HadronNucleonXsc.cc.

1384 {
1385  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1386  G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
1387  // G4double Pcm = Plab * mt / Ecm;
1388  // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
1389 
1390  return Ecm ; // KEcm;
1391 }
G4double
double G4double
Definition: G4Types.hh:76
void G4HadronNucleonXsc::CrossSectionDescription ( std::ostream &  outFile) const

Definition at line 84 of file G4HadronNucleonXsc.cc.

85 {
86  outFile << "G4HadronNucleonXsc calculates the total, inelastic and elastic\n"
87  << "hadron-nucleon cross sections using several different\n"
88  << "parameterizations within the Glauber-Gribov approach. It is\n"
89  << "valid for all incident gammas and long-lived hadrons at\n"
90  << "energies above 30 keV. This is a cross section component which\n"
91  << "is to be used to build a cross section data set.\n";
92 }
outFile
std::ofstream outFile
Definition: GammaRayTel.cc:68
virtual void G4HadronNucleonXsc::DumpPhysicsTable ( const G4ParticleDefinition &  )
inlinevirtual

Definition at line 65 of file G4HadronNucleonXsc.hh.

References G4cout, and G4endl.

66  {G4cout << "G4HadronNucleonXsc: uses parametrisation"<<G4endl;}
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double G4HadronNucleonXsc::GetCoulombBarrier ( const G4DynamicParticle *  aParticle,
const G4ParticleDefinition *  nucleon 
)

Definition at line 1413 of file G4HadronNucleonXsc.cc.

References python.hepunit::fermi, python.hepunit::fine_structure_const, G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and python.hepunit::hbarc.

Referenced by GetHadronNucleonXscNS(), and GetKaonNucleonXscGG().

1415 {
1416  G4double ratio;
1417 
1418  G4double tR = 0.895*fermi, pR;
1419 
1420  if ( aParticle->GetDefinition() == theProton ) pR = 0.895*fermi;
1421  else if( aParticle->GetDefinition() == thePiPlus ) pR = 0.663*fermi;
1422  else if( aParticle->GetDefinition() == theKPlus ) pR = 0.340*fermi;
1423  else pR = 0.500*fermi;
1424 
1425  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
1426  G4double tZ = nucleon->GetPDGCharge();
1427 
1428  G4double pTkin = aParticle->GetKineticEnergy();
1429 
1430  G4double pM = aParticle->GetDefinition()->GetPDGMass();
1431  G4double tM = nucleon->GetPDGMass();
1432 
1433  G4double pElab = pTkin + pM;
1434 
1435  G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
1436 
1437  G4double totTcm = totEcm - pM -tM;
1438 
1439  G4double bC = fine_structure_const*hbarc*pZ*tZ;
1440  bC /= pR + tR;
1441  bC /= 2.; // 4., 2. parametrisation cof ??? vmg
1442 
1443  // G4cout<<"pTkin = "<<pTkin/GeV<<"; pPlab = "
1444  // <<pPlab/GeV<<"; bC = "<<bC/GeV<<"; pTcm = "<<pTcm/GeV<<G4endl;
1445 
1446  if( totTcm <= bC ) ratio = 0.;
1447  else ratio = 1. - bC/totTcm;
1448 
1449  // if(ratio < DBL_MIN) ratio = DBL_MIN;
1450  if( ratio < 0.) ratio = 0.;
1451 
1452  // G4cout <<"ratio = "<<ratio<<G4endl;
1453  return ratio;
1454 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4double G4HadronNucleonXsc::GetElasticHadronNucleonXsc ( )
inline

Definition at line 91 of file G4HadronNucleonXsc.hh.

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

91 { return fElasticXsc; };
G4double G4HadronNucleonXsc::GetHadronNucleonXscEL ( const G4DynamicParticle *  aParticle,
const G4ParticleDefinition *  nucleon 
)

Definition at line 143 of file G4HadronNucleonXsc.cc.

References CalcMandelstamS(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), python.hepunit::GeV, CLHEP::Hep3Vector::mag(), and python.hepunit::millibarn.

145 {
146  G4double xsection;
147 
148 
149  G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
150 
151  G4double proj_mass = aParticle->GetMass();
152  G4double proj_momentum = aParticle->GetMomentum().mag();
153  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
154 
155  sMand /= GeV*GeV; // in GeV for parametrisation
156  proj_momentum /= GeV;
157 
158  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
159 
160  G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
161 
162 
163  if(theParticle == theGamma && pORn )
164  {
165  xsection = (0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
166  }
167  else if(theParticle == theNeutron && pORn ) // as proton ???
168  {
169  xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
170  }
171  else if(theParticle == theProton && pORn )
172  {
173  xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
174 
175  // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
176  // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
177  }
178  else if(theParticle == theAProton && pORn )
179  {
180  xsection = ( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
181  }
182  else if(theParticle == thePiPlus && pORn )
183  {
184  xsection = (13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
185  }
186  else if(theParticle == thePiMinus && pORn )
187  {
188  // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
189  xsection = (13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
190  }
191  else if(theParticle == theKPlus && pORn )
192  {
193  xsection = (11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
194  }
195  else if(theParticle == theKMinus && pORn )
196  {
197  xsection = (11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
198  }
199  else // as proton ???
200  {
201  xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
202  }
203  xsection *= millibarn;
204 
205  fTotalXsc = xsection;
206  fInelasticXsc = 0.83*xsection;
207  fElasticXsc = fTotalXsc - fInelasticXsc;
208  if (fElasticXsc < 0.)fElasticXsc = 0.;
209 
210  return xsection;
211 }
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
python.hepunit.millibarn
int millibarn
Definition: hepunit.py:40
G4DynamicParticle::GetMass
G4double GetMass() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4HadronNucleonXsc::CalcMandelstamS
G4double CalcMandelstamS(const G4double, const G4double, const G4double)
Definition: G4HadronNucleonXsc.cc:1398
G4double G4HadronNucleonXsc::GetHadronNucleonXscNS ( const G4DynamicParticle *  aParticle,
const G4ParticleDefinition *  nucleon 
)

pi- ////////////////////////////////////////////

Definition at line 365 of file G4HadronNucleonXsc.cc.

References CalcMandelstamS(), GetCoulombBarrier(), G4DynamicParticle::GetDefinition(), GetHadronNucleonXscPDG(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), G4ParticleDefinition::GetPDGCharge(), G4DynamicParticle::GetTotalEnergy(), python.hepunit::GeV, LE, CLHEP::Hep3Vector::mag(), python.hepunit::millibarn, neutron, G4InuclParticleNames::proton, G4InuclParticleNames::s0, and G4InuclParticleNames::sp.

Referenced by G4BGGNucleonInelasticXS::BuildPhysicsTable(), G4GlauberGribovCrossSection::GetIsoCrossSection(), G4BGGNucleonInelasticXS::GetIsoCrossSection(), G4ComponentGGHadronNucleusXsc::GetIsoCrossSection(), G4GGNuclNuclCrossSection::GetZandACrossSection(), and G4ComponentGGNuclNuclXsc::GetZandACrossSection().

367 {
368  G4double xsection(0);
369 
370  G4double A0, B0;
371  G4double hpXsc(0);
372  G4double hnXsc(0);
373 
374 
375  G4double tM = 0.939*GeV; // ~mean neutron and proton ???
376 
377  G4double pM = aParticle->GetMass();
378  G4double pE = aParticle->GetTotalEnergy();
379  G4double pLab = aParticle->GetMomentum().mag();
380 
381  G4double sMand = CalcMandelstamS ( pM , tM , pLab );
382 
383  sMand /= GeV*GeV; // in GeV for parametrisation
384  pLab /= GeV;
385  pE /= GeV;
386  pM /= GeV;
387 
388  G4double logP = std::log(pLab);
389 
390 
391  // General PDG fit constants
392 
393  G4double s0 = 5.38*5.38; // in Gev^2
394  G4double eta1 = 0.458;
395  G4double eta2 = 0.458;
396  G4double B = 0.308;
397  G4double minLogP = 3.5; // min of (lnP-minLogP)^2
398  G4double cofLogE = .0557; // elastic (lnP-minLogP)^2
399  G4double cofLogT = .3; // total (lnP-minLogP)^2
400  G4double pMin = .1; // fast LE calculation
401  G4double pMax = 1000.; // fast HE calculation
402 
403 
404  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
405 
406  G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
407  G4bool proton = (nucleon == theProton);
408  G4bool neutron = (nucleon == theNeutron);
409 
410  if( theParticle == theNeutron && pORn )
411  {
412  if( pLab >= 373.)
413  {
414  xsection = GetHadronNucleonXscPDG(aParticle, nucleon)/millibarn;
415 
416  fElasticXsc = 6.5 + 0.308*std::pow(std::log(sMand/400.),1.65) + 9.19*std::pow(sMand,-0.458);
417 
418  fTotalXsc = xsection;
419 
420  }
421  else if( pLab >= 100.)
422  {
423  B0 = 7.5;
424  A0 = 100. - B0*std::log(3.0e7);
425 
426  xsection = A0 + B0*std::log(pE) - 11
427  // + 103*std::pow(2*0.93827*pE + pM*pM+0.93827*0.93827,-0.165); // mb
428  + 103*std::pow(sMand,-0.165); // mb
429 
430  fElasticXsc = 5.53 + 0.308*std::pow(std::log(sMand/28.9),1.1) + 9.19*std::pow(sMand,-0.458);
431 
432  fTotalXsc = xsection;
433  }
434  else if( pLab >= 10.)
435  {
436  B0 = 7.5;
437  A0 = 100. - B0*std::log(3.0e7);
438 
439  xsection = A0 + B0*std::log(pE) - 11
440  + 103*std::pow(2*0.93827*pE + pM*pM+
441  0.93827*0.93827,-0.165); // mb
442  fTotalXsc = xsection;
443  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
444  }
445  else // pLab < 10 GeV/c
446  {
447  if( neutron ) // nn to be pp
448  {
449  if( pLab < 0.4 )
450  {
451  hnXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
452  fElasticXsc = hnXsc;
453  }
454  else if( pLab < 0.73 )
455  {
456  hnXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
457  fElasticXsc = hnXsc;
458  }
459  else if( pLab < 1.05 )
460  {
461  hnXsc = 23 + 40*(std::log(pLab/0.73))*
462  (std::log(pLab/0.73));
463  fElasticXsc = 23 + 20*(std::log(pLab/0.73))*
464  (std::log(pLab/0.73));
465  }
466  else // 1.05 - 10 GeV/c
467  {
468  hnXsc = 39.0+75*(pLab - 1.2)/(std::pow(pLab,3.0) + 0.15);
469 
470  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
471  }
472  fTotalXsc = hnXsc;
473  }
474  if( proton ) // pn to be np
475  {
476  if( pLab < 0.02 )
477  {
478  hpXsc = 4100+30*std::pow(std::log(1.3/pLab),3.6); // was as pLab < 0.8
479  fElasticXsc = hpXsc;
480  }
481  else if( pLab < 0.8 )
482  {
483  hpXsc = 33+30*std::pow(std::log(pLab/1.3),4.0);
484  fElasticXsc = hpXsc;
485  }
486  else if( pLab < 1.05 )
487  {
488  hpXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
489  fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
490  }
491  else if( pLab < 1.4 )
492  {
493  hpXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
494  fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
495  }
496  else // 1.4 < pLab < 10. )
497  {
498  hpXsc = 33.3 + 20.8*(std::pow(pLab,2.0) - 1.35)/(std::pow(pLab,2.50) + 0.95);
499 
500  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
501  }
502  fTotalXsc = hpXsc;
503  }
504  }
505  }
506  else if( theParticle == theProton && pORn ) ////// proton //////////////////////////////////////////////
507  {
508  if( pLab >= 373.) // pdg due to TOTEM data
509  {
510  xsection = GetHadronNucleonXscPDG(aParticle, nucleon)/millibarn;
511 
512  fElasticXsc = 6.5 + 0.308*std::pow(std::log(sMand/400.),1.65) + 9.19*std::pow(sMand,-0.458);
513 
514  fTotalXsc = xsection;
515  }
516  else if( pLab >= 100.)
517  {
518  B0 = 7.5;
519  A0 = 100. - B0*std::log(3.0e7);
520 
521  xsection = A0 + B0*std::log(pE) - 11 + 103*std::pow(sMand,-0.165); // mb
522 
523  fElasticXsc = 5.53 + 0.308*std::pow(std::log(sMand/28.9),1.1) + 9.19*std::pow(sMand,-0.458);
524 
525  fTotalXsc = xsection;
526  }
527  else if( pLab >= 10.)
528  {
529  B0 = 7.5;
530  A0 = 100. - B0*std::log(3.0e7);
531 
532  xsection = A0 + B0*std::log(pE) - 11 + 103*std::pow(sMand,-0.165); // mb
533 
534  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
535 
536  fTotalXsc = xsection;
537  }
538  else
539  {
540  // pp
541 
542  if( proton )
543  {
544  if( pLab < 0.4 )
545  {
546  hpXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
547  fElasticXsc = hpXsc;
548  }
549  else if( pLab < 0.73 )
550  {
551  hpXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
552  fElasticXsc = hpXsc;
553  }
554  else if( pLab < 1.05 )
555  {
556  hpXsc = 23 + 40*(std::log(pLab/0.73))*
557  (std::log(pLab/0.73));
558  fElasticXsc = 23 + 20*(std::log(pLab/0.73))*
559  (std::log(pLab/0.73));
560  }
561  else // 1.05 - 10 GeV/c
562  {
563  hpXsc = 39.0+75*(pLab - 1.2)/(std::pow(pLab,3.0) + 0.15);
564 
565  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
566  }
567  fTotalXsc = hpXsc;
568  }
569  if( neutron ) // pn to be np
570  {
571  if( pLab < 0.02 )
572  {
573  hnXsc = 4100+30*std::pow(std::log(1.3/pLab),3.6); // was as pLab < 0.8
574  fElasticXsc = hnXsc;
575  }
576  else if( pLab < 0.8 )
577  {
578  hnXsc = 33+30*std::pow(std::log(pLab/1.3),4.0);
579  fElasticXsc = hnXsc;
580  }
581  else if( pLab < 1.05 )
582  {
583  hnXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
584  fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
585  }
586  else if( pLab < 1.4 )
587  {
588  hnXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
589  fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
590  }
591  else // 1.4 < pLab < 10. )
592  {
593  hnXsc = 33.3 + 20.8*(std::pow(pLab,2.0) - 1.35)/(std::pow(pLab,2.50) + 0.95);
594 
595  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
596  }
597  fTotalXsc = hnXsc;
598  }
599  }
600  }
601  else if( theParticle == theAProton && pORn ) /////////////////// p_bar ///////////////////////////
602  {
603  if( proton )
604  {
605  xsection = 35.45 + B*std::pow(std::log(sMand/s0),2.)
606  + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2);
607  }
608  if( neutron ) // ???
609  {
610  xsection = 35.80 + B*std::pow(std::log(sMand/s0),2.)
611  + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2);
612  }
613  fTotalXsc = xsection;
614  }
615  else if( theParticle == thePiPlus && pORn ) // pi+ /////////////////////////////////////////////
616  {
617  if( proton ) // pi+ p
618  {
619  if( pLab < 0.28 )
620  {
621  hpXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
622  fElasticXsc = hpXsc;
623  }
624  else if( pLab < 0.4 )
625  {
626  hpXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
627  fElasticXsc = hpXsc;
628  }
629  else if( pLab < 0.68 )
630  {
631  hpXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
632  fElasticXsc = hpXsc;
633  }
634  else if( pLab < 0.85 )
635  {
636  G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
637  hpXsc = Ex4 + 14.9;
638  fElasticXsc = hpXsc*std::exp(-3.*(pLab - 0.68));
639  }
640  else if( pLab < 1.15 )
641  {
642  G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
643  hpXsc = Ex4 + 14.9;
644 
645  fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
646  }
647  else if( pLab < 1.4) // ns original
648  {
649  G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
650  G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
651  hpXsc = Ex1 + Ex2 + 27.5;
652  fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
653  }
654  else if( pLab < 2.0 ) // ns original
655  {
656  G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
657  G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
658  hpXsc = Ex1 + Ex2 + 27.5;
659  fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
660  }
661  else if( pLab < 3.5 ) // ns original
662  {
663  G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
664  G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
665  hpXsc = Ex1 + Ex2 + 27.5;
666  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
667  }
668  else if( pLab < 200. ) // my
669  {
670  hpXsc = 10.6 + 2.*std::log(pE) + 25*std::pow(pE, -0.43 ); // ns original
671  // hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
672  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
673  }
674  else // pLab > 100 // my
675  {
676  hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
677  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
678  }
679  fTotalXsc = hpXsc;
680  }
681  if( neutron ) // pi+ n = pi- p??
682  {
683  if( pLab < 0.28 )
684  {
685  hnXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
686  fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
687  }
688  else if( pLab < 0.395676 ) // first peak
689  {
690  hnXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
691  fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
692  }
693  else if( pLab < 0.5 )
694  {
695  hnXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
696  fElasticXsc = 0.37*hnXsc;
697  }
698  else if( pLab < 0.65 )
699  {
700  hnXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
701  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
702  }
703  else if( pLab < 0.72 )
704  {
705  hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
706  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
707  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
708  }
709  else if( pLab < 0.88 )
710  {
711  hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
712  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
713  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
714  }
715  else if( pLab < 1.03 )
716  {
717  hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
718  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
719  fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
720  }
721  else if( pLab < 1.15 )
722  {
723  hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
724  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
725  fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
726  }
727  else if( pLab < 1.3 )
728  {
729  hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
730  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
731  fElasticXsc = 3. + 13./pLab;
732  }
733  else if( pLab < 2.6 ) // < 3.0) // ns original
734  {
735  hnXsc = 36.1 + 0.079-4.313*std::log(pLab)+
736  3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
737  1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
738  fElasticXsc = 3. + 13./pLab;
739  }
740  else if( pLab < 20. ) // < 3.0) // ns original
741  {
742  hnXsc = 36.1 + 0.079 - 4.313*std::log(pLab)+
743  3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
744  1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
745  fElasticXsc = 3. + 13./pLab;
746  }
747  else // mb
748  {
749  hnXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
750  fElasticXsc = 3. + 13./pLab;
751  }
752  fTotalXsc = hnXsc;
753  }
754  }
755  else if( theParticle == thePiMinus && pORn ) /// pi- ////////////////////////////////////////////
756  {
757  if( neutron ) // pi- n = pi+ p??
758  {
759  if( pLab < 0.28 )
760  {
761  hnXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
762  fElasticXsc = hnXsc;
763  }
764  else if( pLab < 0.4 )
765  {
766  hnXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
767  fElasticXsc = hnXsc;
768  }
769  else if( pLab < 0.68 )
770  {
771  hnXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
772  fElasticXsc = hnXsc;
773  }
774  else if( pLab < 0.85 )
775  {
776  G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
777  hnXsc = Ex4 + 14.9;
778  fElasticXsc = hnXsc*std::exp(-3.*(pLab - 0.68));
779  }
780  else if( pLab < 1.15 )
781  {
782  G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
783  hnXsc = Ex4 + 14.9;
784 
785  fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
786  }
787  else if( pLab < 1.4) // ns original
788  {
789  G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
790  G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
791  hnXsc = Ex1 + Ex2 + 27.5;
792  fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
793  }
794  else if( pLab < 2.0 ) // ns original
795  {
796  G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
797  G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
798  hnXsc = Ex1 + Ex2 + 27.5;
799  fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
800  }
801  else if( pLab < 3.5 ) // ns original
802  {
803  G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
804  G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
805  hnXsc = Ex1 + Ex2 + 27.5;
806  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
807  }
808  else if( pLab < 200. ) // my
809  {
810  hnXsc = 10.6 + 2.*std::log(pE) + 25*std::pow(pE, -0.43 ); // ns original
811  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
812  }
813  else // pLab > 100 // my
814  {
815  hnXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
816  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
817  }
818  fTotalXsc = hnXsc;
819  }
820  if( proton ) // pi- p
821  {
822  if( pLab < 0.28 )
823  {
824  hpXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
825  fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
826  }
827  else if( pLab < 0.395676 ) // first peak
828  {
829  hpXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
830  fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
831  }
832  else if( pLab < 0.5 )
833  {
834  hpXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
835  fElasticXsc = 0.37*hpXsc;
836  }
837  else if( pLab < 0.65 )
838  {
839  hpXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
840  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
841  }
842  else if( pLab < 0.72 )
843  {
844  hpXsc = 36.1+
845  10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
846  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
847  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
848  }
849  else if( pLab < 0.88 )
850  {
851  hpXsc = 36.1+
852  10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
853  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
854  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
855  }
856  else if( pLab < 1.03 )
857  {
858  hpXsc = 36.1+
859  10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
860  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
861  fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
862  }
863  else if( pLab < 1.15 )
864  {
865  hpXsc = 36.1+
866  10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
867  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
868  fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
869  }
870  else if( pLab < 1.3 )
871  {
872  hpXsc = 36.1+
873  10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
874  24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
875  fElasticXsc = 3. + 13./pLab;
876  }
877  else if( pLab < 2.6 ) // < 3.0) // ns original
878  {
879  hpXsc = 36.1+0.079-4.313*std::log(pLab)+
880  3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
881  1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
882  fElasticXsc = 3. +13./pLab; // *std::log(pLab*6.79);
883  }
884  else // mb
885  {
886  hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
887  fElasticXsc = 3. + 13./pLab;
888  }
889  fTotalXsc = hpXsc;
890  }
891  }
892  else if( (theParticle == theKMinus || theParticle == theK0S) && proton ) // Kmp/K0p /////////////////////////////////
893  {
894  if( pLab < pMin)
895  {
896  G4double psp = pLab*std::sqrt(pLab);
897  fElasticXsc = 5.2/psp;
898  fTotalXsc = 14./psp;
899  }
900  else if( pLab > pMax )
901  {
902  G4double ld = std::log(pLab) - minLogP;
903  G4double ld2 = ld*ld;
904  fElasticXsc = cofLogE*ld2 + 2.23;
905  fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
906  }
907  else
908  {
909  G4double ld = std::log(pLab) - minLogP;
910  G4double ld2 = ld*ld;
911  G4double sp = std::sqrt(pLab);
912  G4double psp = pLab*sp;
913  G4double p2 = pLab*pLab;
914  G4double p4 = p2*p2;
915  G4double lm = pLab - .39;
916  G4double md = lm*lm + .000356;
917 
918  G4double lh1 = pLab - 0.78;
919  G4double hd1 = lh1*lh1 + .00166;
920 
921  G4double lh = pLab - 1.01;
922  G4double hd = lh*lh + .011;
923 
924  G4double lh2 = pLab - 1.63;
925  G4double hd2 = lh2*lh2 + .007;
926 
927  fElasticXsc = 5.2/psp + (1.1*cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4)
928  + .004/md + 0.005/hd1+ 0.01/hd2 +.15/hd; // small peaks were added
929 
930  fTotalXsc = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
931  + .006/md + 0.01/hd1+ 0.02/hd2 + .20/hd ;
932  }
933  }
934  else if( (theParticle == theKMinus || theParticle == theK0S) && neutron ) // Kmn/K0n //////////////////////////////
935  {
936  if( pLab > pMax )
937  {
938  G4double ld = std::log(pLab) - minLogP;
939  G4double ld2 = ld*ld;
940  fElasticXsc = cofLogE*ld2 + 2.23;
941  fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
942  }
943  else
944  {
945 
946  G4double lh = pLab - 0.98;
947  G4double hd = lh*lh + .021;
948 
949  G4double LogPlab = std::log( pLab );
950  G4double sqrLogPlab = LogPlab * LogPlab;
951 
952  fElasticXsc = // 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md
953  5.0 + 8.1*std::pow(pLab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab + .15/hd;
954  fTotalXsc = // 14./psp +
955  // (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
956  25.2 + 0. *std::pow(pLab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab
957  // + .006/md + 0.01/hd1+ 0.02/hd2
958  + 0.30/hd ;
959  }
960  }
961  else if( (theParticle == theKPlus || theParticle == theK0L) && proton ) // Kpp/aKp ////////////////////////
962  {
963  if( pLab < pMin )
964  {
965  G4double lr = pLab - .38;
966  G4double lm = pLab - 1.;
967  G4double md = lm*lm + .392;
968  fElasticXsc = .7/(lr*lr + .076) + 2./md;
969  fTotalXsc = .7/(lr*lr + .076) + 2.6/md;
970  }
971  else if( pLab > pMax )
972  {
973  G4double ld = std::log(pLab) - minLogP;
974  G4double ld2 = ld*ld;
975  fElasticXsc = cofLogE*ld2 + 2.23;
976  fTotalXsc = cofLogT*ld2 + 19.2;
977  }
978  else
979  {
980  G4double ld = std::log(pLab) - minLogP;
981  G4double ld2 = ld*ld;
982  G4double lr = pLab - .38;
983  G4double LE = .7/(lr*lr + .076);
984  G4double sp = std::sqrt(pLab);
985  G4double p2 = pLab*pLab;
986  G4double p4 = p2*p2;
987  G4double lm = pLab - 1.;
988  G4double md = lm*lm + .392;
989  fElasticXsc = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
990  fTotalXsc = LE + (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 2.6/md;
991  }
992  }
993  else if( (theParticle == theKPlus || theParticle == theK0L) && neutron ) // Kpn/aKn ///////////////////////
994  {
995  if( pLab < pMin )
996  {
997  G4double lm = pLab - 0.94;
998  G4double md = lm*lm + .392;
999  fElasticXsc = 2./md;
1000  fTotalXsc = 4.6/md;
1001  }
1002  else if( pLab > pMax )
1003  {
1004  G4double ld = std::log(pLab) - minLogP;
1005  G4double ld2 = ld*ld;
1006  fElasticXsc = cofLogE*ld2 + 2.23;
1007  fTotalXsc = cofLogT*ld2 + 19.2;
1008  }
1009  else
1010  {
1011  G4double ld = std::log(pLab) - minLogP;
1012  G4double ld2 = ld*ld;
1013  G4double sp = std::sqrt(pLab);
1014  G4double p2 = pLab*pLab;
1015  G4double p4 = p2*p2;
1016  G4double lm = pLab - 0.94;
1017  G4double md = lm*lm + .392;
1018  fElasticXsc = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
1019  fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 4.6/md;
1020  }
1021  }
1022  else if( theParticle == theSMinus && pORn )
1023  {
1024  xsection = 35.20 + B*std::pow(std::log(sMand/s0),2.)
1025  - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2);
1026  }
1027  else if( theParticle == theGamma && pORn ) // modify later on
1028  {
1029  xsection = 0.0 + B*std::pow(std::log(sMand/s0),2.)
1030  + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2);
1031  fTotalXsc = xsection;
1032  }
1033  else // other then p,n,pi+,pi-,K+,K- as proton ???
1034  {
1035  if( proton )
1036  {
1037  xsection = 35.45 + B*std::pow(std::log(sMand/s0),2.)
1038  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2);
1039  }
1040  if( neutron )
1041  {
1042  xsection += 35.80 + B*std::pow(std::log(sMand/s0),2.)
1043  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2);
1044  }
1045  fTotalXsc = xsection;
1046  }
1047  fTotalXsc *= millibarn; // parametrised in mb
1048  fElasticXsc *= millibarn; // parametrised in mb
1049 
1050  if( proton && aParticle->GetDefinition()->GetPDGCharge() > 0. )
1051  {
1052  G4double cB = GetCoulombBarrier(aParticle, nucleon);
1053  fTotalXsc *= cB;
1054  fElasticXsc *= cB;
1055  }
1056  fInelasticXsc = fTotalXsc - fElasticXsc;
1057  if( fInelasticXsc < 0. ) fInelasticXsc = 0.;
1058 
1059  // G4cout<<fTotalXsc/millibarn<<"; "<<fElasticXsc/millibarn<<"; "<<fInelasticXsc/millibarn<<G4endl;
1060 
1061  return fTotalXsc;
1062 }
G4HadronNucleonXsc::GetCoulombBarrier
G4double GetCoulombBarrier(const G4DynamicParticle *aParticle, const G4ParticleDefinition *nucleon)
Definition: G4HadronNucleonXsc.cc:1413
G4DynamicParticle::GetTotalEnergy
G4double GetTotalEnergy() const
neutron
Definition: G4DataQuestionaire.hh:35
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
python.hepunit.millibarn
int millibarn
Definition: hepunit.py:40
G4DynamicParticle::GetMass
G4double GetMass() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4HadronNucleonXsc::GetHadronNucleonXscPDG
G4double GetHadronNucleonXscPDG(const G4DynamicParticle *, const G4ParticleDefinition *)
Definition: G4HadronNucleonXsc.cc:221
LE
Definition: Evaluator.cc:66
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
CLHEP::Hep3Vector::mag
double mag() const
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4HadronNucleonXsc::CalcMandelstamS
G4double CalcMandelstamS(const G4double, const G4double, const G4double)
Definition: G4HadronNucleonXsc.cc:1398
G4double G4HadronNucleonXsc::GetHadronNucleonXscPDG ( const G4DynamicParticle *  aParticle,
const G4ParticleDefinition *  nucleon 
)

Definition at line 221 of file G4HadronNucleonXsc.cc.

References CalcMandelstamS(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), python.hepunit::GeV, CLHEP::Hep3Vector::mag(), python.hepunit::millibarn, neutron, G4InuclParticleNames::proton, and G4InuclParticleNames::s0.

Referenced by G4BGGNucleonInelasticXS::BuildPhysicsTable(), G4BGGPionElasticXS::BuildPhysicsTable(), G4BGGNucleonElasticXS::BuildPhysicsTable(), G4BGGPionInelasticXS::BuildPhysicsTable(), G4NeutronElasticXS::GetElementCrossSection(), G4NeutronInelasticXS::GetElementCrossSection(), GetHadronNucleonXscNS(), G4BGGPionElasticXS::GetIsoCrossSection(), G4BGGNucleonInelasticXS::GetIsoCrossSection(), G4BGGNucleonElasticXS::GetIsoCrossSection(), and G4BGGPionInelasticXS::GetIsoCrossSection().

223 {
224  G4double xsection(0);
225  G4int Zt=1, Nt=1, At=1;
226 
227  G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
228 
229  G4double proj_mass = aParticle->GetMass();
230  G4double proj_momentum = aParticle->GetMomentum().mag();
231 
232  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
233 
234  sMand /= GeV*GeV; // in GeV for parametrisation
235 
236  // General PDG fit constants
237 
238  G4double s0 = 5.38*5.38; // in Gev^2
239  G4double eta1 = 0.458;
240  G4double eta2 = 0.458;
241  G4double B = 0.308;
242 
243  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
244 
245  G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
246  G4bool proton = (nucleon == theProton);
247  G4bool neutron = (nucleon == theNeutron);
248 
249  if(theParticle == theNeutron) // proton-neutron fit
250  {
251  if ( proton )
252  {
253  xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
254  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));// on p
255  }
256  if ( neutron )
257  {
258  xsection = Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
259  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // on n pp for nn
260  }
261  }
262  else if(theParticle == theProton)
263  {
264  if ( proton )
265  {
266  xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
267  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
268  }
269  if ( neutron )
270  {
271  xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
272  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
273  }
274  }
275  else if(theParticle == theAProton)
276  {
277  if ( proton )
278  {
279  xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
280  + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
281  }
282  if ( neutron )
283  {
284  xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
285  + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
286  }
287  }
288  else if(theParticle == thePiPlus && pORn )
289  {
290  xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
291  + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
292  }
293  else if(theParticle == thePiMinus && pORn )
294  {
295  xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
296  + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
297  }
298  else if(theParticle == theKPlus)
299  {
300  if ( proton )
301  {
302  xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
303  + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
304  }
305  if ( neutron )
306  {
307  xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
308  + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
309  }
310  }
311  else if(theParticle == theKMinus)
312  {
313  if ( proton )
314  {
315  xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
316  + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
317  }
318  if ( neutron )
319  {
320  xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
321  + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2) );
322  }
323  }
324  else if(theParticle == theSMinus && pORn )
325  {
326  xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
327  - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2) );
328  }
329  else if(theParticle == theGamma && pORn ) // modify later on
330  {
331  xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
332  + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2) );
333 
334  }
335  else // as proton ???
336  {
337  if ( proton )
338  {
339  xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
340  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2) );
341  }
342  if ( neutron )
343  {
344  xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
345  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
346  }
347  }
348  xsection *= millibarn; // parametrised in mb
349 
350  fTotalXsc = xsection;
351  fInelasticXsc = 0.75*xsection;
352  fElasticXsc = fTotalXsc - fInelasticXsc;
353  if (fElasticXsc < 0.) fElasticXsc = 0.;
354 
355  return xsection;
356 }
neutron
Definition: G4DataQuestionaire.hh:35
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.millibarn
int millibarn
Definition: hepunit.py:40
G4DynamicParticle::GetMass
G4double GetMass() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4HadronNucleonXsc::CalcMandelstamS
G4double CalcMandelstamS(const G4double, const G4double, const G4double)
Definition: G4HadronNucleonXsc.cc:1398
G4double G4HadronNucleonXsc::GetHadronNucleonXscVU ( const G4DynamicParticle *  aParticle,
const G4ParticleDefinition *  nucleon 
)

Definition at line 1234 of file G4HadronNucleonXsc.cc.

References G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMomentum(), G4ParticleDefinition::GetPDGEncoding(), G4DynamicParticle::GetTotalEnergy(), python.hepunit::GeV, CLHEP::Hep3Vector::mag(), python.hepunit::millibarn, neutron, and G4InuclParticleNames::proton.

1236 {
1237  G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
1238  G4int absPDGcode = std::abs(PDGcode);
1239  G4double Elab = aParticle->GetTotalEnergy();
1240  // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
1241  G4double Plab = aParticle->GetMomentum().mag();
1242  // std::sqrt(Elab * Elab - 0.88);
1243 
1244  Elab /= GeV;
1245  Plab /= GeV;
1246 
1247  G4double LogPlab = std::log( Plab );
1248  G4double sqrLogPlab = LogPlab * LogPlab;
1249 
1250  G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
1251  G4bool proton = (nucleon == theProton);
1252  G4bool neutron = (nucleon == theNeutron);
1253 
1254 
1255  if( absPDGcode > 1000 && pORn ) //------Projectile is baryon -
1256  {
1257  if(proton)
1258  {
1259  fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
1260  fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
1261  }
1262  if(neutron)
1263  {
1264  fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
1265  fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
1266  }
1267  }
1268  else if( PDGcode == 211 && pORn ) //------Projectile is PionPlus ----
1269  {
1270  if(proton)
1271  {
1272  fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
1273  fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
1274  }
1275  if(neutron)
1276  {
1277  fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
1278  fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
1279  }
1280  }
1281  else if( PDGcode == -211 && pORn ) //------Projectile is PionMinus ----
1282  {
1283  if(proton)
1284  {
1285  fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
1286  fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
1287  }
1288  if(neutron)
1289  {
1290  fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
1291  fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
1292  }
1293  }
1294  else if( PDGcode == 111 && pORn ) //------Projectile is PionZero --
1295  {
1296  if(proton)
1297  {
1298  fTotalXsc = (16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
1299  33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
1300 
1301  fElasticXsc = ( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
1302  1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1303 
1304  }
1305  if(neutron)
1306  {
1307  fTotalXsc = (33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
1308  16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1309  fElasticXsc = ( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
1310  0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1311  }
1312  }
1313  else if( PDGcode == 321 && pORn ) //------Projectile is KaonPlus --
1314  {
1315  if(proton)
1316  {
1317  fTotalXsc = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
1318  fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
1319  }
1320  if(neutron)
1321  {
1322  fTotalXsc = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
1323  fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
1324  }
1325  }
1326  else if( PDGcode ==-321 && pORn ) //------Projectile is KaonMinus ----
1327  {
1328  if(proton)
1329  {
1330  fTotalXsc = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66*sqrLogPlab - 5.6*LogPlab;
1331  fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29*sqrLogPlab - 2.4*LogPlab;
1332  }
1333  if(neutron)
1334  {
1335  fTotalXsc = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab;
1336  fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab;
1337  }
1338  }
1339  else if( PDGcode == 311 && pORn ) //------Projectile is KaonZero -----
1340  {
1341  if(proton)
1342  {
1343  fTotalXsc = ( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
1344  32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
1345  fElasticXsc = ( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
1346  7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
1347  }
1348  if(neutron)
1349  {
1350  fTotalXsc = ( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
1351  25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
1352  fElasticXsc = ( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
1353  5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
1354  }
1355  }
1356  else //------Projectile is undefined, Nucleon assumed
1357  {
1358  if(proton)
1359  {
1360  fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
1361  fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
1362  }
1363  if(neutron)
1364  {
1365  fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
1366  fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
1367  }
1368  }
1369  fTotalXsc *= millibarn;
1370  fElasticXsc *= millibarn;
1371  fInelasticXsc = fTotalXsc - fElasticXsc;
1372  if (fInelasticXsc < 0.) fInelasticXsc = 0.;
1373 
1374  return fTotalXsc;
1375 }
G4DynamicParticle::GetTotalEnergy
G4double GetTotalEnergy() const
neutron
Definition: G4DataQuestionaire.hh:35
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.millibarn
int millibarn
Definition: hepunit.py:40
G4bool
bool G4bool
Definition: G4Types.hh:79
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4double G4HadronNucleonXsc::GetInelasticHadronNucleonXsc ( )
inline

Definition at line 92 of file G4HadronNucleonXsc.hh.

Referenced by G4BGGNucleonInelasticXS::BuildPhysicsTable(), G4BGGPionInelasticXS::BuildPhysicsTable(), G4NeutronInelasticXS::GetElementCrossSection(), G4GlauberGribovCrossSection::GetIsoCrossSection(), G4BGGNucleonInelasticXS::GetIsoCrossSection(), G4BGGPionInelasticXS::GetIsoCrossSection(), G4ComponentGGHadronNucleusXsc::GetIsoCrossSection(), G4GGNuclNuclCrossSection::GetZandACrossSection(), and G4ComponentGGNuclNuclXsc::GetZandACrossSection().

92 { return fInelasticXsc; };
G4double G4HadronNucleonXsc::GetKaonNucleonXscGG ( const G4DynamicParticle *  aParticle,
const G4ParticleDefinition *  nucleon 
)

Definition at line 1069 of file G4HadronNucleonXsc.cc.

References GetCoulombBarrier(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMomentum(), G4ParticleDefinition::GetPDGCharge(), python.hepunit::GeV, LE, CLHEP::Hep3Vector::mag(), python.hepunit::millibarn, neutron, G4InuclParticleNames::proton, and G4InuclParticleNames::sp.

Referenced by G4ComponentGGHadronNucleusXsc::GetIsoCrossSection().

1071 {
1072  G4double pLab = aParticle->GetMomentum().mag();
1073 
1074  pLab /= GeV;
1075  G4double LogPlab = std::log( pLab );
1076  G4double sqrLogPlab = LogPlab * LogPlab;
1077 
1078  G4double minLogP = 3.5; // min of (lnP-minLogP)^2
1079  G4double cofLogE = .0557; // elastic (lnP-minLogP)^2
1080  G4double cofLogT = .3; // total (lnP-minLogP)^2
1081  G4double pMin = .1; // fast LE calculation
1082  G4double pMax = 1000.; // fast HE calculation
1083 
1084  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1085 
1086  G4bool proton = (nucleon == theProton);
1087  G4bool neutron = (nucleon == theNeutron);
1088 
1089  if( (theParticle == theKMinus || theParticle == theK0S) && proton ) // (K-,K0)on p ////////////////////////////
1090  {
1091 
1092  if( pLab < pMin)
1093  {
1094  G4double psp = pLab*std::sqrt(pLab);
1095  fElasticXsc = 5.2/psp;
1096  fTotalXsc = 14./psp;
1097  }
1098  else if( pLab > pMax )
1099  {
1100  G4double ld = std::log(pLab) - minLogP;
1101  G4double ld2 = ld*ld;
1102  fElasticXsc = cofLogE*ld2 + 2.23;
1103  fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
1104  }
1105  else
1106  {
1107  G4double ld = std::log(pLab) - minLogP;
1108  G4double ld2 = ld*ld;
1109  G4double sp = std::sqrt(pLab);
1110  G4double psp = pLab*sp;
1111  G4double p2 = pLab*pLab;
1112  G4double p4 = p2*p2;
1113 
1114  G4double lh = pLab - 0.98;
1115  G4double hd = lh*lh + .045;
1116 
1117 
1118  fElasticXsc = 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) // + .004/md
1119  + .15/hd;
1120  fTotalXsc = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
1121  // + .006/md + 0.01/hd1 + 0.02/hd2
1122  + .60/hd;
1123  }
1124  }
1125  else if( (theParticle == theKMinus || theParticle == theK0S) && neutron ) // Kmn/K0n /////////////////////////////
1126  {
1127  if( pLab > pMax )
1128  {
1129  G4double ld = std::log(pLab) - minLogP;
1130  G4double ld2 = ld*ld;
1131  fElasticXsc = cofLogE*ld2 + 2.23;
1132  fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
1133  }
1134  else
1135  {
1136 
1137  G4double lh = pLab - 0.98;
1138  G4double hd = lh*lh + .045;
1139 
1140  fElasticXsc = // 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md
1141  5.0 + 8.1*std::pow(pLab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab + .15/hd;
1142  fTotalXsc = // 14./psp +
1143  // (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
1144  25.2 + 0. *std::pow(pLab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab
1145  // + .006/md + 0.01/hd1+ 0.02/hd2
1146  + 0.60/hd ;
1147  }
1148  }
1149  else if( (theParticle == theKPlus || theParticle == theK0L) && proton ) // Kpp/aKp //////////////////////
1150  {
1151  if( pLab < pMin )
1152  {
1153  G4double lr = pLab - .38;
1154  G4double lm = pLab - 1.;
1155  G4double md = lm*lm + .392;
1156  fElasticXsc = .7/(lr*lr + .076) + 2./md;
1157  fTotalXsc = // .7/(lr*lr + .076) +
1158  2.6/md;
1159  }
1160  else if( pLab > pMax )
1161  {
1162  G4double ld = std::log(pLab) - minLogP;
1163  G4double ld2 = ld*ld;
1164  fElasticXsc = cofLogE*ld2 + 2.23;
1165  fTotalXsc = cofLogT*ld2 + 19.2;
1166  }
1167  else
1168  {
1169  G4double ld = std::log(pLab) - minLogP;
1170  G4double ld2 = ld*ld;
1171  G4double lr = pLab - .38;
1172  G4double LE = .7/(lr*lr + .076);
1173  G4double sp = std::sqrt(pLab);
1174  G4double p2 = pLab*pLab;
1175  G4double p4 = p2*p2;
1176  G4double lm = pLab - 0.8;
1177  G4double md = lm*lm + .652;
1178  fElasticXsc = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
1179  fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 7.6/md; // + LE;
1180  }
1181  }
1182  else if( (theParticle == theKPlus || theParticle == theK0L) && neutron ) // Kpn/aKn //////////////////////////////////
1183  {
1184  if( pLab < pMin )
1185  {
1186  G4double lm = pLab - 0.94;
1187  G4double md = lm*lm + .392;
1188  fElasticXsc = 2./md;
1189  fTotalXsc = 4.6/md;
1190  }
1191  else if( pLab > pMax )
1192  {
1193  G4double ld = std::log(pLab) - minLogP;
1194  G4double ld2 = ld*ld;
1195  fElasticXsc = cofLogE*ld2 + 2.23;
1196  fTotalXsc = cofLogT*ld2 + 19.2;
1197  }
1198  else
1199  {
1200  G4double ld = std::log(pLab) - minLogP;
1201  G4double ld2 = ld*ld;
1202  G4double sp = std::sqrt(pLab);
1203  G4double p2 = pLab*pLab;
1204  G4double p4 = p2*p2;
1205  G4double lm = pLab - 0.8;
1206  G4double md = lm*lm + .652;
1207  fElasticXsc = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
1208  fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 7.6/md;
1209  }
1210  }
1211  fTotalXsc *= millibarn; // parametrised in mb
1212  fElasticXsc *= millibarn; // parametrised in mb
1213 
1214  if( proton && aParticle->GetDefinition()->GetPDGCharge() > 0. )
1215  {
1216  G4double cB = GetCoulombBarrier(aParticle, nucleon);
1217  fTotalXsc *= cB;
1218  fElasticXsc *= cB;
1219  }
1220  fInelasticXsc = fTotalXsc - fElasticXsc;
1221  if( fInelasticXsc < 0. ) fInelasticXsc = 0.;
1222 
1223  // G4cout<<fTotalXsc/millibarn<<"; "<<fElasticXsc/millibarn<<"; "<<fInelasticXsc/millibarn<<G4endl;
1224 
1225  return fTotalXsc;
1226 }
G4HadronNucleonXsc::GetCoulombBarrier
G4double GetCoulombBarrier(const G4DynamicParticle *aParticle, const G4ParticleDefinition *nucleon)
Definition: G4HadronNucleonXsc.cc:1413
neutron
Definition: G4DataQuestionaire.hh:35
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
python.hepunit.millibarn
int millibarn
Definition: hepunit.py:40
G4bool
bool G4bool
Definition: G4Types.hh:79
LE
Definition: Evaluator.cc:66
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
CLHEP::Hep3Vector::mag
double mag() const
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4double G4HadronNucleonXsc::GetKmNeutronTotXscVector ( G4double  logEnergy)
inline

Definition at line 99 of file G4HadronNucleonXsc.hh.

References G4PhysicsVector::Value().

Referenced by G4GlauberGribovCrossSection::GetKaonNucleonXscVector(), and G4ComponentGGHadronNucleusXsc::GetKaonNucleonXscVector().

99 { return fKmNeutronTotXscVector.Value(logEnergy); };
G4PhysicsVector::Value
G4double Value(G4double theEnergy, size_t &lastidx) const
Definition: G4PhysicsVector.cc:511
G4double G4HadronNucleonXsc::GetKmProtonTotXscVector ( G4double  logEnergy)
inline

Definition at line 98 of file G4HadronNucleonXsc.hh.

References G4PhysicsVector::Value().

Referenced by G4GlauberGribovCrossSection::GetKaonNucleonXscVector(), and G4ComponentGGHadronNucleusXsc::GetKaonNucleonXscVector().

98 { return fKmProtonTotXscVector.Value(logEnergy); };
G4PhysicsVector::Value
G4double Value(G4double theEnergy, size_t &lastidx) const
Definition: G4PhysicsVector.cc:511
G4double G4HadronNucleonXsc::GetKpNeutronTotXscVector ( G4double  logEnergy)
inline

Definition at line 97 of file G4HadronNucleonXsc.hh.

References G4PhysicsVector::Value().

Referenced by G4GlauberGribovCrossSection::GetKaonNucleonXscVector(), and G4ComponentGGHadronNucleusXsc::GetKaonNucleonXscVector().

97 { return fKpNeutronTotXscVector.Value(logEnergy); };
G4PhysicsVector::Value
G4double Value(G4double theEnergy, size_t &lastidx) const
Definition: G4PhysicsVector.cc:511
G4double G4HadronNucleonXsc::GetKpProtonTotXscVector ( G4double  logEnergy)
inline

Definition at line 96 of file G4HadronNucleonXsc.hh.

References G4PhysicsVector::Value().

Referenced by G4GlauberGribovCrossSection::GetKaonNucleonXscVector(), and G4ComponentGGHadronNucleusXsc::GetKaonNucleonXscVector().

96 { return fKpProtonTotXscVector.Value(logEnergy); };
G4PhysicsVector::Value
G4double Value(G4double theEnergy, size_t &lastidx) const
Definition: G4PhysicsVector.cc:511
G4double G4HadronNucleonXsc::GetTotalHadronNucleonXsc ( )
inline

Definition at line 90 of file G4HadronNucleonXsc.hh.

90 { return fTotalXsc; };
void G4HadronNucleonXsc::InitialiseKaonNucleonTotXsc ( )

Definition at line 1465 of file G4HadronNucleonXsc.cc.

References python.hepunit::millibarn, G4LPhysicsFreeVector::PutValues(), and G4PhysicsVector::SetSpline().

Referenced by G4HadronNucleonXsc().

1466 {
1467  G4int i = 0, iMax;
1468  G4double tmpxsc[106];
1469 
1470  // Kp-proton tot xsc
1471 
1472  iMax = 66;
1473  fKpProtonTotXscVector = G4LPhysicsFreeVector(iMax, fKpProtonTotTkin[0], fKpProtonTotTkin[iMax-1]);
1474  fKpProtonTotXscVector.SetSpline(true);
1475 
1476  for( i = 0; i < iMax; i++)
1477  {
1478  tmpxsc[i] = 0.;
1479 
1480  if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKpProtonTotXsc[i];
1481  else tmpxsc[i] = 0.5*(fKpProtonTotXsc[i-1]+fKpProtonTotXsc[i+1]);
1482 
1483  fKpProtonTotXscVector.PutValues(size_t(i), fKpProtonTotTkin[i], tmpxsc[i]*millibarn);
1484  }
1485 
1486  // Kp-neutron tot xsc
1487 
1488  iMax = 75;
1489  fKpNeutronTotXscVector = G4LPhysicsFreeVector(iMax, fKpNeutronTotTkin[0], fKpNeutronTotTkin[iMax-1]);
1490  fKpNeutronTotXscVector.SetSpline(true);
1491 
1492  for( i = 0; i < iMax; i++)
1493  {
1494  tmpxsc[i] = 0.;
1495  if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKpNeutronTotXsc[i];
1496  else tmpxsc[i] = 0.5*(fKpNeutronTotXsc[i-1]+fKpNeutronTotXsc[i+1]);
1497 
1498  fKpNeutronTotXscVector.PutValues(size_t(i), fKpNeutronTotTkin[i], tmpxsc[i]*millibarn);
1499  }
1500 
1501  // Km-proton tot xsc
1502 
1503  iMax = 106;
1504  fKmProtonTotXscVector = G4LPhysicsFreeVector(iMax, fKmProtonTotTkin[0], fKmProtonTotTkin[iMax-1]);
1505  fKmProtonTotXscVector.SetSpline(true);
1506 
1507  for( i = 0; i < iMax; i++)
1508  {
1509  tmpxsc[i] = 0.;
1510 
1511  if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKmProtonTotXsc[i];
1512  else tmpxsc[i] = 0.5*(fKmProtonTotXsc[i-1]+fKmProtonTotXsc[i+1]);
1513 
1514  fKmProtonTotXscVector.PutValues(size_t(i), fKmProtonTotTkin[i], tmpxsc[i]*millibarn);
1515  }
1516 
1517  // Km-neutron tot xsc
1518 
1519  iMax = 68;
1520  fKmNeutronTotXscVector = G4LPhysicsFreeVector(iMax, fKmNeutronTotTkin[0], fKmNeutronTotTkin[iMax-1]);
1521  fKmNeutronTotXscVector.SetSpline(true);
1522 
1523  for( i = 0; i < iMax; i++)
1524  {
1525  tmpxsc[i] = 0.;
1526  if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKmNeutronTotXsc[i];
1527  else tmpxsc[i] = 0.5*(fKmNeutronTotXsc[i-1]+fKmNeutronTotXsc[i+1]);
1528 
1529  fKmNeutronTotXscVector.PutValues(size_t(i), fKmNeutronTotTkin[i], tmpxsc[i]*millibarn);
1530  }
1531 }
G4LPhysicsFreeVector::PutValues
void PutValues(size_t binNumber, G4double binValue, G4double dataValue)
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.millibarn
int millibarn
Definition: hepunit.py:40
G4PhysicsVector::SetSpline
void SetSpline(G4bool)
G4double
double G4double
Definition: G4Types.hh:76
G4bool G4HadronNucleonXsc::IsApplicable ( const G4DynamicParticle *  aDP,
const G4Element *  anElement 
)
virtual

Definition at line 95 of file G4HadronNucleonXsc.cc.

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

97 {
98  G4int Z = G4lrint(anElement->GetZ());
99  G4int A = G4lrint(anElement->GetN());
100  return IsIsoApplicable(aDP, Z, A);
101 }
G4Element::GetN
G4double GetN() const
Definition: G4Element.hh:134
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4int
int G4int
Definition: G4Types.hh:78
G4HadronNucleonXsc::IsIsoApplicable
virtual G4bool IsIsoApplicable(const G4DynamicParticle *aDP, G4int Z, G4int A)
Definition: G4HadronNucleonXsc.cc:107
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
G4bool G4HadronNucleonXsc::IsIsoApplicable ( const G4DynamicParticle *  aDP,
G4int  Z,
G4int  A 
)
virtual

Definition at line 107 of file G4HadronNucleonXsc.cc.

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

Referenced by IsApplicable().

109 {
110  G4bool applicable = false;
111  // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
112  G4double kineticEnergy = aDP->GetKineticEnergy();
113 
114  const G4ParticleDefinition* theParticle = aDP->GetDefinition();
115 
116  if ( ( kineticEnergy >= fLowerLimit &&
117  Z > 1 && // >= He
118  ( theParticle == theAProton ||
119  theParticle == theGamma ||
120  theParticle == theKPlus ||
121  theParticle == theKMinus ||
122  theParticle == theSMinus) ) ||
123 
124  ( kineticEnergy >= 0.1*fLowerLimit &&
125  Z > 1 && // >= He
126  ( theParticle == theProton ||
127  theParticle == theNeutron ||
128  theParticle == thePiPlus ||
129  theParticle == thePiMinus ) ) ) applicable = true;
130 
131  return applicable;
132 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4double
double G4double
Definition: G4Types.hh:76
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:32 for Geant4.10 by   doxygen 1.8.7