67 G4bool onlyOnceInit =
true;
68 G4bool onlyOnceDestroy =
true;
81 mass2GeV2= massGeV*massGeV;
82 AtomicWeight =
G4lrint(AWeight);
84 DefineNucleusParameters(AWeight);
86 limitQ2 = 35./(R1*R1);
94 TableQ2[ii] = TableQ2[ii-1]+dQ2;
104 G4double plab2= eGeV[kk]*(eGeV[kk] + 2.0*massGeV);
105 G4double Q2m = 4.0*plab2*massA2/(mass2GeV2 + massA2 + 2.*massA2*elab);
110 TableCrossSec[ONQ2*kk] = 0.0;
122 switch (AtomicWeight)
219 if(A < 100 && A > 3) Pnucl = 0.176 + 0.00275*
A;
225 if(A >= 100) Aeff = 0.7;
226 else if(A < 100 && A > 75) Aeff = 1.5 - 0.008*
A;
244 if ( onlyOnceInit ) {
245 for (
int i = 0 ; i<
NHADRONS ; ++i) {
246 for (
int j = 0 ; j<
ZMAX ; ++j ) {
251 onlyOnceInit =
false;
300 G4int ic[
NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311,
301 3122,3222,3112,3212,3312,3322,3334,
302 -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334};
304 G4int id[
NHADRONS] = {2,3,6,0,4,5,4,4,4,5,
308 G4int id1[
NHADRONS] = {3,4,1,0,5,6,5,5,5,6,
326 outFile <<
"G4ElasticHadrNucleusHE is a hadron-nucleus elastic scattering\n"
327 <<
"model developed by N. Starkov which uses a Glauber model\n"
328 <<
"parameterization to calculate the final state. It is valid\n"
329 <<
"for all hadrons with incident energies above 1 GeV.\n";
342 if ( onlyOnceDestroy ) {
354 onlyOnceDestroy =
false;
378 G4cout<<
" G4ElasticHadrNucleusHE::SampleT: "
380 <<
" at Z= " << Z <<
" A= " << N
381 <<
" plab(GeV)= " << plab
386 for( idx = 0 ; idx <
NHADRONS; idx++)
393 if( idx >= NHADRONS )
return Q2;
424 G4cout<<
" G4ElasticHadrNucleusHE::SampleT: new record " << idx
443 G4cout<<
" SampleT: Q2(GeV^2)= "<<Q2<<
" t/tmax= " << Q2/Q2max <<
G4endl;
483 G4cout<<
"Q2_2: ekin plab "<<ekin<<
" "<<plab<<
" tmax "<<tmax<<
G4endl;
487 for( NumbOnE = 0; NumbOnE <
NENERGY-1; NumbOnE++ )
493 if(NumbOnE >= NENERGY-1) { NumbOnE = NENERGY-2; }
515 G4cout<<
"1 plab T "<<plab<<
" "<<T<<
" sigTot B ReIm "
524 G4cout<<
" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
525 <<
" length= " << length
526 <<
" Q2max= " << Q2max
527 <<
" ekin= " << ekin <<
G4endl;
551 for( iNumbQ2 = 1; iNumbQ2<length; ++iNumbQ2)
553 if(Rand <= pElD->TableCrossSec[index+iNumbQ2])
break;
555 if(iNumbQ2 >= ONQ2) { iNumbQ2 = ONQ2 - 1; }
556 Q2 =
GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
558 if(tmax < Q2max) Q2 *= tmax/Q2max;
561 G4cout<<
" HadrNucleusQ2_2(2): Q2= "<<Q2<<
" iNumbQ2= " << iNumbQ2
562 <<
" rand= " << Rand <<
G4endl;
584 G4cout <<
"GetQ2_2 kk= " << kk <<
" X2= " << X2 <<
" X3= " << X3
585 <<
" F2= " << F2 <<
" F3= " << F3 <<
" Rndm= " << ranUni <<
G4endl;
587 if(kk == 1 || kk == 0)
601 G4double D0 = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3;
604 G4cout <<
" X1= " << X1 <<
" F1= " << F1 <<
" D0= "
607 if(std::abs(D0) < 0.00000001)
609 ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
613 G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1;
614 G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*
F22;
615 G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22
616 -X1*F2*F32-X2*F3*F12-X3*F1*
F22;
617 ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0;
636 for(ii = 1; ii<
ONQ2; ii++)
641 for(jj = 0; jj<20; jj++)
646 curSum += curSec*aSimp;
648 if(aSimp > 3) aSimp = 2;
653 <<
" Diff "<<curSec<<
" totSum "<<totSum<<
G4endl;
663 G4cout<<
" GetHeavy: Q2 dQ2 totSum "<<Q2l<<
" "<<
dQ2<<
" "<<totSum
665 <<curSec<<
" totSum "<< totSum<<
" DTot "
706 G4cout <<
"GetFq2: Stot= " << Stot <<
" Bhad= " << Bhad
708 <<
" Asq= " << Asq <<
G4endl;
723 G4double FiH = std::asin(HadrReIm/Rho2);
727 G4cout <<
"UnucRho2= " << UnucRho2 <<
" FiH= " << FiH <<
" NN2= " << NN2
728 <<
" Norm= " << Norm <<
G4endl;
740 G4int i1, i2, j1, j2;
742 for(i1 = 1; i1<= Nucleus; i1++)
744 N1 = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2;
749 for(i2 = 1; i2<=Nucleus; i2++)
751 N2 = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2;
754 for(j2=0; j2<= i2; j2++)
757 exp2 = 1/(j2/R22B+(i2-j2)/R12B);
760 for(j1=0; j1<=i1; j1++)
762 exp1 = 1/(j1/R22B+(i1-j1)/R12B);
765 Prod3 = Prod3+N4*exp1*exp2*
768 Prod2 = Prod2 +Prod3*N5*
SetBinom[i2][j2];
770 Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2));
772 if (std::fabs(Prod2*N2/Prod1)<prec)
break;
774 Prod0 = Prod0 + Prod1*N1;
775 if(std::fabs(N1*Prod1/Prod0) < prec)
break;
782 G4cout <<
"GetLightFq2 Z= " << Z <<
" A= " << Nucleus
783 <<
" Q2= " << Q2 <<
" Res= " << Prod0 <<
G4endl;
838 G4double R23Ap = R22*R2/R22Ap*Pnuclp;
842 G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);
846 std::sqrt((R12+R22)/2)));
865 G4double Prod1, Tot1=0, medTot, DTot1, DmedTot;
868 for( i=1; i<=NWeight; i++)
870 N = -N*Unucl*(NWeight-i+1)/i*Rho2;
872 Prod1 =
G4Exp(-theQ2/i*R12B/4)/i*R12B;
875 for(
G4int l=1; l<=i; l++)
877 exp1 = l/R22B+(i-l)/R12B;
878 N4 = -N4*(i-l+1)/l*N2;
879 Prod1 = Prod1+N4/exp1*
G4Exp(-theQ2/(exp1*4));
880 medTot = medTot+N4/exp1;
883 ReElasticAmpl0 = ReElasticAmpl0+Prod1*N*std::sin(FiH*i);
884 ImElasticAmpl0 = ImElasticAmpl0+Prod1*N*std::cos(FiH*i);
885 Tot1 = Tot1+medTot*N*std::cos(FiH*i);
886 if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001)
break;
889 ImElasticAmpl0 = ImElasticAmpl0*
pi/2.568;
890 ReElasticAmpl0 = ReElasticAmpl0*
pi/2.568;
891 Tot1 = Tot1*
twopi/2.568;
901 Din1 = 0.5*(C1*
G4Exp(-theQ2/8*R12Ap)/2*R12Ap-
902 C2/R12ApdR22Ap*
G4Exp(-theQ2/(4*R12ApdR22Ap))+
903 C3*R22Ap/2*
G4Exp(-theQ2/8*R22Ap));
905 DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2);
915 for( i = 1; i<= NWeight-2; i++)
917 N1p = -N1p*UnuclScr*(NWeight-i-1)/i*Rho2;
921 for(
G4int l = 0; l<=i; l++)
923 if(l == 0) BinCoeff = 1;
924 else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l;
926 exp1 = l/R22B+(i-l)/R12B;
928 exp2p = exp1+R12ApdR22Ap;
931 Din2 = Din2 + N2p*BinCoeff*
932 (C1/exp1p*
G4Exp(-theQ2/(4*exp1p))-
933 C2/exp2p*
G4Exp(-theQ2/(4*exp2p))+
934 C3/exp3p*
G4Exp(-theQ2/(4*exp3p)));
936 DmedTot = DmedTot + N2p*BinCoeff*
937 (C1/exp1p-C2/exp2p+C3/exp3p);
942 Din1 = Din1+Din2*N1p/((i+2)*(i+1))*std::cos(FiH*i);
943 DTot1 = DTot1+DmedTot*N1p/((i+2)*(i+1))*std::cos(FiH*i);
945 if(std::fabs(Din2*N1p/Din1) < 0.000001)
break;
948 Din1 = -Din1*NWeight*(NWeight-1)*4/(Normp*Normp);
950 DTot1 = DTot1*NWeight*(NWeight-1)*4/(Normp*Normp);
957 G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
958 (ImElasticAmpl0+Din1)*
959 (ImElasticAmpl0+Din1))/
twopi;
964 aAIm = ImElasticAmpl0;
986 <<
" E(GeV)= " <<
HadrEnergy <<
" sqrS= " << sqrS
1001 TotP = TotN = 7.5*logE - 40.12525 + 103*
G4Exp(-
G4Log(sHadr)*0.165);
1017 TotN = 33.0 + 25.5*A0*A0;
1022 TotN = 33.0 + 30.*A0*A0*A0*A0;
1036 TotP = 23.0 + 40.*A0*A0;
1104 HadrTot = 5.2+5.2*logE + 123.2/sqrS;
1162 TotP = 88*(logE+0.2877)*(logE+0.2877)+14.0;
1179 TotN = 36.1+0.079 - 4.313*logE + 3.*
G4Exp(-x*x) +
1186 TotN = 36.1 + 10.*
G4Exp(-x*x) + 24*
G4Exp(-y*y);
1191 TotN = 26. + 110.*x*
x;
1196 TotN = 28.0 + 40.*
G4Exp(-x*x);
1212 else HadrSlope = 1.0+1.76*logS - 2.84/sqrS;
1221 HadrTot = 10+1.8*logE + 25./sqrS;
1243 static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0};
1244 static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002};
1245 static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0};
1246 static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25};
1247 static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8};
1250 static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0};
1251 static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01};
1252 static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003};
1253 static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5};
1254 static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8};
1257 static const G4double EnP[2]={1.5,4.0};
1258 static const G4double C0P[2]={0.001,0.0005};
1259 static const G4double C1P[2]={0.003,0.001};
1260 static const G4double B0P[2]={2.5,4.5};
1261 static const G4double B1P[2]={1.0,4.0};
1264 static const G4double EnPP[4]={1.0,2.0,3.0,4.0};
1265 static const G4double C0PP[4]={0.0,0.0,0.0,0.0};
1266 static const G4double C1PP[4]={0.15,0.08,0.02,0.01};
1267 static const G4double B0PP[4]={1.5,2.8,3.8,3.8};
1268 static const G4double B1PP[4]={0.8,1.6,3.6,4.6};
1271 static const G4double EnPPN[4]={1.0,2.0,3.0,4.0};
1272 static const G4double C0PPN[4]={0.0,0.0,0.0,0.0};
1273 static const G4double C1PPN[4]={0.0,0.0,0.0,0.0};
1274 static const G4double B0PPN[4]={1.5,2.8,3.8,3.8};
1275 static const G4double B1PPN[4]={0.8,1.6,3.6,4.6};
1278 static const G4double EnK[4]={1.4,2.33,3.0,5.0};
1279 static const G4double C0K[4]={0.0,0.0,0.0,0.0};
1280 static const G4double C1K[4]={0.01,0.007,0.005,0.003};
1281 static const G4double B0K[4]={1.5,2.0,3.8,3.8};
1282 static const G4double B1K[4]={1.6,1.6,1.6,1.6};
1285 static const G4double EnKM[2]={1.4,4.0};
1286 static const G4double C0KM[2]={0.006,0.002};
1287 static const G4double C1KM[2]={0.00,0.00};
1288 static const G4double B0KM[2]={2.5,3.5};
1289 static const G4double B1KM[2]={1.6,1.6};
1360 G4cout<<
"1 GetKin.: HadronName MomentumH "
1387 G4cout<<
"3 GetKin. : NumberH "<<NumberH
1388 <<
" Bound.P[NumberH] "<<
BoundaryP[NumberH]
1411 <<
" Fdistr "<<Fdistr<<
G4endl;
1423 const G4int maxNumberOfLoops = 10000;
1424 G4int loopCounter = -1;
1425 while ( (std::fabs(delta) > 0.0001) &&
1426 ++loopCounter < maxNumberOfLoops )
1431 DDD0 = (DDD0+DDD1)*0.5;
1436 DDD0 = (DDD0+DDD2)*0.5;
1440 if ( loopCounter >= maxNumberOfLoops ) {
1480 for(N = 0; N < 240; N++)
1484 for( M = 0; M <=
N; M++ )
1488 if ( ( N > 0 ) && ( N > M ) && M > 0 )
static G4ElasticData * SetOfElasticData[NHADRONS][ZMAX]
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
#define G4MUTEXDESTROY(mutex)
virtual ~G4ElasticHadrNucleusHE()
void DefineNucleusParameters(G4double A)
static constexpr double MeV
G4double GetFt(G4double Q2)
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4int HadronType[NHADRONS]
static const G4int NHADRONS
const G4String & GetParticleName() const
G4double HadrNucDifferCrSec(G4int Z, G4int Nucleus, G4double Q2)
static const G4int NENERGY
static G4Proton * Proton()
G4double GetQ2_2(G4int N, G4double *Q, G4double *F, G4double R)
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGMass() const
G4double G4Log(G4double x)
void InterpolateHN(G4int n, const G4double EnP[], const G4double C0P[], const G4double C1P[], const G4double B0P[], const G4double B1P[])
static constexpr double proton_mass_c2
static constexpr double amu_c2
G4int HadronType1[NHADRONS]
#define G4MUTEX_INITIALIZER
G4double GetDistrFun(G4double Q2)
#define G4MUTEXINIT(mutex)
G4int HadronCode[NHADRONS]
static constexpr double twopi
G4double TableCrossSec[NQTABLE]
G4double HadronProtonQ2(const G4ParticleDefinition *aHadron, G4double inLabMom)
double A(double temperature)
G4ElasticData(const G4ParticleDefinition *h, G4int Z, G4double A, G4double *eGeV)
G4double SetBinom[240][240]
G4double GetLightFq2(G4int Z, G4int A, G4double Q)
G4double HadronNucleusQ2_2(G4ElasticData *pElD, G4int Z, G4double plabGeV, G4double tmax)
void GetKinematics(const G4ParticleDefinition *aHadron, G4double MomentumH)
static G4Mutex eldata_m[NHADRONS][ZMAX]
G4double LowEdgeEnergy[NENERGY]
G4ElasticHadrNucleusHE(const G4String &name="hElasticGlauber")
G4int GetPDGEncoding() const
G4GLOB_DLL std::ostream G4cout
G4double GetQ2(G4double Ran)
G4NistManager * nistManager
static constexpr double pi
void DefineHadronValues(G4int Z)
virtual void ModelDescription(std::ostream &) const
static constexpr double GeV
G4double SampleT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
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 lowestEnergyLimit
static G4NistManager * Instance()