68 mass2GeV2= massGeV*massGeV;
69 AtomicWeight =
G4int(AWeight + 0.5);
71 DefineNucleusParameters(AWeight);
73 limitQ2 = 35./(R1*R1);
81 TableQ2[ii] = TableQ2[ii-1]+dQ2;
84 massA = AWeight*amu_c2/
GeV;
91 G4double plab2= eGeV[kk]*(eGeV[kk] + 2.0*massGeV);
92 G4double Q2m = 4.0*plab2*massA2/(mass2GeV2 + massA2 + 2.*massA2*elab);
97 TableCrossSec[ONQ2*kk] = 0.0;
109 switch (AtomicWeight)
203 R1 = 4.45*std::pow(A - 1.,0.309)*0.9;
204 R2 = 2.3 *std::pow(A, 0.36);
206 if(A < 100 && A > 3) Pnucl = 0.176 + 0.00275*
A;
212 if(A >= 100) Aeff = 0.7;
213 else if(A < 100 && A > 75) Aeff = 1.5 - 0.008*
A;
273 G4int ic[
NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311,
274 3122,3222,3112,3212,3312,3322,3334,
275 -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334};
277 G4int id[
NHADRONS] = {2,3,6,0,4,5,4,4,4,5,
281 G4int id1[
NHADRONS] = {3,4,1,0,5,6,5,5,5,6,
298 char* dirName = getenv(
"G4PhysListDocDir");
303 outFile.open(pathName);
304 outFile <<
"<html>\n";
305 outFile <<
"<head>\n";
307 outFile <<
"<title>Description of G4ElasticHadrNucleusHE Model</title>\n";
308 outFile <<
"</head>\n";
309 outFile <<
"<body>\n";
311 outFile <<
"G4ElasticHadrNucleusHE is a hadron-nucleus elastic scattering\n"
312 <<
"model developed by N. Starkov which uses a Glauber model\n"
313 <<
"parameterization to calculate the final state. It is valid\n"
314 <<
"for all hadrons with incident energies above 1 GeV.\n";
316 outFile <<
"</body>\n";
317 outFile <<
"</html>\n";
331 for(
G4int k = 0; k < 93; k++)
356 G4cout<<
" G4ElasticHadrNucleusHE::SampleT: "
358 <<
" at Z= " << Z <<
" A= " << N
359 <<
" plab(GeV)= " << plab
364 for( idx = 0 ; idx <
NHADRONS; idx++)
371 if( idx >= NHADRONS )
return Q2;
401 G4cout<<
" G4ElasticHadrNucleusHE::SampleT: new record " << idx
420 G4cout<<
" SampleT: Q2(GeV^2)= "<<Q2<<
" t/tmax= " << Q2/Q2max <<
G4endl;
456 G4cout<<
"Q2_2: ekin plab "<<ekin<<
" "<<plab<<
" tmax "<<tmax<<
G4endl;
460 for( NumbOnE = 0; NumbOnE <
NENERGY-1; NumbOnE++ )
487 G4cout<<
"1 plab T "<<plab<<
" "<<T<<
" sigTot B ReIm "
496 G4cout<<
" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
497 <<
" length= " << length
498 <<
" Q2max= " << Q2max
499 <<
" ekin= " << ekin <<
G4endl;
523 for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ )
525 if(Rand <= pElD->TableCrossSec[index+iNumbQ2])
break;
527 Q2 =
GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
529 if(tmax < Q2max) Q2 *= tmax/Q2max;
532 G4cout<<
" HadrNucleusQ2_2(2): Q2= "<<Q2<<
" iNumbQ2= " << iNumbQ2
533 <<
" rand= " << Rand <<
G4endl;
555 G4cout <<
"GetQ2_2 kk= " << kk <<
" X2= " << X2 <<
" X3= " << X3
556 <<
" F2= " << F2 <<
" F3= " << F3 <<
" Rndm= " << ranUni <<
G4endl;
558 if(kk == 1 || kk == 0)
572 G4double D0 = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3;
575 G4cout <<
" X1= " << X1 <<
" F1= " << F1 <<
" D0= "
578 if(std::abs(D0) < 0.00000001)
580 ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
584 G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1;
585 G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*
F22;
586 G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22
587 -X1*F2*F32-X2*F3*F12-X3*F1*
F22;
588 ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0;
607 for(ii = 1; ii<
ONQ2; ii++)
612 for(jj = 0; jj<20; jj++)
617 curSum += curSec*aSimp;
619 if(aSimp > 3) aSimp = 2;
624 <<
" Diff "<<curSec<<
" totSum "<<totSum<<
G4endl;
634 G4cout<<
" GetHeavy: Q2 dQ2 totSum "<<Q2l<<
" "<<
dQ2<<
" "<<totSum
636 <<curSec<<
" totSum "<< totSum<<
" DTot "
677 G4cout <<
"GetFq2: Stot= " << Stot <<
" Bhad= " << Bhad
679 <<
" Asq= " << Asq <<
G4endl;
691 G4double Unucl = Stot/twopi/Norm*R13;
694 G4double FiH = std::asin(HadrReIm/Rho2);
698 G4cout <<
"UnucRho2= " << UnucRho2 <<
" FiH= " << FiH <<
" NN2= " << NN2
699 <<
" Norm= " << Norm <<
G4endl;
711 G4int i1, i2, j1, j2;
713 for(i1 = 1; i1<= Nucleus; i1++)
715 N1 = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2;
720 for(i2 = 1; i2<=Nucleus; i2++)
722 N2 = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2;
725 for(j2=0; j2<= i2; j2++)
728 exp2 = 1/(j2/R22B+(i2-j2)/R12B);
731 for(j1=0; j1<=i1; j1++)
733 exp1 = 1/(j1/R22B+(i1-j1)/R12B);
736 Prod3 = Prod3+N4*exp1*exp2*
737 (1-std::exp(-Q2*dddd/4))/dddd*4*
SetBinom[i1][j1];
739 Prod2 = Prod2 +Prod3*N5*
SetBinom[i2][j2];
741 Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2));
743 if (std::fabs(Prod2*N2/Prod1)<prec)
break;
745 Prod0 = Prod0 + Prod1*N1;
746 if(std::fabs(N1*Prod1/Prod0) < prec)
break;
753 G4cout <<
"GetLightFq2 Z= " << Z <<
" A= " << Nucleus
754 <<
" Q2= " << Q2 <<
" Res= " << Prod0 <<
G4endl;
809 G4double R23Ap = R22*R2/R22Ap*Pnuclp;
813 G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);
817 std::sqrt((R12+R22)/2)/4));
824 G4double Unucl = Stot/twopi/Norm*R13;
825 G4double UnuclScr = Stot/twopi/Normp*R13Ap;
836 G4double Prod1, Tot1=0, medTot, DTot1, DmedTot;
839 for( i=1; i<=NWeight; i++)
841 N = -N*Unucl*(NWeight-i+1)/i*Rho2;
843 Prod1 = std::exp(-theQ2/i*R12B/4)/i*R12B;
846 for(
G4int l=1; l<=i; l++)
848 exp1 = l/R22B+(i-l)/R12B;
849 N4 = -N4*(i-l+1)/l*N2;
850 Prod1 = Prod1+N4/exp1*std::exp(-theQ2/exp1/4);
851 medTot = medTot+N4/exp1;
854 ReElasticAmpl0 = ReElasticAmpl0+Prod1*N*std::sin(FiH*i);
855 ImElasticAmpl0 = ImElasticAmpl0+Prod1*N*std::cos(FiH*i);
856 Tot1 = Tot1+medTot*N*std::cos(FiH*i);
857 if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001)
break;
860 ImElasticAmpl0 = ImElasticAmpl0*
pi/2.568;
861 ReElasticAmpl0 = ReElasticAmpl0*
pi/2.568;
862 Tot1 = Tot1*twopi/2.568;
865 G4double C2 = 2*R23Ap*R13Ap/2*DDSec2p;
872 Din1 = 0.5*(C1*std::exp(-theQ2/8*R12Ap)/2*R12Ap-
873 C2/R12ApdR22Ap*std::exp(-theQ2/4/R12ApdR22Ap)+
874 C3*R22Ap/2*std::exp(-theQ2/8*R22Ap));
876 DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2);
886 for( i = 1; i<= NWeight-2; i++)
888 N1p = -N1p*UnuclScr*(NWeight-i-1)/i*Rho2;
892 for(
G4int l = 0; l<=i; l++)
894 if(l == 0) BinCoeff = 1;
895 else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l;
897 exp1 = l/R22B+(i-l)/R12B;
899 exp2p = exp1+R12ApdR22Ap;
902 Din2 = Din2 + N2p*BinCoeff*
903 (C1/exp1p*std::exp(-theQ2/4/exp1p)-
904 C2/exp2p*std::exp(-theQ2/4/exp2p)+
905 C3/exp3p*std::exp(-theQ2/4/exp3p));
907 DmedTot = DmedTot + N2p*BinCoeff*
908 (C1/exp1p-C2/exp2p+C3/exp3p);
913 Din1 = Din1+Din2*N1p/(i+2)/(i+1)*std::cos(FiH*i);
914 DTot1 = DTot1+DmedTot*N1p/(i+2)/(i+1)*std::cos(FiH*i);
916 if(std::fabs(Din2*N1p/Din1) < 0.000001)
break;
919 Din1 = -Din1*NWeight*(NWeight-1)
922 DTot1 = DTot1*NWeight*(NWeight-1)
930 G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
931 (ImElasticAmpl0+Din1)*
932 (ImElasticAmpl0+Din1))*2/4/
pi;
937 aAIm = ImElasticAmpl0;
940 return DiffCrSec2*1.0;
959 <<
" E(GeV)= " <<
HadrEnergy <<
" sqrS= " << sqrS
974 TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165);
990 TotN = 33.0 + 25.5*A0*A0;
995 TotN = 33.0 + 30.*A0*A0*A0*A0;
1009 TotP = 23.0 + 40.*A0*A0;
1030 HadrReIm = 0.13*(logS - 5.8579332)*std::pow(sHadr,-0.18);
1077 HadrTot = 5.2+5.2*logE + 123.2/sqrS;
1081 HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8);
1083 HadrReIm = 0.6*(logS - 5.8579332)*std::pow(sHadr,-0.25);
1130 TotP = 3.2*std::exp(-x*x) + 12.*std::exp(-y*y) + 27.5;
1135 TotP = 88*(logE+0.2877)*(logE+0.2877)+14.0;
1141 TotP = 20. + 180.*std::exp(-x*x);
1146 TotN = 10.6 + 2.*logE + 30.*std::pow(
HadrEnergy,-0.43);
1152 TotN = 36.1+0.079 - 4.313*logE + 3.*std::exp(-x*x) +
1159 TotN = 36.1 + 10.*std::exp(-x*x) + 24*std::exp(-y*y);
1164 TotN = 26. + 110.*x*x;
1169 TotN = 28.0 + 40.*std::exp(-x*x);
1174 HadrReIm = 0.2*(logS - 4.6051702)*std::pow(sHadr,-0.15);
1185 else HadrSlope = 1.0+1.76*logS - 2.84/sqrS;
1187 HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1);
1194 HadrTot = 10+1.8*logE + 25./sqrS;
1196 HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1);
1216 static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0};
1217 static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002};
1218 static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0};
1219 static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25};
1220 static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8};
1223 static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0};
1224 static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01};
1225 static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003};
1226 static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5};
1227 static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8};
1230 static const G4double EnP[2]={1.5,4.0};
1231 static const G4double C0P[2]={0.001,0.0005};
1232 static const G4double C1P[2]={0.003,0.001};
1233 static const G4double B0P[2]={2.5,4.5};
1234 static const G4double B1P[2]={1.0,4.0};
1237 static const G4double EnPP[4]={1.0,2.0,3.0,4.0};
1238 static const G4double C0PP[4]={0.0,0.0,0.0,0.0};
1239 static const G4double C1PP[4]={0.15,0.08,0.02,0.01};
1240 static const G4double B0PP[4]={1.5,2.8,3.8,3.8};
1241 static const G4double B1PP[4]={0.8,1.6,3.6,4.6};
1244 static const G4double EnPPN[4]={1.0,2.0,3.0,4.0};
1245 static const G4double C0PPN[4]={0.0,0.0,0.0,0.0};
1246 static const G4double C1PPN[4]={0.0,0.0,0.0,0.0};
1247 static const G4double B0PPN[4]={1.5,2.8,3.8,3.8};
1248 static const G4double B1PPN[4]={0.8,1.6,3.6,4.6};
1251 static const G4double EnK[4]={1.4,2.33,3.0,5.0};
1252 static const G4double C0K[4]={0.0,0.0,0.0,0.0};
1253 static const G4double C1K[4]={0.01,0.007,0.005,0.003};
1254 static const G4double B0K[4]={1.5,2.0,3.8,3.8};
1255 static const G4double B1K[4]={1.6,1.6,1.6,1.6};
1258 static const G4double EnKM[2]={1.4,4.0};
1259 static const G4double C0KM[2]={0.006,0.002};
1260 static const G4double C1KM[2]={0.00,0.00};
1261 static const G4double B0KM[2]={2.5,3.5};
1262 static const G4double B1KM[2]={1.6,1.6};
1333 G4cout<<
"1 GetKin.: HadronName MomentumH "
1360 G4cout<<
"3 GetKin. : NumberH "<<NumberH
1361 <<
" Bound.P[NumberH] "<<
BoundaryP[NumberH]
1384 <<
" Fdistr "<<Fdistr<<
G4endl;
1396 while(std::fabs(delta) > 0.0001)
1401 DDD0 = (DDD0+DDD1)*0.5;
1406 DDD0 = (DDD0+DDD2)*0.5;
1447 for(N = 0; N < 240; N++)
1451 for( M = 0; M <= N; M++ )
1455 if ( ( N > 0 ) && ( N > M ) && M > 0 )
G4int HadronType[NHADRONS]
G4double SampleT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4double GetQ2(G4double Ran)
G4NistManager * nistManager
void DefineHadronValues(G4int Z)
G4double GetDistrFun(G4double Q2)
G4double HadronProtonQ2(const G4ParticleDefinition *aHadron, G4double inLabMom)
G4int GetPDGEncoding() const
G4double LowEdgeEnergy[NENERGY]
const G4String & GetModelName() const
static G4NistManager * Instance()
const G4String & GetParticleName() const
void GetKinematics(const G4ParticleDefinition *aHadron, G4double MomentumH)
G4double SetBinom[240][240]
G4GLOB_DLL std::ostream G4cout
G4ElasticHadrNucleusHE(const G4String &name="hElasticGlauber")
G4double TableCrossSec[NQTABLE]
G4double HadrNucDifferCrSec(G4int Z, G4int Nucleus, G4double Q2)
G4double lowestEnergyLimit
static G4Proton * Proton()
virtual void Description() const
static const G4double A[nN]
static const G4int NENERGY
virtual ~G4ElasticHadrNucleusHE()
G4double GetPDGMass() const
G4ElasticData(const G4ParticleDefinition *h, G4int Z, G4double A, G4double *eGeV)
G4int HadronType1[NHADRONS]
G4double GetLightFq2(G4int Z, G4int A, G4double Q)
G4double GetAtomicMassAmu(const G4String &symb) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetHeavyFq2(G4int Z, G4int Nucleus, G4double *LineFq2)
G4double GetQ2_2(G4int N, G4double *Q, G4double *F, G4double R)
G4double HadronNucleusQ2_2(G4ElasticData *pElD, G4int Z, G4double plabGeV, G4double tmax)
void InterpolateHN(G4int n, const G4double EnP[], const G4double C0P[], const G4double C1P[], const G4double B0P[], const G4double B1P[])
G4double GetFt(G4double Q2)
void DefineNucleusParameters(G4double A)
static const G4int NHADRONS
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4ElasticData * SetOfElasticData[NHADRONS][93]
G4int HadronCode[NHADRONS]