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
00030 #include "G4GGNuclNuclCrossSection.hh"
00031
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "G4ParticleTable.hh"
00035 #include "G4IonTable.hh"
00036 #include "G4ParticleDefinition.hh"
00037 #include "G4HadTmpUtil.hh"
00038 #include "G4HadronNucleonXsc.hh"
00039
00040
00041 #include "G4CrossSectionFactory.hh"
00042
00043 G4_DECLARE_XS_FACTORY(G4GGNuclNuclCrossSection);
00044
00045
00046 G4GGNuclNuclCrossSection::G4GGNuclNuclCrossSection()
00047 : G4VCrossSectionDataSet(Default_Name()),
00048 fUpperLimit(100000*GeV), fLowerLimit(0.1*MeV),
00049 fRadiusConst(1.08*fermi),
00050 fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
00051 fDiffractionXsc(0.0), fHadronNucleonXsc(0.0)
00052 {
00053 theProton = G4Proton::Proton();
00054 theNeutron = G4Neutron::Neutron();
00055 hnXsc = new G4HadronNucleonXsc();
00056 }
00057
00058
00059 G4GGNuclNuclCrossSection::~G4GGNuclNuclCrossSection()
00060 {
00061 delete hnXsc;
00062 }
00063
00064 void
00065 G4GGNuclNuclCrossSection::CrossSectionDescription(std::ostream& outFile) const
00066 {
00067 outFile << "G4GGNuclNuclCrossSection calculates total, inelastic and\n"
00068 << "elastic cross sections for nucleus-nucleus collisions using\n"
00069 << "the Glauber model with Gribov corrections. It is valid for\n"
00070 << "all incident energies above 100 keV./n";
00071 }
00072
00073 G4bool
00074 G4GGNuclNuclCrossSection::IsElementApplicable(const G4DynamicParticle*,
00075 G4int, const G4Material*)
00076 {
00077 G4bool applicable = true;
00078
00079
00080
00081 return applicable;
00082 }
00083
00085
00086
00087
00088
00089
00090
00091
00092 G4double G4GGNuclNuclCrossSection::
00093 GetElementCrossSection(const G4DynamicParticle* aParticle, G4int Z,
00094 const G4Material*)
00095 {
00096 G4int A = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
00097 return GetZandACrossSection(aParticle, Z, A);
00098 }
00099
00101
00102
00103
00104
00105
00106
00107
00108
00109 G4double G4GGNuclNuclCrossSection::
00110 GetZandACrossSection(const G4DynamicParticle* aParticle,
00111 G4int tZ, G4int tA)
00112 {
00113 G4double xsection;
00114 G4double sigma;
00115 G4double cofInelastic = 2.4;
00116 G4double cofTotal = 2.0;
00117 G4double nucleusSquare;
00118 G4double cB;
00119 G4double ratio;
00120
00121 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
00122 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
00123
00124 G4double pTkin = aParticle->GetKineticEnergy();
00125 pTkin /= pA;
00126
00127 G4double pN = pA - pZ;
00128 if( pN < 0. ) pN = 0.;
00129
00130 G4double tN = tA - tZ;
00131 if( tN < 0. ) tN = 0.;
00132
00133 G4double tR = GetNucleusRadius( G4double(tZ),G4double(tA) );
00134 G4double pR = GetNucleusRadius(pZ,pA);
00135
00136 cB = GetCoulombBarier(aParticle, G4double(tZ), G4double(tA), pR, tR);
00137
00138 if ( cB > 0. )
00139 {
00140 G4DynamicParticle* dProton = new G4DynamicParticle(theProton,
00141 G4ParticleMomentum(1.,0.,0.),
00142 pTkin);
00143
00144 G4DynamicParticle* dNeutron = new G4DynamicParticle(theNeutron,
00145 G4ParticleMomentum(1.,0.,0.),
00146 pTkin);
00147
00148 sigma = (pZ*tZ+pN*tN)*hnXsc->GetHadronNucleonXscNS(dProton, theProton);
00149
00150 G4double ppInXsc = hnXsc->GetInelasticHadronNucleonXsc();
00151
00152 sigma += (pZ*tN+pN*tZ)*hnXsc->GetHadronNucleonXscNS(dNeutron, theProton);
00153
00154 G4double npInXsc = hnXsc->GetInelasticHadronNucleonXsc();
00155
00156 delete dProton;
00157 delete dNeutron;
00158
00159
00160
00161
00162
00163 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );
00164
00165 ratio = sigma/nucleusSquare;
00166 xsection = nucleusSquare*std::log( 1. + ratio );
00167 fTotalXsc = xsection;
00168 fTotalXsc *= cB;
00169
00170 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00171
00172 fInelasticXsc *= cB;
00173 fElasticXsc = fTotalXsc - fInelasticXsc;
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 sigma = (pZ*tZ+pN*tN)*ppInXsc + (pZ*tN+pN*tZ)*npInXsc;
00187
00188 ratio = sigma/nucleusSquare;
00189 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00190
00191 if (fElasticXsc < 0.) fElasticXsc = 0.;
00192 }
00193 else
00194 {
00195 fInelasticXsc = 0.;
00196 fTotalXsc = 0.;
00197 fElasticXsc = 0.;
00198 fProductionXsc = 0.;
00199 }
00200
00201 return fInelasticXsc;
00202 }
00203
00205
00206
00207
00208 G4double G4GGNuclNuclCrossSection::
00209 GetCoulombBarier(const G4DynamicParticle* aParticle, G4double tZ, G4double tA,
00210 G4double pR, G4double tR)
00211 {
00212 G4double ratio;
00213 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
00214
00215 G4double pTkin = aParticle->GetKineticEnergy();
00216
00217 G4double pM = aParticle->GetDefinition()->GetPDGMass();
00218
00219 G4double tM = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( G4int(tZ), G4int(tA) );
00220 G4double pElab = pTkin + pM;
00221 G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
00222
00223
00224 G4double totTcm = totEcm - pM -tM;
00225
00226 G4double bC = fine_structure_const*hbarc*pZ*tZ;
00227 bC /= pR + tR;
00228 bC /= 2.;
00229
00230
00231
00232
00233 if( totTcm <= bC ) ratio = 0.;
00234 else ratio = 1. - bC/totTcm;
00235
00236
00237 if( ratio < 0.) ratio = 0.;
00238
00239
00240 return ratio;
00241 }
00242
00243
00245
00246
00247
00248 G4double G4GGNuclNuclCrossSection::
00249 GetRatioSD(const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
00250 {
00251 G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
00252
00253 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
00254 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
00255
00256 G4double pTkin = aParticle->GetKineticEnergy();
00257 pTkin /= pA;
00258
00259 G4double pN = pA - pZ;
00260 if( pN < 0. ) pN = 0.;
00261
00262 G4double tN = tA - tZ;
00263 if( tN < 0. ) tN = 0.;
00264
00265 G4double tR = GetNucleusRadius(tZ,tA);
00266 G4double pR = GetNucleusRadius(pZ,pA);
00267
00268 sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
00269 (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
00270
00271 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );
00272 ratio = sigma/nucleusSquare;
00273 fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
00274 G4double difratio = ratio/(1.+ratio);
00275
00276 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
00277
00278 if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
00279 else ratio = 0.;
00280
00281 return ratio;
00282 }
00283
00285
00286
00287
00288 G4double G4GGNuclNuclCrossSection::
00289 GetRatioQE(const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
00290 {
00291 G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
00292
00293 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
00294 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
00295
00296 G4double pTkin = aParticle->GetKineticEnergy();
00297 pTkin /= pA;
00298
00299 G4double pN = pA - pZ;
00300 if( pN < 0. ) pN = 0.;
00301
00302 G4double tN = tA - tZ;
00303 if( tN < 0. ) tN = 0.;
00304
00305 G4double tR = GetNucleusRadius(tZ,tA);
00306 G4double pR = GetNucleusRadius(pZ,pA);
00307
00308 sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
00309 (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
00310
00311 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );
00312 ratio = sigma/nucleusSquare;
00313 fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
00314
00315
00316 ratio = sigma/nucleusSquare;
00317 fProductionXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
00318
00319 if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
00320 else ratio = 0.;
00321 if ( ratio < 0. ) ratio = 0.;
00322
00323 return ratio;
00324 }
00325
00327
00328
00329
00330
00331
00332
00333 G4double
00334 G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
00335 const G4Element* anElement)
00336 {
00337 G4int At = G4lrint(anElement->GetN());
00338 G4int Zt = G4lrint(anElement->GetZ());
00339 return GetHadronNucleonXsc(aParticle, At, Zt);
00340 }
00341
00342
00343
00344
00346
00347
00348
00349
00350
00351
00352 G4double
00353 G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
00354 G4int At, G4int Zt)
00355 {
00356 G4double xsection = 0.;
00357
00358 G4double targ_mass = G4ParticleTable::GetParticleTable()->
00359 GetIonTable()->GetIonMass(Zt, At);
00360 targ_mass = 0.939*GeV;
00361
00362 G4double proj_mass = aParticle->GetMass();
00363 G4double proj_momentum = aParticle->GetMomentum().mag();
00364 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00365
00366 sMand /= GeV*GeV;
00367 proj_momentum /= GeV;
00368 const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
00369
00370 if(pParticle == theNeutron)
00371 {
00372 xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00373 }
00374 else if(pParticle == theProton)
00375 {
00376 xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00377 }
00378
00379 xsection *= millibarn;
00380 return xsection;
00381 }
00382
00383
00384
00386
00387
00388
00389
00390
00391 G4double
00392 G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(G4ParticleDefinition* pParticle,
00393 G4double sMand,
00394 G4ParticleDefinition* tParticle)
00395 {
00396 G4double xsection = 0.;
00397
00398 G4bool proton = (tParticle == theProton);
00399 G4bool neutron = (tParticle == theNeutron);
00400
00401
00402
00403 G4double s0 = 5.38*5.38;
00404 G4double eta1 = 0.458;
00405 G4double eta2 = 0.458;
00406 G4double B = 0.308;
00407
00408
00409
00410 if(pParticle == theNeutron)
00411 {
00412 if ( proton )
00413 {
00414 xsection = ( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00415 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00416 }
00417 if ( neutron )
00418 {
00419 xsection = (35.45 + B*std::pow(std::log(sMand/s0),2.)
00420 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00421 }
00422 }
00423 else if(pParticle == theProton)
00424 {
00425 if ( proton )
00426 {
00427 xsection = (35.45 + B*std::pow(std::log(sMand/s0),2.)
00428 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00429
00430 }
00431 if ( neutron )
00432 {
00433 xsection = (35.80 + B*std::pow(std::log(sMand/s0),2.)
00434 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00435 }
00436 }
00437 xsection *= millibarn;
00438 return xsection;
00439 }
00440
00441
00443
00444
00445
00446
00447
00448 G4double
00449 G4GGNuclNuclCrossSection::GetHadronNucleonXscNS(G4ParticleDefinition* pParticle,
00450 G4double pTkin,
00451 G4ParticleDefinition* tParticle)
00452 {
00453 G4double xsection(0);
00454
00455 G4double A0, B0;
00456 G4double hpXscv(0);
00457 G4double hnXscv(0);
00458
00459 G4double targ_mass = tParticle->GetPDGMass();
00460 G4double proj_mass = pParticle->GetPDGMass();
00461
00462 G4double proj_energy = proj_mass + pTkin;
00463 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
00464
00465 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00466
00467 sMand /= GeV*GeV;
00468 proj_momentum /= GeV;
00469 proj_energy /= GeV;
00470 proj_mass /= GeV;
00471
00472
00473
00474
00475
00476
00477
00478
00479 if( proj_momentum >= 373.)
00480 {
00481 return GetHadronNucleonXscPDG(pParticle,sMand,tParticle);
00482 }
00483 else if( proj_momentum >= 10. )
00484
00485 {
00486
00487
00488
00489 if (proj_momentum >= 10.) {
00490 B0 = 7.5;
00491 A0 = 100. - B0*std::log(3.0e7);
00492
00493 xsection = A0 + B0*std::log(proj_energy) - 11
00494 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
00495 0.93827*0.93827,-0.165);
00496 }
00497 }
00498 else
00499 {
00500 if(pParticle == tParticle)
00501 {
00502 if( proj_momentum < 0.73 )
00503 {
00504 hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
00505 }
00506 else if( proj_momentum < 1.05 )
00507 {
00508 hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
00509 (std::log(proj_momentum/0.73));
00510 }
00511 else
00512 {
00513 hnXscv = 39.0 +
00514 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
00515 }
00516 xsection = hnXscv;
00517 }
00518 else
00519 {
00520 if( proj_momentum < 0.8 )
00521 {
00522 hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
00523 }
00524 else if( proj_momentum < 1.4 )
00525 {
00526 hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
00527 }
00528 else
00529 {
00530 hpXscv = 33.3+
00531 20.8*(std::pow(proj_momentum,2.0)-1.35)/
00532 (std::pow(proj_momentum,2.50)+0.95);
00533 }
00534 xsection = hpXscv;
00535 }
00536 }
00537 xsection *= millibarn;
00538 return xsection;
00539 }
00540
00542
00543
00544
00545 G4double
00546 G4GGNuclNuclCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
00547 G4int At, G4int Zt)
00548 {
00549 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
00550 G4int absPDGcode = std::abs(PDGcode);
00551 G4double Elab = aParticle->GetTotalEnergy();
00552
00553 G4double Plab = aParticle->GetMomentum().mag();
00554
00555
00556 Elab /= GeV;
00557 Plab /= GeV;
00558
00559 G4double LogPlab = std::log( Plab );
00560 G4double sqrLogPlab = LogPlab * LogPlab;
00561
00562
00563
00564 G4double NumberOfTargetProtons = Zt;
00565 G4double NumberOfTargetNucleons = At;
00566 G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
00567
00568 if(NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.;
00569
00570 G4double Xtotal = 0., Xelastic = 0., Xinelastic =0.;
00571
00572 if( absPDGcode > 1000 )
00573 {
00574 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
00575 0.522*sqrLogPlab - 4.51*LogPlab;
00576
00577 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
00578 0.513*sqrLogPlab - 4.27*LogPlab;
00579
00580 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
00581 0.169*sqrLogPlab - 1.85*LogPlab;
00582
00583 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
00584 0.169*sqrLogPlab - 1.85*LogPlab;
00585
00586 Xtotal = ( NumberOfTargetProtons * XtotPP +
00587 NumberOfTargetNeutrons * XtotPN );
00588
00589 Xelastic = ( NumberOfTargetProtons * XelPP +
00590 NumberOfTargetNeutrons * XelPN );
00591 }
00592
00593 Xinelastic = Xtotal - Xelastic;
00594 if(Xinelastic < 0.) Xinelastic = 0.;
00595
00596 return Xinelastic*= millibarn;
00597 }
00598
00600
00601
00602
00603 G4double
00604 G4GGNuclNuclCrossSection::GetNucleusRadius(const G4DynamicParticle* ,
00605 const G4Element* anElement)
00606 {
00607 G4double At = anElement->GetN();
00608 G4double oneThird = 1.0/3.0;
00609 G4double cubicrAt = std::pow (At, oneThird);
00610
00611 G4double R;
00612 R = fRadiusConst*cubicrAt;
00613
00614 G4double meanA = 21.;
00615 G4double tauA1 = 40.;
00616 G4double tauA2 = 10.;
00617 G4double tauA3 = 5.;
00618
00619 G4double a1 = 0.85;
00620 G4double b1 = 1. - a1;
00621
00622 G4double b2 = 0.3;
00623 G4double b3 = 4.;
00624
00625 if (At > 20.)
00626 {
00627 R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
00628 }
00629 else if (At > 3.5)
00630 {
00631 R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
00632 }
00633 else
00634 {
00635 R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
00636 }
00637
00638 return R;
00639 }
00640
00642
00643
00644
00645 G4double
00646 G4GGNuclNuclCrossSection::GetNucleusRadius(G4double Zt, G4double At)
00647 {
00648 G4double R;
00649 R = GetNucleusRadiusDE(Zt,At);
00650
00651
00652 return R;
00653 }
00654
00656
00657 G4double
00658 G4GGNuclNuclCrossSection::GetNucleusRadiusGG(G4double At)
00659 {
00660 G4double oneThird = 1.0/3.0;
00661 G4double cubicrAt = std::pow (At, oneThird);
00662
00663 G4double R;
00664 R = fRadiusConst*cubicrAt;
00665
00666 G4double meanA = 20.;
00667 G4double tauA = 20.;
00668
00669 if ( At > 20.)
00670 {
00671 R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
00672 }
00673 else
00674 {
00675 R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) );
00676 }
00677
00678 return R;
00679 }
00680
00682
00683
00684
00685 G4double
00686 G4GGNuclNuclCrossSection::GetNucleusRadiusDE(G4double Z, G4double A)
00687 {
00688
00689
00690 G4double R, r0, a11, a12, a13, a2, a3;
00691
00692 a11 = 1.26;
00693 a12 = 1.;
00694 a13 = 1.12;
00695 a2 = 1.1;
00696 a3 = 1.;
00697
00698
00699
00700 if (A < 50.)
00701 {
00702 if (std::abs(A-1.) < 0.5) return 0.89*fermi;
00703 else if(std::abs(A-2.) < 0.5) return 2.13*fermi;
00704 else if(std::abs(Z-1.) < 0.5 && std::abs(A-3.) < 0.5) return 1.80*fermi;
00705
00706 else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96*fermi;
00707 else if(std::abs(Z-2.) < 0.5 && std::abs(A-4.) < 0.5) return 1.68*fermi;
00708
00709 else if(std::abs(Z-3.) < 0.5) return 2.40*fermi;
00710 else if(std::abs(Z-4.) < 0.5) return 2.51*fermi;
00711
00712 else if( 10. < A && A <= 16. ) r0 = a11*( 1 - std::pow(A, -2./3.) )*fermi;
00713 else if( 15. < A && A <= 20. ) r0 = a12*( 1 - std::pow(A, -2./3.) )*fermi;
00714 else if( 20. < A && A <= 30. ) r0 = a13*( 1 - std::pow(A, -2./3.) )*fermi;
00715 else r0 = a2*fermi;
00716
00717 R = r0*std::pow( A, 1./3. );
00718 }
00719 else
00720 {
00721 r0 = a3*fermi;
00722
00723 R = r0*std::pow(A, 0.27);
00724 }
00725 return R;
00726 }
00727
00728
00730
00731
00732
00733 G4double
00734 G4GGNuclNuclCrossSection::GetNucleusRadiusRMS(G4double Z, G4double A)
00735 {
00736
00737 if (std::abs(A-1.) < 0.5) return 0.89*fermi;
00738 else if(std::abs(A-2.) < 0.5) return 2.13*fermi;
00739 else if(std::abs(Z-1.) < 0.5 && std::abs(A-3.) < 0.5) return 1.80*fermi;
00740
00741 else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96*fermi;
00742 else if(std::abs(Z-2.) < 0.5 && std::abs(A-4.) < 0.5) return 1.68*fermi;
00743
00744 else if(std::abs(Z-3.) < 0.5) return 2.40*fermi;
00745 else if(std::abs(Z-4.) < 0.5) return 2.51*fermi;
00746
00747 else return 1.24*std::pow(A, 0.28 )*fermi;
00748 }
00749
00750
00752
00753
00754
00755 G4double G4GGNuclNuclCrossSection::CalculateEcmValue(const G4double mp,
00756 const G4double mt,
00757 const G4double Plab)
00758 {
00759 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
00760 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
00761
00762
00763
00764 return Ecm ;
00765 }
00766
00767
00769
00770
00771
00772 G4double G4GGNuclNuclCrossSection::CalcMandelstamS(const G4double mp,
00773 const G4double mt,
00774 const G4double Plab)
00775 {
00776 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
00777 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
00778
00779 return sMand;
00780 }
00781
00782
00783