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