50 fRadiusConst(1.08*
fermi),
51 fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
69 outFile <<
"G4GGNuclNuclCrossSection calculates total, inelastic and\n"
70 <<
"elastic cross sections for nucleus-nucleus collisions using\n"
71 <<
"the Glauber model with Gribov corrections. It is valid for\n"
72 <<
"all incident energies above 100 keV./n";
130 if( pN < 0. ) pN = 0.;
133 if( tN < 0. ) tN = 0.;
165 nucleusSquare = cofTotal*
pi*( pR*pR + tR*tR );
167 ratio = sigma/nucleusSquare;
168 xsection = nucleusSquare*std::log( 1. + ratio );
169 fTotalXsc = xsection;
172 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
175 fElasticXsc = fTotalXsc - fInelasticXsc;
188 sigma = (pZ*tZ+pN*tN)*ppInXsc + (pZ*tN+pN*tZ)*npInXsc;
190 ratio = sigma/nucleusSquare;
191 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
193 if (fElasticXsc < 0.) fElasticXsc = 0.;
203 return fInelasticXsc;
223 G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
235 if( totTcm <= bC ) ratio = 0.;
236 else ratio = 1. - bC/totTcm;
239 if( ratio < 0.) ratio = 0.;
253 G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
262 if( pN < 0. ) pN = 0.;
265 if( tN < 0. ) tN = 0.;
273 nucleusSquare = cofTotal*
pi*( pR*pR + tR*tR );
274 ratio = sigma/nucleusSquare;
275 fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
276 G4double difratio = ratio/(1.+ratio);
278 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
280 if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
293 G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
302 if( pN < 0. ) pN = 0.;
305 if( tN < 0. ) tN = 0.;
313 nucleusSquare = cofTotal*
pi*( pR*pR + tR*tR );
314 ratio = sigma/nucleusSquare;
315 fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
318 ratio = sigma/nucleusSquare;
319 fProductionXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
321 if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
323 if ( ratio < 0. ) ratio = 0.;
361 GetIonTable()->GetIonMass(Zt, At);
362 targ_mass = 0.939*
GeV;
369 proj_momentum /=
GeV;
372 if(pParticle == theNeutron)
374 xsection =
G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
376 else if(pParticle == theProton)
378 xsection =
G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
412 if(pParticle == theNeutron)
416 xsection = ( 35.80 + B*std::pow(std::log(sMand/s0),2.)
417 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
421 xsection = (35.45 + B*std::pow(std::log(sMand/s0),2.)
422 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
425 else if(pParticle == theProton)
429 xsection = (35.45 + B*std::pow(std::log(sMand/s0),2.)
430 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
435 xsection = (35.80 + B*std::pow(std::log(sMand/s0),2.)
436 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
464 G4double proj_energy = proj_mass + pTkin;
465 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
470 proj_momentum /=
GeV;
481 if( proj_momentum >= 373.)
485 else if( proj_momentum >= 10. )
491 if (proj_momentum >= 10.) {
493 A0 = 100. - B0*std::log(3.0e7);
495 xsection = A0 + B0*std::log(proj_energy) - 11
496 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
497 0.93827*0.93827,-0.165);
502 if(pParticle == tParticle)
504 if( proj_momentum < 0.73 )
506 hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
508 else if( proj_momentum < 1.05 )
510 hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
511 (std::log(proj_momentum/0.73));
516 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
522 if( proj_momentum < 0.8 )
524 hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
526 else if( proj_momentum < 1.4 )
528 hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
533 20.8*(std::pow(proj_momentum,2.0)-1.35)/
534 (std::pow(proj_momentum,2.50)+0.95);
552 G4int absPDGcode = std::abs(PDGcode);
561 G4double LogPlab = std::log( Plab );
562 G4double sqrLogPlab = LogPlab * LogPlab;
566 G4double NumberOfTargetProtons = Zt;
567 G4double NumberOfTargetNucleons = At;
568 G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
570 if(NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.;
572 G4double Xtotal = 0., Xelastic = 0., Xinelastic =0.;
574 if( absPDGcode > 1000 )
576 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
577 0.522*sqrLogPlab - 4.51*LogPlab;
579 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
580 0.513*sqrLogPlab - 4.27*LogPlab;
582 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
583 0.169*sqrLogPlab - 1.85*LogPlab;
585 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
586 0.169*sqrLogPlab - 1.85*LogPlab;
588 Xtotal = ( NumberOfTargetProtons * XtotPP +
589 NumberOfTargetNeutrons * XtotPN );
591 Xelastic = ( NumberOfTargetProtons * XelPP +
592 NumberOfTargetNeutrons * XelPN );
595 Xinelastic = Xtotal - Xelastic;
596 if(Xinelastic < 0.) Xinelastic = 0.;
611 G4double cubicrAt = std::pow (At, oneThird);
614 R = fRadiusConst*cubicrAt;
629 R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
633 R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
637 R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
663 G4double cubicrAt = std::pow (At, oneThird);
666 R = fRadiusConst*cubicrAt;
673 R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
677 R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) );
692 G4double R, r0, a11, a12, a13, a2, a3;
704 if (std::abs(A-1.) < 0.5)
return 0.89*
fermi;
705 else if(std::abs(A-2.) < 0.5)
return 2.13*
fermi;
706 else if(std::abs(Z-1.) < 0.5 && std::abs(A-3.) < 0.5)
return 1.80*
fermi;
708 else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5)
return 1.96*
fermi;
709 else if(std::abs(Z-2.) < 0.5 && std::abs(A-4.) < 0.5)
return 1.68*
fermi;
711 else if(std::abs(Z-3.) < 0.5)
return 2.40*
fermi;
712 else if(std::abs(Z-4.) < 0.5)
return 2.51*
fermi;
714 else if( 10. < A && A <= 16. ) r0 = a11*( 1 - std::pow(A, -2./3.) )*
fermi;
715 else if( 15. < A && A <= 20. ) r0 = a12*( 1 - std::pow(A, -2./3.) )*
fermi;
716 else if( 20. < A && A <= 30. ) r0 = a13*( 1 - std::pow(A, -2./3.) )*
fermi;
719 R = r0*std::pow( A, 1./3. );
725 R = r0*std::pow(A, 0.27);
739 if (std::abs(A-1.) < 0.5)
return 0.89*
fermi;
740 else if(std::abs(A-2.) < 0.5)
return 2.13*
fermi;
741 else if(std::abs(Z-1.) < 0.5 && std::abs(A-3.) < 0.5)
return 1.80*
fermi;
743 else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5)
return 1.96*
fermi;
744 else if(std::abs(Z-2.) < 0.5 && std::abs(A-4.) < 0.5)
return 1.68*
fermi;
746 else if(std::abs(Z-3.) < 0.5)
return 2.40*
fermi;
747 else if(std::abs(Z-4.) < 0.5)
return 2.51*
fermi;
749 else return 1.24*std::pow(A, 0.28 )*
fermi;
761 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
762 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
778 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
779 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *)
G4double CalcMandelstamS(const G4double, const G4double, const G4double)
G4double GetHadronNucleonXscPDG(G4ParticleDefinition *, G4double sMand, G4ParticleDefinition *)
G4double GetNucleusRadiusRMS(G4double Z, G4double A)
G4double GetRatioQE(const G4DynamicParticle *, G4double At, G4double Zt)
G4int GetPDGEncoding() const
G4ParticleDefinition * GetDefinition() const
static G4NistManager * Instance()
G4double GetHadronNucleonXscNS(const G4DynamicParticle *, const G4ParticleDefinition *)
G4double GetNucleusRadiusGG(G4double At)
G4IonTable * GetIonTable() const
G4double GetHNinelasticXscVU(const G4DynamicParticle *, G4int At, G4int Zt)
G4double GetRatioSD(const G4DynamicParticle *, G4double At, G4double Zt)
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
static G4Proton * Proton()
static G4Neutron * Neutron()
G4GGNuclNuclCrossSection()
virtual ~G4GGNuclNuclCrossSection()
G4double CalculateEcmValue(const G4double, const G4double, const G4double)
#define G4_DECLARE_XS_FACTORY(cross_section)
G4double GetZandACrossSection(const G4DynamicParticle *, G4int Z, G4int A)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4double GetCoulombBarier(const G4DynamicParticle *, G4double Z, G4double A, G4double pR, G4double tR)
G4double GetNucleusRadiusDE(G4double Z, G4double A)
virtual void CrossSectionDescription(std::ostream &) const
virtual G4bool IsElementApplicable(const G4DynamicParticle *, G4int Z, const G4Material *)
G4double GetHadronNucleonXsc(const G4DynamicParticle *, const G4Element *)
G4double GetNucleusRadius(const G4DynamicParticle *, const G4Element *)
G4double GetPDGCharge() const
G4ThreeVector G4ParticleMomentum
G4double GetHadronNucleonXscNS(G4ParticleDefinition *, G4double pTkin, G4ParticleDefinition *)
G4double GetInelasticHadronNucleonXsc()
G4ThreeVector GetMomentum() const
G4int GetBaryonNumber() const