00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "G4HadronNucleonXsc.hh"
00030
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4ParticleTable.hh"
00034 #include "G4IonTable.hh"
00035 #include "G4ParticleDefinition.hh"
00036 #include "G4HadTmpUtil.hh"
00037
00038
00039 G4HadronNucleonXsc::G4HadronNucleonXsc()
00040 : fUpperLimit( 10000 * GeV ),
00041 fLowerLimit( 0.03 * MeV ),
00042 fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fHadronNucleonXsc(0.0)
00043 {
00044 theGamma = G4Gamma::Gamma();
00045 theProton = G4Proton::Proton();
00046 theNeutron = G4Neutron::Neutron();
00047 theAProton = G4AntiProton::AntiProton();
00048 theANeutron = G4AntiNeutron::AntiNeutron();
00049 thePiPlus = G4PionPlus::PionPlus();
00050 thePiMinus = G4PionMinus::PionMinus();
00051 thePiZero = G4PionZero::PionZero();
00052 theKPlus = G4KaonPlus::KaonPlus();
00053 theKMinus = G4KaonMinus::KaonMinus();
00054 theK0S = G4KaonZeroShort::KaonZeroShort();
00055 theK0L = G4KaonZeroLong::KaonZeroLong();
00056 theL = G4Lambda::Lambda();
00057 theAntiL = G4AntiLambda::AntiLambda();
00058 theSPlus = G4SigmaPlus::SigmaPlus();
00059 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
00060 theSMinus = G4SigmaMinus::SigmaMinus();
00061 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
00062 theS0 = G4SigmaZero::SigmaZero();
00063 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
00064 theXiMinus = G4XiMinus::XiMinus();
00065 theXi0 = G4XiZero::XiZero();
00066 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
00067 theAXi0 = G4AntiXiZero::AntiXiZero();
00068 theOmega = G4OmegaMinus::OmegaMinus();
00069 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
00070 theD = G4Deuteron::Deuteron();
00071 theT = G4Triton::Triton();
00072 theA = G4Alpha::Alpha();
00073 theHe3 = G4He3::He3();
00074
00075 InitialiseKaonNucleonTotXsc();
00076 }
00077
00078
00079 G4HadronNucleonXsc::~G4HadronNucleonXsc()
00080 {}
00081
00082 void G4HadronNucleonXsc::CrossSectionDescription(std::ostream& outFile) const
00083 {
00084 outFile << "G4HadronNucleonXsc calculates the total, inelastic and elastic\n"
00085 << "hadron-nucleon cross sections using several different\n"
00086 << "parameterizations within the Glauber-Gribov approach. It is\n"
00087 << "valid for all incident gammas and long-lived hadrons at\n"
00088 << "energies above 30 keV. This is a cross section component which\n"
00089 << "is to be used to build a cross section data set.\n";
00090 }
00091
00092 G4bool
00093 G4HadronNucleonXsc::IsApplicable(const G4DynamicParticle* aDP,
00094 const G4Element* anElement)
00095 {
00096 G4int Z = G4lrint(anElement->GetZ());
00097 G4int A = G4lrint(anElement->GetN());
00098 return IsIsoApplicable(aDP, Z, A);
00099 }
00100
00102
00103
00104 G4bool
00105 G4HadronNucleonXsc::IsIsoApplicable(const G4DynamicParticle* aDP,
00106 G4int Z, G4int)
00107 {
00108 G4bool applicable = false;
00109
00110 G4double kineticEnergy = aDP->GetKineticEnergy();
00111
00112 const G4ParticleDefinition* theParticle = aDP->GetDefinition();
00113
00114 if ( ( kineticEnergy >= fLowerLimit &&
00115 Z > 1 &&
00116 ( theParticle == theAProton ||
00117 theParticle == theGamma ||
00118 theParticle == theKPlus ||
00119 theParticle == theKMinus ||
00120 theParticle == theSMinus) ) ||
00121
00122 ( kineticEnergy >= 0.1*fLowerLimit &&
00123 Z > 1 &&
00124 ( theParticle == theProton ||
00125 theParticle == theNeutron ||
00126 theParticle == thePiPlus ||
00127 theParticle == thePiMinus ) ) ) applicable = true;
00128
00129 return applicable;
00130 }
00131
00132
00134
00135
00136
00137
00138
00139
00140 G4double
00141 G4HadronNucleonXsc::GetHadronNucleonXscEL(const G4DynamicParticle* aParticle,
00142 const G4ParticleDefinition* nucleon )
00143 {
00144 G4double xsection;
00145
00146
00147 G4double targ_mass = 0.939*GeV;
00148
00149 G4double proj_mass = aParticle->GetMass();
00150 G4double proj_momentum = aParticle->GetMomentum().mag();
00151 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00152
00153 sMand /= GeV*GeV;
00154 proj_momentum /= GeV;
00155
00156 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00157
00158 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
00159
00160
00161 if(theParticle == theGamma && pORn )
00162 {
00163 xsection = (0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
00164 }
00165 else if(theParticle == theNeutron && pORn )
00166 {
00167 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00168 }
00169 else if(theParticle == theProton && pORn )
00170 {
00171 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00172
00173
00174
00175 }
00176 else if(theParticle == theAProton && pORn )
00177 {
00178 xsection = ( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
00179 }
00180 else if(theParticle == thePiPlus && pORn )
00181 {
00182 xsection = (13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
00183 }
00184 else if(theParticle == thePiMinus && pORn )
00185 {
00186
00187 xsection = (13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
00188 }
00189 else if(theParticle == theKPlus && pORn )
00190 {
00191 xsection = (11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
00192 }
00193 else if(theParticle == theKMinus && pORn )
00194 {
00195 xsection = (11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
00196 }
00197 else
00198 {
00199 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00200 }
00201 xsection *= millibarn;
00202
00203 fTotalXsc = xsection;
00204 fInelasticXsc = 0.83*xsection;
00205 fElasticXsc = fTotalXsc - fInelasticXsc;
00206 if (fElasticXsc < 0.)fElasticXsc = 0.;
00207
00208 return xsection;
00209 }
00210
00211
00213
00214
00215
00216
00217
00218 G4double
00219 G4HadronNucleonXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
00220 const G4ParticleDefinition* nucleon )
00221 {
00222 G4double xsection(0);
00223 G4int Zt=1, Nt=1, At=1;
00224
00225 G4double targ_mass = 0.939*GeV;
00226
00227 G4double proj_mass = aParticle->GetMass();
00228 G4double proj_momentum = aParticle->GetMomentum().mag();
00229
00230 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00231
00232 sMand /= GeV*GeV;
00233
00234
00235
00236 G4double s0 = 5.38*5.38;
00237 G4double eta1 = 0.458;
00238 G4double eta2 = 0.458;
00239 G4double B = 0.308;
00240
00241 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00242
00243 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
00244 G4bool proton = (nucleon == theProton);
00245 G4bool neutron = (nucleon == theNeutron);
00246
00247 if(theParticle == theNeutron)
00248 {
00249 if ( proton )
00250 {
00251 xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00252 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00253 }
00254 if ( neutron )
00255 {
00256 xsection = Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
00257 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00258 }
00259 }
00260 else if(theParticle == theProton)
00261 {
00262 if ( proton )
00263 {
00264 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
00265 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00266 }
00267 if ( neutron )
00268 {
00269 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00270 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00271 }
00272 }
00273 else if(theParticle == theAProton)
00274 {
00275 if ( proton )
00276 {
00277 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
00278 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
00279 }
00280 if ( neutron )
00281 {
00282 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00283 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
00284 }
00285 }
00286 else if(theParticle == thePiPlus && pORn )
00287 {
00288 xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
00289 + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
00290 }
00291 else if(theParticle == thePiMinus && pORn )
00292 {
00293 xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
00294 + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
00295 }
00296 else if(theParticle == theKPlus)
00297 {
00298 if ( proton )
00299 {
00300 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
00301 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
00302 }
00303 if ( neutron )
00304 {
00305 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
00306 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
00307 }
00308 }
00309 else if(theParticle == theKMinus)
00310 {
00311 if ( proton )
00312 {
00313 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
00314 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
00315 }
00316 if ( neutron )
00317 {
00318 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
00319 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2) );
00320 }
00321 }
00322 else if(theParticle == theSMinus && pORn )
00323 {
00324 xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
00325 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2) );
00326 }
00327 else if(theParticle == theGamma && pORn )
00328 {
00329 xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
00330 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2) );
00331
00332 }
00333 else
00334 {
00335 if ( proton )
00336 {
00337 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
00338 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2) );
00339 }
00340 if ( neutron )
00341 {
00342 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00343 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00344 }
00345 }
00346 xsection *= millibarn;
00347
00348 fTotalXsc = xsection;
00349 fInelasticXsc = 0.75*xsection;
00350 fElasticXsc = fTotalXsc - fInelasticXsc;
00351 if (fElasticXsc < 0.) fElasticXsc = 0.;
00352
00353 return xsection;
00354 }
00355
00356
00358
00359
00360
00361
00362 G4double
00363 G4HadronNucleonXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
00364 const G4ParticleDefinition* nucleon )
00365 {
00366 G4double xsection(0);
00367
00368 G4double A0, B0;
00369 G4double hpXsc(0);
00370 G4double hnXsc(0);
00371
00372
00373 G4double tM = 0.939*GeV;
00374
00375 G4double pM = aParticle->GetMass();
00376 G4double pE = aParticle->GetTotalEnergy();
00377 G4double pLab = aParticle->GetMomentum().mag();
00378
00379 G4double sMand = CalcMandelstamS ( pM , tM , pLab );
00380
00381 sMand /= GeV*GeV;
00382 pLab /= GeV;
00383 pE /= GeV;
00384 pM /= GeV;
00385
00386 G4double logP = std::log(pLab);
00387
00388
00389
00390
00391 G4double s0 = 5.38*5.38;
00392 G4double eta1 = 0.458;
00393 G4double eta2 = 0.458;
00394 G4double B = 0.308;
00395
00396
00397 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00398
00399 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
00400 G4bool proton = (nucleon == theProton);
00401 G4bool neutron = (nucleon == theNeutron);
00402
00403 if( theParticle == theNeutron && pORn )
00404 {
00405 if( pLab >= 373.)
00406 {
00407 xsection = GetHadronNucleonXscPDG(aParticle, nucleon)/millibarn;
00408
00409 fElasticXsc = 6.5 + 0.308*std::pow(std::log(sMand/400.),1.65) + 9.19*std::pow(sMand,-0.458);
00410
00411 fTotalXsc = xsection;
00412
00413 }
00414 else if( pLab >= 100.)
00415 {
00416 B0 = 7.5;
00417 A0 = 100. - B0*std::log(3.0e7);
00418
00419 xsection = A0 + B0*std::log(pE) - 11
00420
00421 + 103*std::pow(sMand,-0.165);
00422
00423 fElasticXsc = 5.53 + 0.308*std::pow(std::log(sMand/28.9),1.1) + 9.19*std::pow(sMand,-0.458);
00424
00425 fTotalXsc = xsection;
00426 }
00427 else if( pLab >= 10.)
00428 {
00429 B0 = 7.5;
00430 A0 = 100. - B0*std::log(3.0e7);
00431
00432 xsection = A0 + B0*std::log(pE) - 11
00433 + 103*std::pow(2*0.93827*pE + pM*pM+
00434 0.93827*0.93827,-0.165);
00435 fTotalXsc = xsection;
00436 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
00437 }
00438 else
00439 {
00440 if( neutron )
00441 {
00442 if( pLab < 0.4 )
00443 {
00444 hnXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
00445 fElasticXsc = hnXsc;
00446 }
00447 else if( pLab < 0.73 )
00448 {
00449 hnXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
00450 fElasticXsc = hnXsc;
00451 }
00452 else if( pLab < 1.05 )
00453 {
00454 hnXsc = 23 + 40*(std::log(pLab/0.73))*
00455 (std::log(pLab/0.73));
00456 fElasticXsc = 23 + 20*(std::log(pLab/0.73))*
00457 (std::log(pLab/0.73));
00458 }
00459 else
00460 {
00461 hnXsc = 39.0+75*(pLab - 1.2)/(std::pow(pLab,3.0) + 0.15);
00462
00463 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
00464 }
00465 fTotalXsc = hnXsc;
00466 }
00467 if( proton )
00468 {
00469 if( pLab < 0.02 )
00470 {
00471 hpXsc = 4100+30*std::pow(std::log(1.3/pLab),3.6);
00472 fElasticXsc = hpXsc;
00473 }
00474 else if( pLab < 0.8 )
00475 {
00476 hpXsc = 33+30*std::pow(std::log(pLab/1.3),4.0);
00477 fElasticXsc = hpXsc;
00478 }
00479 else if( pLab < 1.05 )
00480 {
00481 hpXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
00482 fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
00483 }
00484 else if( pLab < 1.4 )
00485 {
00486 hpXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
00487 fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
00488 }
00489 else
00490 {
00491 hpXsc = 33.3 + 20.8*(std::pow(pLab,2.0) - 1.35)/(std::pow(pLab,2.50) + 0.95);
00492
00493 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
00494 }
00495 fTotalXsc = hpXsc;
00496 }
00497 }
00498 }
00499 else if( theParticle == theProton && pORn )
00500 {
00501 if( pLab >= 373.)
00502 {
00503 xsection = GetHadronNucleonXscPDG(aParticle, nucleon)/millibarn;
00504
00505 fElasticXsc = 6.5 + 0.308*std::pow(std::log(sMand/400.),1.65) + 9.19*std::pow(sMand,-0.458);
00506
00507 fTotalXsc = xsection;
00508 }
00509 else if( pLab >= 100.)
00510 {
00511 B0 = 7.5;
00512 A0 = 100. - B0*std::log(3.0e7);
00513
00514 xsection = A0 + B0*std::log(pE) - 11 + 103*std::pow(sMand,-0.165);
00515
00516 fElasticXsc = 5.53 + 0.308*std::pow(std::log(sMand/28.9),1.1) + 9.19*std::pow(sMand,-0.458);
00517
00518 fTotalXsc = xsection;
00519 }
00520 else if( pLab >= 10.)
00521 {
00522 B0 = 7.5;
00523 A0 = 100. - B0*std::log(3.0e7);
00524
00525 xsection = A0 + B0*std::log(pE) - 11 + 103*std::pow(sMand,-0.165);
00526
00527 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
00528
00529 fTotalXsc = xsection;
00530 }
00531 else
00532 {
00533
00534
00535 if( proton )
00536 {
00537 if( pLab < 0.4 )
00538 {
00539 hpXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
00540 fElasticXsc = hpXsc;
00541 }
00542 else if( pLab < 0.73 )
00543 {
00544 hpXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
00545 fElasticXsc = hpXsc;
00546 }
00547 else if( pLab < 1.05 )
00548 {
00549 hpXsc = 23 + 40*(std::log(pLab/0.73))*
00550 (std::log(pLab/0.73));
00551 fElasticXsc = 23 + 20*(std::log(pLab/0.73))*
00552 (std::log(pLab/0.73));
00553 }
00554 else
00555 {
00556 hpXsc = 39.0+75*(pLab - 1.2)/(std::pow(pLab,3.0) + 0.15);
00557
00558 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
00559 }
00560 fTotalXsc = hpXsc;
00561 }
00562 if( neutron )
00563 {
00564 if( pLab < 0.02 )
00565 {
00566 hnXsc = 4100+30*std::pow(std::log(1.3/pLab),3.6);
00567 fElasticXsc = hnXsc;
00568 }
00569 else if( pLab < 0.8 )
00570 {
00571 hnXsc = 33+30*std::pow(std::log(pLab/1.3),4.0);
00572 fElasticXsc = hnXsc;
00573 }
00574 else if( pLab < 1.05 )
00575 {
00576 hnXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
00577 fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
00578 }
00579 else if( pLab < 1.4 )
00580 {
00581 hnXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
00582 fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
00583 }
00584 else
00585 {
00586 hnXsc = 33.3 + 20.8*(std::pow(pLab,2.0) - 1.35)/(std::pow(pLab,2.50) + 0.95);
00587
00588 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
00589 }
00590 fTotalXsc = hnXsc;
00591 }
00592 }
00593 }
00594 else if( theParticle == theAProton && pORn )
00595 {
00596 if( proton )
00597 {
00598 xsection = 35.45 + B*std::pow(std::log(sMand/s0),2.)
00599 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2);
00600 }
00601 if( neutron )
00602 {
00603 xsection = 35.80 + B*std::pow(std::log(sMand/s0),2.)
00604 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2);
00605 }
00606 fTotalXsc = xsection;
00607 }
00608 else if( theParticle == thePiPlus && pORn )
00609 {
00610 if( proton )
00611 {
00612 if( pLab < 0.28 )
00613 {
00614 hpXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
00615 fElasticXsc = hpXsc;
00616 }
00617 else if( pLab < 0.4 )
00618 {
00619 hpXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
00620 fElasticXsc = hpXsc;
00621 }
00622 else if( pLab < 0.68 )
00623 {
00624 hpXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
00625 fElasticXsc = hpXsc;
00626 }
00627 else if( pLab < 0.85 )
00628 {
00629 G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
00630 hpXsc = Ex4 + 14.9;
00631 fElasticXsc = hpXsc*std::exp(-3.*(pLab - 0.68));
00632 }
00633 else if( pLab < 1.15 )
00634 {
00635 G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
00636 hpXsc = Ex4 + 14.9;
00637
00638 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
00639 }
00640 else if( pLab < 1.4)
00641 {
00642 G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
00643 G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
00644 hpXsc = Ex1 + Ex2 + 27.5;
00645 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
00646 }
00647 else if( pLab < 2.0 )
00648 {
00649 G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
00650 G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
00651 hpXsc = Ex1 + Ex2 + 27.5;
00652 fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
00653 }
00654 else if( pLab < 3.5 )
00655 {
00656 G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
00657 G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
00658 hpXsc = Ex1 + Ex2 + 27.5;
00659 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
00660 }
00661 else if( pLab < 200. )
00662 {
00663 hpXsc = 10.6 + 2.*std::log(pE) + 25*std::pow(pE, -0.43 );
00664
00665 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
00666 }
00667 else
00668 {
00669 hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
00670 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
00671 }
00672 fTotalXsc = hpXsc;
00673 }
00674 if( neutron )
00675 {
00676 if( pLab < 0.28 )
00677 {
00678 hnXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
00679 fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
00680 }
00681 else if( pLab < 0.395676 )
00682 {
00683 hnXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
00684 fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
00685 }
00686 else if( pLab < 0.5 )
00687 {
00688 hnXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
00689 fElasticXsc = 0.37*hnXsc;
00690 }
00691 else if( pLab < 0.65 )
00692 {
00693 hnXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
00694 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
00695 }
00696 else if( pLab < 0.72 )
00697 {
00698 hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00699 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00700 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
00701 }
00702 else if( pLab < 0.88 )
00703 {
00704 hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00705 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00706 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
00707 }
00708 else if( pLab < 1.03 )
00709 {
00710 hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00711 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00712 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
00713 }
00714 else if( pLab < 1.15 )
00715 {
00716 hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00717 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00718 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
00719 }
00720 else if( pLab < 1.3 )
00721 {
00722 hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00723 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00724 fElasticXsc = 3. + 13./pLab;
00725 }
00726 else if( pLab < 2.6 )
00727 {
00728 hnXsc = 36.1 + 0.079-4.313*std::log(pLab)+
00729 3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
00730 1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
00731 fElasticXsc = 3. + 13./pLab;
00732 }
00733 else if( pLab < 20. )
00734 {
00735 hnXsc = 36.1 + 0.079 - 4.313*std::log(pLab)+
00736 3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
00737 1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
00738 fElasticXsc = 3. + 13./pLab;
00739 }
00740 else
00741 {
00742 hnXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
00743 fElasticXsc = 3. + 13./pLab;
00744 }
00745 fTotalXsc = hnXsc;
00746 }
00747 }
00748 else if( theParticle == thePiMinus && pORn )
00749 {
00750 if( neutron )
00751 {
00752 if( pLab < 0.28 )
00753 {
00754 hnXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
00755 fElasticXsc = hnXsc;
00756 }
00757 else if( pLab < 0.4 )
00758 {
00759 hnXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
00760 fElasticXsc = hnXsc;
00761 }
00762 else if( pLab < 0.68 )
00763 {
00764 hnXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
00765 fElasticXsc = hnXsc;
00766 }
00767 else if( pLab < 0.85 )
00768 {
00769 G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
00770 hnXsc = Ex4 + 14.9;
00771 fElasticXsc = hnXsc*std::exp(-3.*(pLab - 0.68));
00772 }
00773 else if( pLab < 1.15 )
00774 {
00775 G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
00776 hnXsc = Ex4 + 14.9;
00777
00778 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
00779 }
00780 else if( pLab < 1.4)
00781 {
00782 G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
00783 G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
00784 hnXsc = Ex1 + Ex2 + 27.5;
00785 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
00786 }
00787 else if( pLab < 2.0 )
00788 {
00789 G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
00790 G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
00791 hnXsc = Ex1 + Ex2 + 27.5;
00792 fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
00793 }
00794 else if( pLab < 3.5 )
00795 {
00796 G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
00797 G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
00798 hnXsc = Ex1 + Ex2 + 27.5;
00799 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
00800 }
00801 else if( pLab < 200. )
00802 {
00803 hnXsc = 10.6 + 2.*std::log(pE) + 25*std::pow(pE, -0.43 );
00804 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
00805 }
00806 else
00807 {
00808 hnXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
00809 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
00810 }
00811 fTotalXsc = hnXsc;
00812 }
00813 if( proton )
00814 {
00815 if( pLab < 0.28 )
00816 {
00817 hpXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
00818 fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
00819 }
00820 else if( pLab < 0.395676 )
00821 {
00822 hpXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
00823 fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
00824 }
00825 else if( pLab < 0.5 )
00826 {
00827 hpXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
00828 fElasticXsc = 0.37*hpXsc;
00829 }
00830 else if( pLab < 0.65 )
00831 {
00832 hpXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
00833 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
00834 }
00835 else if( pLab < 0.72 )
00836 {
00837 hpXsc = 36.1+
00838 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00839 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00840 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
00841 }
00842 else if( pLab < 0.88 )
00843 {
00844 hpXsc = 36.1+
00845 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00846 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00847 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
00848 }
00849 else if( pLab < 1.03 )
00850 {
00851 hpXsc = 36.1+
00852 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00853 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00854 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
00855 }
00856 else if( pLab < 1.15 )
00857 {
00858 hpXsc = 36.1+
00859 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00860 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00861 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
00862 }
00863 else if( pLab < 1.3 )
00864 {
00865 hpXsc = 36.1+
00866 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
00867 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
00868 fElasticXsc = 3. + 13./pLab;
00869 }
00870 else if( pLab < 2.6 )
00871 {
00872 hpXsc = 36.1+0.079-4.313*std::log(pLab)+
00873 3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
00874 1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
00875 fElasticXsc = 3. +13./pLab;
00876 }
00877 else
00878 {
00879 hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
00880 fElasticXsc = 3. + 13./pLab;
00881 }
00882 fTotalXsc = hpXsc;
00883 }
00884 }
00885 else if( theParticle == theKPlus && pORn )
00886 {
00887 if( proton )
00888 {
00889 xsection = 17.91 + B*std::pow(std::log(sMand/s0),2.)
00890 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2);
00891 }
00892 if( neutron )
00893 {
00894 xsection = 17.87 + B*std::pow(std::log(sMand/s0),2.)
00895 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2);
00896 }
00897 fTotalXsc = xsection;
00898 }
00899 else if( theParticle == theKMinus && pORn )
00900 {
00901 if( proton )
00902 {
00903 xsection = 17.91 + B*std::pow(std::log(sMand/s0),2.)
00904 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2);
00905 }
00906 if( neutron )
00907 {
00908 xsection = 17.87 + B*std::pow(std::log(sMand/s0),2.)
00909 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2);
00910 }
00911 fTotalXsc = xsection;
00912 }
00913 else if( theParticle == theSMinus && pORn )
00914 {
00915 xsection = 35.20 + B*std::pow(std::log(sMand/s0),2.)
00916 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2);
00917 }
00918 else if( theParticle == theGamma && pORn )
00919 {
00920 xsection = 0.0 + B*std::pow(std::log(sMand/s0),2.)
00921 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2);
00922 fTotalXsc = xsection;
00923 }
00924 else
00925 {
00926 if( proton )
00927 {
00928 xsection = 35.45 + B*std::pow(std::log(sMand/s0),2.)
00929 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2);
00930 }
00931 if( neutron )
00932 {
00933 xsection += 35.80 + B*std::pow(std::log(sMand/s0),2.)
00934 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2);
00935 }
00936 fTotalXsc = xsection;
00937 }
00938 fTotalXsc *= millibarn;
00939 fElasticXsc *= millibarn;
00940
00941 if( proton && aParticle->GetDefinition()->GetPDGCharge() > 0. )
00942 {
00943 G4double cB = GetCoulombBarrier(aParticle, nucleon);
00944 fTotalXsc *= cB;
00945 fElasticXsc *= cB;
00946 }
00947 fInelasticXsc = fTotalXsc - fElasticXsc;
00948 if( fInelasticXsc < 0. ) fInelasticXsc = 0.;
00949
00950
00951
00952 return fTotalXsc;
00953 }
00954
00956
00957
00958
00959
00960 G4double
00961 G4HadronNucleonXsc::GetHadronNucleonXscVU(const G4DynamicParticle* aParticle,
00962 const G4ParticleDefinition* nucleon )
00963 {
00964 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
00965 G4int absPDGcode = std::abs(PDGcode);
00966 G4double Elab = aParticle->GetTotalEnergy();
00967
00968 G4double Plab = aParticle->GetMomentum().mag();
00969
00970
00971 Elab /= GeV;
00972 Plab /= GeV;
00973
00974 G4double LogPlab = std::log( Plab );
00975 G4double sqrLogPlab = LogPlab * LogPlab;
00976
00977 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
00978 G4bool proton = (nucleon == theProton);
00979 G4bool neutron = (nucleon == theNeutron);
00980
00981
00982 if( absPDGcode > 1000 && pORn )
00983 {
00984 if(proton)
00985 {
00986 fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
00987 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
00988 }
00989 if(neutron)
00990 {
00991 fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
00992 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
00993 }
00994 }
00995 else if( PDGcode == 211 && pORn )
00996 {
00997 if(proton)
00998 {
00999 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
01000 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
01001 }
01002 if(neutron)
01003 {
01004 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
01005 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
01006 }
01007 }
01008 else if( PDGcode == -211 && pORn )
01009 {
01010 if(proton)
01011 {
01012 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
01013 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
01014 }
01015 if(neutron)
01016 {
01017 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
01018 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
01019 }
01020 }
01021 else if( PDGcode == 111 && pORn )
01022 {
01023 if(proton)
01024 {
01025 fTotalXsc = (16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab +
01026 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2;
01027
01028 fElasticXsc = ( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab +
01029 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2;
01030
01031 }
01032 if(neutron)
01033 {
01034 fTotalXsc = (33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab +
01035 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2;
01036 fElasticXsc = ( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab +
01037 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2;
01038 }
01039 }
01040 else if( PDGcode == 321 && pORn )
01041 {
01042 if(proton)
01043 {
01044 fTotalXsc = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
01045 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
01046 }
01047 if(neutron)
01048 {
01049 fTotalXsc = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
01050 fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
01051 }
01052 }
01053 else if( PDGcode ==-321 && pORn )
01054 {
01055 if(proton)
01056 {
01057 fTotalXsc = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66*sqrLogPlab - 5.6*LogPlab;
01058 fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29*sqrLogPlab - 2.4*LogPlab;
01059 }
01060 if(neutron)
01061 {
01062 fTotalXsc = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab;
01063 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab;
01064 }
01065 }
01066 else if( PDGcode == 311 && pORn )
01067 {
01068 if(proton)
01069 {
01070 fTotalXsc = ( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab +
01071 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2;
01072 fElasticXsc = ( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab +
01073 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2;
01074 }
01075 if(neutron)
01076 {
01077 fTotalXsc = ( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab +
01078 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2;
01079 fElasticXsc = ( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab +
01080 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2;
01081 }
01082 }
01083 else
01084 {
01085 if(proton)
01086 {
01087 fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
01088 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
01089 }
01090 if(neutron)
01091 {
01092 fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
01093 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
01094 }
01095 }
01096 fTotalXsc *= millibarn;
01097 fElasticXsc *= millibarn;
01098 fInelasticXsc = fTotalXsc - fElasticXsc;
01099 if (fInelasticXsc < 0.) fInelasticXsc = 0.;
01100
01101 return fTotalXsc;
01102 }
01103
01105
01106
01107
01108 G4double G4HadronNucleonXsc::CalculateEcmValue( const G4double mp ,
01109 const G4double mt ,
01110 const G4double Plab )
01111 {
01112 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01113 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
01114
01115
01116
01117 return Ecm ;
01118 }
01119
01120
01122
01123
01124
01125 G4double G4HadronNucleonXsc::CalcMandelstamS( const G4double mp ,
01126 const G4double mt ,
01127 const G4double Plab )
01128 {
01129 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01130 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
01131
01132 return sMand;
01133 }
01134
01135
01137
01138
01139
01140 G4double G4HadronNucleonXsc::GetCoulombBarrier(const G4DynamicParticle* aParticle,
01141 const G4ParticleDefinition* nucleon )
01142 {
01143 G4double ratio;
01144
01145 G4double tR = 0.895*fermi, pR;
01146
01147 if ( aParticle->GetDefinition() == theProton ) pR = 0.895*fermi;
01148 else if( aParticle->GetDefinition() == thePiPlus ) pR = 0.663*fermi;
01149 else if( aParticle->GetDefinition() == theKPlus ) pR = 0.340*fermi;
01150 else pR = 0.500*fermi;
01151
01152 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
01153 G4double tZ = nucleon->GetPDGCharge();
01154
01155 G4double pTkin = aParticle->GetKineticEnergy();
01156
01157 G4double pM = aParticle->GetDefinition()->GetPDGMass();
01158 G4double tM = nucleon->GetPDGMass();
01159
01160 G4double pElab = pTkin + pM;
01161
01162 G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
01163
01164 G4double totTcm = totEcm - pM -tM;
01165
01166 G4double bC = fine_structure_const*hbarc*pZ*tZ;
01167 bC /= pR + tR;
01168 bC /= 2.;
01169
01170
01171
01172
01173 if( totTcm <= bC ) ratio = 0.;
01174 else ratio = 1. - bC/totTcm;
01175
01176
01177 if( ratio < 0.) ratio = 0.;
01178
01179
01180 return ratio;
01181 }
01182
01183
01184
01185
01186
01188
01189
01190
01191
01192 void G4HadronNucleonXsc::InitialiseKaonNucleonTotXsc()
01193 {
01194 G4int i = 0, iMax;
01195 G4double tmpxsc[106];
01196
01197
01198
01199 iMax = 66;
01200 fKpProtonTotXscVector = G4LPhysicsFreeVector(iMax, fKpProtonTotTkin[0], fKpProtonTotTkin[iMax-1]);
01201 fKpProtonTotXscVector.SetSpline(true);
01202
01203 for( i = 0; i < iMax; i++)
01204 {
01205 tmpxsc[i] = 0.;
01206
01207 if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKpProtonTotXsc[i];
01208 else tmpxsc[i] = 0.5*(fKpProtonTotXsc[i-1]+fKpProtonTotXsc[i+1]);
01209
01210 fKpProtonTotXscVector.PutValues(size_t(i), fKpProtonTotTkin[i], tmpxsc[i]*millibarn);
01211 }
01212
01213
01214
01215 iMax = 75;
01216 fKpNeutronTotXscVector = G4LPhysicsFreeVector(iMax, fKpNeutronTotTkin[0], fKpNeutronTotTkin[iMax-1]);
01217 fKpNeutronTotXscVector.SetSpline(true);
01218
01219 for( i = 0; i < iMax; i++)
01220 {
01221 tmpxsc[i] = 0.;
01222 if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKpNeutronTotXsc[i];
01223 else tmpxsc[i] = 0.5*(fKpNeutronTotXsc[i-1]+fKpNeutronTotXsc[i+1]);
01224
01225 fKpNeutronTotXscVector.PutValues(size_t(i), fKpNeutronTotTkin[i], tmpxsc[i]*millibarn);
01226 }
01227
01228
01229
01230 iMax = 106;
01231 fKmProtonTotXscVector = G4LPhysicsFreeVector(iMax, fKmProtonTotTkin[0], fKmProtonTotTkin[iMax-1]);
01232 fKmProtonTotXscVector.SetSpline(true);
01233
01234 for( i = 0; i < iMax; i++)
01235 {
01236 tmpxsc[i] = 0.;
01237
01238 if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKmProtonTotXsc[i];
01239 else tmpxsc[i] = 0.5*(fKmProtonTotXsc[i-1]+fKmProtonTotXsc[i+1]);
01240
01241 fKmProtonTotXscVector.PutValues(size_t(i), fKmProtonTotTkin[i], tmpxsc[i]*millibarn);
01242 }
01243
01244
01245
01246 iMax = 68;
01247 fKmNeutronTotXscVector = G4LPhysicsFreeVector(iMax, fKmNeutronTotTkin[0], fKmNeutronTotTkin[iMax-1]);
01248 fKmNeutronTotXscVector.SetSpline(true);
01249
01250 for( i = 0; i < iMax; i++)
01251 {
01252 tmpxsc[i] = 0.;
01253 if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKmNeutronTotXsc[i];
01254 else tmpxsc[i] = 0.5*(fKmNeutronTotXsc[i-1]+fKmNeutronTotXsc[i+1]);
01255
01256 fKmNeutronTotXscVector.PutValues(size_t(i), fKmNeutronTotTkin[i], tmpxsc[i]*millibarn);
01257 }
01258 }
01259
01261
01262
01263
01264 const G4double G4HadronNucleonXsc::fKpProtonTotXsc[66] = {
01265 0.000000e+00, 1.592400e-01, 3.184700e-01, 7.961800e-01, 1.433120e+00, 2.070060e+00,
01266 2.866240e+00, 3.582800e+00, 4.378980e+00, 5.015920e+00, 5.573250e+00, 6.449040e+00,
01267 7.404460e+00, 8.200640e+00, 8.837580e+00, 9.713380e+00, 1.027070e+01, 1.090764e+01,
01268 1.130573e+01, 1.170382e+01, 1.242038e+01, 1.281847e+01, 1.321656e+01, 1.337580e+01,
01269 1.345541e+01, 1.329618e+01, 1.265924e+01, 1.242038e+01, 1.250000e+01, 1.305732e+01,
01270 1.369427e+01, 1.425159e+01, 1.544586e+01, 1.648089e+01, 1.751592e+01, 1.791401e+01,
01271 1.791401e+01, 1.775478e+01, 1.751592e+01, 1.735669e+01, 1.719745e+01, 1.711783e+01,
01272 1.703822e+01, 1.695860e+01, 1.695860e+01, 1.695860e+01, 1.695860e+01, 1.687898e+01,
01273 1.687898e+01, 1.703822e+01, 1.719745e+01, 1.735669e+01, 1.751592e+01, 1.767516e+01,
01274 1.783439e+01, 1.799363e+01, 1.815287e+01, 1.839172e+01, 1.855095e+01, 1.871019e+01,
01275 1.886943e+01, 1.918790e+01, 1.942675e+01, 1.966561e+01, 2.006369e+01, 2.054140e+01
01276 };
01277
01278
01279 const G4double G4HadronNucleonXsc::fKpProtonTotTkin[66] = {
01280 2.100000e+00, 2.180770e+00, 2.261540e+00, 2.396150e+00, 2.476920e+00, 2.557690e+00,
01281 2.557690e+00, 2.584620e+00, 2.638460e+00, 2.665380e+00, 2.719230e+00, 2.746150e+00,
01282 2.800000e+00, 2.853850e+00, 2.934620e+00, 3.042310e+00, 3.150000e+00, 3.311540e+00,
01283 3.446150e+00, 3.607690e+00, 3.930770e+00, 4.226920e+00, 4.361540e+00, 4.846150e+00,
01284 4.980770e+00, 5.088460e+00, 5.465380e+00, 5.653850e+00, 5.950000e+00, 6.084620e+00,
01285 6.246150e+00, 6.300000e+00, 6.380770e+00, 6.515380e+00, 6.730770e+00, 6.838460e+00,
01286 7.000000e+00, 7.161540e+00, 7.323080e+00, 7.457690e+00, 7.619230e+00, 7.780770e+00,
01287 7.915380e+00, 8.130770e+00, 8.265380e+00, 8.453850e+00, 8.642310e+00, 8.803850e+00,
01288 9.019230e+00, 9.234620e+00, 9.530770e+00, 9.773080e+00, 1.001538e+01, 1.017692e+01,
01289 1.033846e+01, 1.058077e+01, 1.082308e+01, 1.098462e+01, 1.114615e+01, 1.138846e+01,
01290 1.160385e+01, 1.173846e+01, 1.192692e+01, 1.216923e+01, 1.238461e+01, 1.257308e+01
01291 };
01292
01293 const G4double G4HadronNucleonXsc::fKpNeutronTotXsc[75] = {
01294 3.980900e-01, 3.184700e-01, 3.184700e-01, 3.980900e-01, 3.980900e-01, 3.980900e-01,
01295 3.980900e-01, 3.980900e-01, 3.980900e-01, 4.777100e-01, 3.980900e-01, 3.980900e-01,
01296 4.777100e-01, 5.573200e-01, 1.035030e+00, 1.512740e+00, 2.149680e+00, 2.786620e+00,
01297 3.503180e+00, 4.219750e+00, 5.015920e+00, 5.652870e+00, 6.289810e+00, 7.245220e+00,
01298 8.121020e+00, 8.837580e+00, 9.633760e+00, 1.042994e+01, 1.114650e+01, 1.194268e+01,
01299 1.265924e+01, 1.329618e+01, 1.393312e+01, 1.449045e+01, 1.496815e+01, 1.552548e+01,
01300 1.592357e+01, 1.664013e+01, 1.727707e+01, 1.783439e+01, 1.831210e+01, 1.902866e+01,
01301 1.902866e+01, 1.878981e+01, 1.847134e+01, 1.831210e+01, 1.807325e+01, 1.791401e+01,
01302 1.783439e+01, 1.767516e+01, 1.759554e+01, 1.743631e+01, 1.743631e+01, 1.751592e+01,
01303 1.743631e+01, 1.735669e+01, 1.751592e+01, 1.759554e+01, 1.767516e+01, 1.783439e+01,
01304 1.783439e+01, 1.791401e+01, 1.815287e+01, 1.823248e+01, 1.847134e+01, 1.878981e+01,
01305 1.894905e+01, 1.902866e+01, 1.934713e+01, 1.966561e+01, 1.990446e+01, 2.014331e+01,
01306 2.030255e+01, 2.046178e+01, 2.085987e+01
01307 };
01308
01309 const G4double G4HadronNucleonXsc::fKpNeutronTotTkin[75] = {
01310 2.692000e-02, 1.615400e-01, 2.961500e-01, 4.576900e-01, 6.461500e-01, 7.538500e-01,
01311 8.884600e-01, 1.103850e+00, 1.211540e+00, 1.400000e+00, 1.561540e+00, 1.776920e+00,
01312 1.992310e+00, 2.126920e+00, 2.342310e+00, 2.423080e+00, 2.557690e+00, 2.692310e+00,
01313 2.800000e+00, 2.988460e+00, 3.203850e+00, 3.365380e+00, 3.500000e+00, 3.688460e+00,
01314 3.850000e+00, 4.011540e+00, 4.173080e+00, 4.415380e+00, 4.630770e+00, 4.873080e+00,
01315 5.061540e+00, 5.276920e+00, 5.492310e+00, 5.707690e+00, 5.896150e+00, 6.030770e+00,
01316 6.138460e+00, 6.219230e+00, 6.273080e+00, 6.326920e+00, 6.407690e+00, 6.650000e+00,
01317 6.784620e+00, 7.026920e+00, 7.242310e+00, 7.350000e+00, 7.484620e+00, 7.619230e+00,
01318 7.807690e+00, 7.915380e+00, 8.050000e+00, 8.211540e+00, 8.453850e+00, 8.588460e+00,
01319 8.830770e+00, 9.073080e+00, 9.288460e+00, 9.476920e+00, 9.665380e+00, 9.826920e+00,
01320 1.004231e+01, 1.031154e+01, 1.052692e+01, 1.071538e+01, 1.095769e+01, 1.120000e+01,
01321 1.138846e+01, 1.155000e+01, 1.176538e+01, 1.190000e+01, 1.214231e+01, 1.222308e+01,
01322 1.238461e+01, 1.246538e+01, 1.265385e+01
01323 };
01324
01325 const G4double G4HadronNucleonXsc::fKmProtonTotXsc[106] = {
01326 1.136585e+02, 9.749129e+01, 9.275262e+01, 8.885017e+01, 8.334146e+01, 7.955401e+01,
01327 7.504530e+01, 7.153658e+01, 6.858537e+01, 6.674913e+01, 6.525784e+01, 6.448781e+01,
01328 6.360279e+01, 6.255401e+01, 6.127526e+01, 6.032404e+01, 5.997910e+01, 5.443554e+01,
01329 5.376307e+01, 5.236934e+01, 5.113937e+01, 5.090941e+01, 4.967944e+01, 4.844948e+01,
01330 4.705575e+01, 4.638327e+01, 4.571080e+01, 4.475958e+01, 4.397213e+01, 4.257840e+01,
01331 4.102090e+01, 4.090592e+01, 3.906969e+01, 3.839721e+01, 3.756097e+01, 3.644599e+01,
01332 3.560976e+01, 3.533101e+01, 3.533101e+01, 3.644599e+01, 3.811847e+01, 3.839721e+01,
01333 3.979094e+01, 4.090592e+01, 4.257840e+01, 4.341463e+01, 4.425087e+01, 4.564460e+01,
01334 4.759582e+01, 4.703833e+01, 4.843206e+01, 4.787457e+01, 4.452962e+01, 4.202090e+01,
01335 4.034843e+01, 3.839721e+01, 3.616725e+01, 3.365854e+01, 3.170732e+01, 3.087108e+01,
01336 3.170732e+01, 3.254355e+01, 3.310104e+01, 3.254355e+01, 3.142857e+01, 3.059233e+01,
01337 2.947735e+01, 2.891986e+01, 2.836237e+01, 2.752613e+01, 2.696864e+01, 2.641115e+01,
01338 2.501742e+01, 2.473868e+01, 2.418118e+01, 2.362369e+01, 2.334495e+01, 2.278746e+01,
01339 2.250871e+01, 2.222997e+01, 2.167247e+01, 2.139373e+01, 2.139373e+01, 2.139373e+01,
01340 2.111498e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01,
01341 2.083624e+01, 2.055749e+01, 2.055749e+01, 2.055749e+01, 2.027875e+01, 2.000000e+01,
01342 2.055749e+01, 2.027875e+01, 2.083624e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01,
01343 2.083624e+01, 2.083624e+01, 2.139373e+01, 2.139373e+01
01344 };
01345
01346 const G4double G4HadronNucleonXsc::fKmProtonTotTkin[106] = {
01347 4.017980e+00, 4.125840e+00, 4.179780e+00, 4.251690e+00, 4.287640e+00, 4.341570e+00,
01348 4.395510e+00, 4.467420e+00, 4.503370e+00, 4.575280e+00, 4.683150e+00, 4.737080e+00,
01349 4.773030e+00, 4.826970e+00, 4.880900e+00, 4.916850e+00, 4.952810e+00, 4.988760e+00,
01350 4.988760e+00, 5.006740e+00, 5.006740e+00, 5.042700e+00, 5.078650e+00, 5.114610e+00,
01351 5.132580e+00, 5.150560e+00, 5.186520e+00, 5.204490e+00, 5.276400e+00, 5.348310e+00,
01352 5.366290e+00, 5.384270e+00, 5.456180e+00, 5.564040e+00, 5.600000e+00, 5.671910e+00,
01353 5.743820e+00, 5.833710e+00, 5.905620e+00, 5.977530e+00, 6.085390e+00, 6.085390e+00,
01354 6.157300e+00, 6.175280e+00, 6.211240e+00, 6.229210e+00, 6.247190e+00, 6.337080e+00,
01355 6.391010e+00, 6.516850e+00, 6.462920e+00, 6.498880e+00, 6.570790e+00, 6.606740e+00,
01356 6.660670e+00, 6.678650e+00, 6.696630e+00, 6.732580e+00, 6.804490e+00, 6.876400e+00,
01357 6.948310e+00, 7.020220e+00, 7.074160e+00, 7.182020e+00, 7.235960e+00, 7.289890e+00,
01358 7.397750e+00, 7.523600e+00, 7.631460e+00, 7.757300e+00, 7.901120e+00, 8.062920e+00,
01359 8.260670e+00, 8.386520e+00, 8.530340e+00, 8.674160e+00, 8.817980e+00, 8.943820e+00,
01360 9.087640e+00, 9.267420e+00, 9.429210e+00, 9.573030e+00, 9.698880e+00, 9.896630e+00,
01361 1.002247e+01, 1.016629e+01, 1.031011e+01, 1.048989e+01, 1.063371e+01, 1.077753e+01,
01362 1.095730e+01, 1.108315e+01, 1.120899e+01, 1.135281e+01, 1.149663e+01, 1.162247e+01,
01363 1.174831e+01, 1.187416e+01, 1.200000e+01, 1.212584e+01, 1.221573e+01, 1.234157e+01,
01364 1.239551e+01, 1.250337e+01, 1.261124e+01, 1.273708e+01
01365 };
01366
01367 const G4double G4HadronNucleonXsc::fKmNeutronTotXsc[68] = {
01368 2.621810e+01, 2.741123e+01, 2.868413e+01, 2.963889e+01, 3.067343e+01, 3.178759e+01,
01369 3.282148e+01, 3.417466e+01, 3.536778e+01, 3.552620e+01, 3.544576e+01, 3.496756e+01,
01370 3.433030e+01, 3.401166e+01, 3.313537e+01, 3.257772e+01, 3.178105e+01, 3.138264e+01,
01371 3.074553e+01, 2.970952e+01, 2.891301e+01, 2.827542e+01, 2.787700e+01, 2.715978e+01,
01372 2.660181e+01, 2.612394e+01, 2.564574e+01, 2.516721e+01, 2.421098e+01, 2.365235e+01,
01373 2.317366e+01, 2.261437e+01, 2.237389e+01, 2.205427e+01, 2.181395e+01, 2.165357e+01,
01374 2.149335e+01, 2.133297e+01, 2.109232e+01, 2.093128e+01, 2.069030e+01, 2.052992e+01,
01375 2.028927e+01, 2.012824e+01, 1.996737e+01, 1.996590e+01, 1.988530e+01, 1.964432e+01,
01376 1.948361e+01, 1.940236e+01, 1.940040e+01, 1.931882e+01, 1.947593e+01, 1.947429e+01,
01377 1.939320e+01, 1.939157e+01, 1.946922e+01, 1.962715e+01, 1.970481e+01, 1.970301e+01,
01378 1.993958e+01, 2.009669e+01, 2.025380e+01, 2.033178e+01, 2.049003e+01, 2.064747e+01,
01379 2.080540e+01, 2.096333e+01
01380 };
01381
01382 const G4double G4HadronNucleonXsc::fKmNeutronTotTkin[68] = {
01383 5.708500e+00, 5.809560e+00, 5.896270e+00, 5.954120e+00, 5.997630e+00, 6.041160e+00,
01384 6.142160e+00, 6.171410e+00, 6.272470e+00, 6.344390e+00, 6.416230e+00, 6.459180e+00,
01385 6.487690e+00, 6.501940e+00, 6.544740e+00, 6.573280e+00, 6.616110e+00, 6.644710e+00,
01386 6.658840e+00, 6.744700e+00, 6.773150e+00, 6.830410e+00, 6.859010e+00, 6.916240e+00,
01387 6.973530e+00, 6.987730e+00, 7.030670e+00, 7.102360e+00, 7.173880e+00, 7.288660e+00,
01388 7.374720e+00, 7.547000e+00, 7.690650e+00, 7.791150e+00, 7.920420e+00, 8.020980e+00,
01389 8.107160e+00, 8.207720e+00, 8.365740e+00, 8.523790e+00, 8.710560e+00, 8.811110e+00,
01390 8.969140e+00, 9.127190e+00, 9.270860e+00, 9.400230e+00, 9.486440e+00, 9.673210e+00,
01391 9.802510e+00, 9.946220e+00, 1.011870e+01, 1.029116e+01, 1.047808e+01, 1.062181e+01,
01392 1.075114e+01, 1.089488e+01, 1.106739e+01, 1.118244e+01, 1.135496e+01, 1.151307e+01,
01393 1.171439e+01, 1.190130e+01, 1.208822e+01, 1.223199e+01, 1.231829e+01, 1.247646e+01,
01394 1.259150e+01, 1.270655e+01
01395 };
01396
01397
01398
01399
01400
01401
01402