Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4ElasticHadrNucleusHE Class Reference

#include <G4ElasticHadrNucleusHE.hh>

+ Inheritance diagram for G4ElasticHadrNucleusHE:

Public Member Functions

 G4ElasticHadrNucleusHE (const G4String &name="hElasticGlauber")
 
virtual ~G4ElasticHadrNucleusHE ()
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual void Description () const
 
G4double SampleT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double HadronNucleusQ2_2 (G4ElasticData *pElD, G4int Z, G4double plabGeV, G4double tmax)
 
void DefineHadronValues (G4int Z)
 
G4double GetLightFq2 (G4int Z, G4int A, G4double Q)
 
G4double GetHeavyFq2 (G4int Z, G4int Nucleus, G4double *LineFq2)
 
G4double GetQ2_2 (G4int N, G4double *Q, G4double *F, G4double R)
 
G4double LineInterpol (G4double p0, G4double p2, G4double c1, G4double c2, G4double p)
 
G4double HadrNucDifferCrSec (G4int Z, G4int Nucleus, G4double Q2)
 
void InterpolateHN (G4int n, const G4double EnP[], const G4double C0P[], const G4double C1P[], const G4double B0P[], const G4double B1P[])
 
G4ElasticHadrNucleusHE & operator= (const G4ElasticHadrNucleusHE &right)
 
 G4ElasticHadrNucleusHE (const G4ElasticHadrNucleusHE &)
 
G4double GetBinomCof (G4int n, G4int m)
 
G4double GetFt (G4double Q2)
 
G4double GetDistrFun (G4double Q2)
 
G4double GetQ2 (G4double Ran)
 
G4double HadronProtonQ2 (const G4ParticleDefinition *aHadron, G4double inLabMom)
 
void GetKinematics (const G4ParticleDefinition *aHadron, G4double MomentumH)
 
- Public Member Functions inherited from G4HadronElastic
 G4HadronElastic (const G4String &name="hElasticLHEP")
 
virtual ~G4HadronElastic ()
 
virtual G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
void SetLowestEnergyLimit (G4double value)
 
G4double LowestEnergyLimit () const
 
G4double ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual void Description () const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &, G4Nucleus &)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
const G4HadronicInteraction * GetMyPointer () const
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4String & GetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
G4bool operator== (const G4HadronicInteraction &right) const
 
G4bool operator!= (const G4HadronicInteraction &right) const
 
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void ModelDescription (std::ostream &outFile) const
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 108 of file G4ElasticHadrNucleusHE.hh.

Constructor & Destructor Documentation

◆ G4ElasticHadrNucleusHE() [1/2]

G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE ( const G4String &  name = "hElasticGlauber")

Definition at line 226 of file G4ElasticHadrNucleusHE.cc.

227 : G4HadronElastic(name)
228{
229 dQ2 = hMass = hMass2 = hLabMomentum = hLabMomentum2 = MomentumCM = HadrEnergy
230 = R1 = R2 = Pnucl = Aeff = HadrTot = HadrSlope = HadrReIm = TotP = DDSect2
231 = DDSect3 = ConstU = FmaxT = Slope1 = Slope2 = Coeff1 = Coeff2 = MaxTR
232 = Slope0 = Coeff0 = aAIm = aDIm = Dtot11 = 0.0;
233 NumbN = iHadrCode = iHadron = 0;
234
235 verboseLevel = 0;
236 plabLowLimit = 20.0*MeV;
237 lowestEnergyLimit = 0.0;
238 //Description();
239
240 MbToGeV2 = 2.568;
241 sqMbToGeV = 1.602;
242 Fm2ToGeV2 = 25.68;
243 GeV2 = GeV*GeV;
244 protonM = proton_mass_c2/GeV;
245 protonM2 = protonM*protonM;
246
247 BoundaryP[0]=9.0;BoundaryTG[0]=5.0;BoundaryTL[0]=0.;
248 BoundaryP[1]=20.0;BoundaryTG[1]=1.5;BoundaryTL[1]=0.;
249 BoundaryP[2]=5.0; BoundaryTG[2]=1.0;BoundaryTL[2]=1.5;
250 BoundaryP[3]=8.0; BoundaryTG[3]=3.0;BoundaryTL[3]=0.;
251 BoundaryP[4]=7.0; BoundaryTG[4]=3.0;BoundaryTL[4]=0.;
252 BoundaryP[5]=5.0; BoundaryTG[5]=2.0;BoundaryTL[5]=0.;
253 BoundaryP[6]=5.0; BoundaryTG[6]=1.5;BoundaryTL[6]=3.0;
254
255 Binom();
256 // energy in GeV
257 Energy[0] = 0.4;
258 Energy[1] = 0.6;
259 Energy[2] = 0.8;
260 LowEdgeEnergy[0] = 0.0;
261 LowEdgeEnergy[1] = 0.5;
262 LowEdgeEnergy[2] = 0.7;
263 G4double e = 1.0;
264 G4double f = std::pow(10.,0.1);
265 for(G4int i=3; i<NENERGY; i++) {
266 Energy[i] = e;
267 LowEdgeEnergy[i] = e/f;
268 e *= f*f;
269 }
270 nistManager = G4NistManager::Instance();
271
272 // PDG code for hadrons
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};
276 // internal index
277 G4int id[NHADRONS] = {2,3,6,0,4,5,4,4,4,5,
278 0,0,0,0,0,0,0,
279 1,7,1,1,1,1,1,1,1};
280
281 G4int id1[NHADRONS] = {3,4,1,0,5,6,5,5,5,6,
282 0,0,0,0,0,0,0,
283 2,2,2,2,2,2,2,2,2};
284
285 for(G4int j=0; j<NHADRONS; j++)
286 {
287 HadronCode[j] = ic[j];
288 HadronType[j] = id[j];
289 HadronType1[j] = id1[j];
290
291 for(G4int k = 0; k < 93; k++) { SetOfElasticData[j][k] = 0; }
292 }
293}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68

◆ ~G4ElasticHadrNucleusHE()

G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusHE ( )
virtual

Definition at line 327 of file G4ElasticHadrNucleusHE.cc.

328{
329 for(G4int j = 0; j < NHADRONS; j++)
330 {
331 for(G4int k = 0; k < 93; k++)
332 {
333 if(SetOfElasticData[j][k]) delete SetOfElasticData[j][k];
334 }
335 }
336}

◆ G4ElasticHadrNucleusHE() [2/2]

G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE ( const G4ElasticHadrNucleusHE &  )

Member Function Documentation

◆ DefineHadronValues()

void G4ElasticHadrNucleusHE::DefineHadronValues ( G4int  Z)

Definition at line 948 of file G4ElasticHadrNucleusHE.cc.

949{
950 HadrEnergy = std::sqrt(hMass2 + hLabMomentum2);
951
952 G4double sHadr = 2.*HadrEnergy*protonM+protonM2+hMass2;
953 G4double sqrS = std::sqrt(sHadr);
954 G4double Ecm = 0.5*(sHadr-hMass2+protonM2)/sqrS;
955 MomentumCM = std::sqrt(Ecm*Ecm-protonM2);
956
957 if(verboseLevel>2)
958 G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron
959 << " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS
960 << " plab= " << hLabMomentum
961 <<" E - m "<<HadrEnergy - hMass<< G4endl;
962
963 G4double TotN = 0.0;
964 G4double logE = std::log(HadrEnergy);
965 G4double logS = std::log(sHadr);
966 TotP = 0.0;
967
968 switch (iHadron)
969 {
970 case 0: // proton, neutron
971 case 6:
972
973 if(hLabMomentum > 10)
974 TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165); // mb
975
976 else
977 {
978// ================== neutron ================
979
980//// if(iHadrCode == 2112)
981
982
983 if( hLabMomentum > 1.4 )
984 TotN = 33.3+15.2*(hLabMomentum2-1.35)/
985 (std::pow(hLabMomentum,2.37)+0.95);
986
987 else if(hLabMomentum > 0.8)
988 {
989 G4double A0 = logE + 0.0513;
990 TotN = 33.0 + 25.5*A0*A0;
991 }
992 else
993 {
994 G4double A0 = logE - 0.2634; // log(1.3)
995 TotN = 33.0 + 30.*A0*A0*A0*A0;
996 }
997// ================= proton ===============
998// else if(iHadrCode == 2212)
999 {
1000 if(hLabMomentum >= 1.05)
1001 {
1002 TotP = 39.0+75.*(hLabMomentum-1.2)/
1003 (hLabMomentum2*hLabMomentum+0.15);
1004 }
1005
1006 else if(hLabMomentum >= 0.7)
1007 {
1008 G4double A0 = logE + 0.3147;
1009 TotP = 23.0 + 40.*A0*A0;
1010 }
1011 else
1012 {
1013 TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5);
1014 }
1015 }
1016 }
1017
1018// HadrTot = 0.5*(82*TotP+126*TotN)/104; // $$$$$$$$$$$$$$$$$$
1019 HadrTot = 0.5*(TotP+TotN);
1020// ...................................................
1021 // Proton slope
1022 if(hLabMomentum >= 2.) HadrSlope = 5.44 + 0.88*logS;
1023
1024 else if(hLabMomentum >= 0.5) HadrSlope = 3.73*hLabMomentum-0.37;
1025
1026 else HadrSlope = 1.5;
1027
1028// ...................................................
1029 if(hLabMomentum >= 1.2)
1030 HadrReIm = 0.13*(logS - 5.8579332)*std::pow(sHadr,-0.18);
1031
1032 else if(hLabMomentum >= 0.6)
1033 HadrReIm = -75.5*(std::pow(hLabMomentum,0.25)-0.95)/
1034 (std::pow(3*hLabMomentum,2.2)+1);
1035
1036 else
1037 HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2);
1038// ...................................................
1039 DDSect2 = 2.2; //mb*GeV-2
1040 DDSect3 = 0.6; //mb*GeV-2
1041 // ================== lambda ==================
1042 if( iHadrCode == 3122)
1043 {
1044 HadrTot *= 0.88;
1045 HadrSlope *=0.85;
1046 }
1047 // ================== sigma + ==================
1048 else if( iHadrCode == 3222)
1049 {
1050 HadrTot *=0.81;
1051 HadrSlope *=0.85;
1052 }
1053 // ================== sigma 0,- ==================
1054 else if(iHadrCode == 3112 || iHadrCode == 3212 )
1055 {
1056 HadrTot *=0.88;
1057 HadrSlope *=0.85;
1058 }
1059 // =================== xi =================
1060 else if( iHadrCode == 3312 || iHadrCode == 3322 )
1061 {
1062 HadrTot *=0.77;
1063 HadrSlope *=0.75;
1064 }
1065 // ================= omega =================
1066 else if( iHadrCode == 3334)
1067 {
1068 HadrTot *=0.78;
1069 HadrSlope *=0.7;
1070 }
1071
1072 break;
1073// ===========================================================
1074 case 1: // antiproton
1075 case 7: // antineutron
1076
1077 HadrTot = 5.2+5.2*logE + 123.2/sqrS; // mb
1078 HadrSlope = 8.32+0.57*logS; //(GeV/c)^-2
1079
1080 if( HadrEnergy < 1000 )
1081 HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8);
1082 else
1083 HadrReIm = 0.6*(logS - 5.8579332)*std::pow(sHadr,-0.25);
1084
1085 DDSect2 = 11; //mb*(GeV/c)^-2
1086 DDSect3 = 3; //mb*(GeV/c)^-2
1087 // ================== lambda ==================
1088 if( iHadrCode == -3122)
1089 {
1090 HadrTot *= 0.88;
1091 HadrSlope *=0.85;
1092 }
1093 // ================== sigma + ==================
1094 else if( iHadrCode == -3222)
1095 {
1096 HadrTot *=0.81;
1097 HadrSlope *=0.85;
1098 }
1099 // ================== sigma 0,- ==================
1100 else if(iHadrCode == -3112 || iHadrCode == -3212 )
1101 {
1102 HadrTot *=0.88;
1103 HadrSlope *=0.85;
1104 }
1105 // =================== xi =================
1106 else if( iHadrCode == -3312 || iHadrCode == -3322 )
1107 {
1108 HadrTot *=0.77;
1109 HadrSlope *=0.75;
1110 }
1111 // ================= omega =================
1112 else if( iHadrCode == -3334)
1113 {
1114 HadrTot *=0.78;
1115 HadrSlope *=0.7;
1116 }
1117
1118 break;
1119// -------------------------------------------
1120 case 2: // pi plus, pi minus
1121 case 3:
1122
1123 if(hLabMomentum >= 3.5)
1124 TotP = 10.6+2.*logE + 25.*std::pow(HadrEnergy,-0.43); // mb
1125// =========================================
1126 else if(hLabMomentum >= 1.15)
1127 {
1128 G4double x = (hLabMomentum - 2.55)/0.55;
1129 G4double y = (hLabMomentum - 1.47)/0.225;
1130 TotP = 3.2*std::exp(-x*x) + 12.*std::exp(-y*y) + 27.5;
1131 }
1132// =========================================
1133 else if(hLabMomentum >= 0.4)
1134 {
1135 TotP = 88*(logE+0.2877)*(logE+0.2877)+14.0;
1136 }
1137// =========================================
1138 else
1139 {
1140 G4double x = (hLabMomentum - 0.29)/0.085;
1141 TotP = 20. + 180.*std::exp(-x*x);
1142 }
1143// -------------------------------------------
1144
1145 if(hLabMomentum >= 3.0 )
1146 TotN = 10.6 + 2.*logE + 30.*std::pow(HadrEnergy,-0.43); // mb
1147
1148 else if(hLabMomentum >= 1.3)
1149 {
1150 G4double x = (hLabMomentum - 2.1)/0.4;
1151 G4double y = (hLabMomentum - 1.4)/0.12;
1152 TotN = 36.1+0.079 - 4.313*logE + 3.*std::exp(-x*x) +
1153 1.5*std::exp(-y*y);
1154 }
1155 else if(hLabMomentum >= 0.65)
1156 {
1157 G4double x = (hLabMomentum - 0.72)/0.06;
1158 G4double y = (hLabMomentum - 1.015)/0.075;
1159 TotN = 36.1 + 10.*std::exp(-x*x) + 24*std::exp(-y*y);
1160 }
1161 else if(hLabMomentum >= 0.37)
1162 {
1163 G4double x = std::log(hLabMomentum/0.48);
1164 TotN = 26. + 110.*x*x;
1165 }
1166 else
1167 {
1168 G4double x = (hLabMomentum - 0.29)/0.07;
1169 TotN = 28.0 + 40.*std::exp(-x*x);
1170 }
1171 HadrTot = (TotP+TotN)/2;
1172// ........................................
1173 HadrSlope = 7.28+0.245*logS; // GeV-2
1174 HadrReIm = 0.2*(logS - 4.6051702)*std::pow(sHadr,-0.15);
1175
1176 DDSect2 = 0.7; //mb*GeV-2
1177 DDSect3 = 0.27; //mb*GeV-2
1178
1179 break;
1180// ==========================================================
1181 case 4: // K plus
1182
1183 HadrTot = 10.6+1.8*logE + 9.0*std::pow(HadrEnergy,-0.55); // mb
1184 if(HadrEnergy>100) HadrSlope = 15.0;
1185 else HadrSlope = 1.0+1.76*logS - 2.84/sqrS; // GeV-2
1186
1187 HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1);
1188 DDSect2 = 0.7; //mb*GeV-2
1189 DDSect3 = 0.21; //mb*GeV-2
1190 break;
1191// =========================================================
1192 case 5: // K minus
1193
1194 HadrTot = 10+1.8*logE + 25./sqrS; // mb
1195 HadrSlope = 6.98+0.127*logS; // GeV-2
1196 HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1);
1197 DDSect2 = 0.7; //mb*GeV-2
1198 DDSect3 = 0.27; //mb*GeV-2
1199 break;
1200 }
1201// =========================================================
1202 if(verboseLevel>2)
1203 G4cout << "HadrTot= " << HadrTot << " HadrSlope= " << HadrSlope
1204 << " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2
1205 << " DDSect3= " << DDSect3 << G4endl;
1206
1207 if(Z != 1) return;
1208
1209 // Scattering of protons
1210
1211 Coeff0 = Coeff1 = Coeff2 = 0.0;
1212 Slope0 = Slope1 = 1.0;
1213 Slope2 = 5.0;
1214
1215 // data for iHadron=0
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};
1221
1222 // data for iHadron=6,7
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};
1228
1229 // data for iHadron=1
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};
1235
1236 // data for iHadron=2
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};
1242
1243 // data for iHadron=3
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};
1249
1250 // data for iHadron=4
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};
1256
1257 // data for iHadron=5
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};
1263
1264 switch(iHadron)
1265 {
1266 case 0 :
1267
1268 if(hLabMomentum <BoundaryP[0])
1269 InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0);
1270
1271 Coeff2 = 0.8/hLabMomentum/hLabMomentum;
1272 break;
1273
1274 case 6 :
1275
1276 if(hLabMomentum < BoundaryP[1])
1277 InterpolateHN(5,EnN,C0N,C1N,B0N,B1N);
1278
1279 Coeff2 = 0.8/hLabMomentum/hLabMomentum;
1280 break;
1281
1282 case 1 :
1283 case 7 :
1284 if(hLabMomentum < BoundaryP[2])
1285 InterpolateHN(2,EnP,C0P,C1P,B0P,B1P);
1286 break;
1287
1288 case 2 :
1289
1290 if(hLabMomentum < BoundaryP[3])
1291 InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1PP);
1292
1293 Coeff2 = 0.02/hLabMomentum;
1294 break;
1295
1296 case 3 :
1297
1298 if(hLabMomentum < BoundaryP[4])
1299 InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN,B1PPN);
1300
1301 Coeff2 = 0.02/hLabMomentum;
1302 break;
1303
1304 case 4 :
1305
1306 if(hLabMomentum < BoundaryP[5])
1307 InterpolateHN(4,EnK,C0K,C1K,B0K,B1K);
1308
1309 if(hLabMomentum < 1) Coeff2 = 0.34;
1310 else Coeff2 = 0.34/hLabMomentum2/hLabMomentum;
1311 break;
1312
1313 case 5 :
1314 if(hLabMomentum < BoundaryP[6])
1315 InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1KM);
1316
1317 if(hLabMomentum < 1) Coeff2 = 0.01;
1318 else Coeff2 = 0.01/hLabMomentum2/hLabMomentum;
1319 break;
1320 }
1321
1322 if(verboseLevel > 2)
1323 G4cout<<" HadrVal : Plasb "<<hLabMomentum
1324 <<" iHadron "<<iHadron<<" HadrTot "<<HadrTot<<G4endl;
1325}
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4ElasticHadrNucleusHE::InterpolateHN
void InterpolateHN(G4int n, const G4double EnP[], const G4double C0P[], const G4double C1P[], const G4double B0P[], const G4double B1P[])
Definition: G4ElasticHadrNucleusHE.hh:247

Referenced by GetKinematics(), and HadronNucleusQ2_2().

◆ Description()

void G4ElasticHadrNucleusHE::Description ( ) const
virtual

Reimplemented from G4HadronElastic.

Definition at line 296 of file G4ElasticHadrNucleusHE.cc.

297{
298 char* dirName = getenv("G4PhysListDocDir");
299 if (dirName) {
300 std::ofstream outFile;
301 G4String outFileName = GetModelName() + ".html";
302 G4String pathName = G4String(dirName) + "/" + outFileName;
303 outFile.open(pathName);
304 outFile << "<html>\n";
305 outFile << "<head>\n";
306
307 outFile << "<title>Description of G4ElasticHadrNucleusHE Model</title>\n";
308 outFile << "</head>\n";
309 outFile << "<body>\n";
310
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";
315
316 outFile << "</body>\n";
317 outFile << "</html>\n";
318 outFile.close();
319 }
320}
G4HadronicInteraction::GetModelName
const G4String & GetModelName() const
Definition: G4HadronicInteraction.hh:119

◆ GetBinomCof()

G4double G4ElasticHadrNucleusHE::GetBinomCof ( G4int  n,
G4int  m 
)
inline

Definition at line 269 of file G4ElasticHadrNucleusHE.hh.

270{
271 if ( numN >= numM && numN <= 240) return SetBinom[numN][numM];
272 else return 0.;
273}

◆ GetDistrFun()

G4double G4ElasticHadrNucleusHE::GetDistrFun ( G4double  Q2)
inline

Definition at line 278 of file G4ElasticHadrNucleusHE.hh.

279{
280 return GetFt(Q2)/FmaxT;
281}

Referenced by GetQ2().

◆ GetFt()

G4double G4ElasticHadrNucleusHE::GetFt ( G4double  Q2)

Definition at line 1369 of file G4ElasticHadrNucleusHE.cc.

1370{
1371 G4double Fdistr=0;
1372 G4double SqrQ2 = std::sqrt(Q2);
1373
1374 Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*std::exp(ConstU))
1375 /HadrSlope*(1-std::exp(-HadrSlope*Q2))
1376 + Coeff0*(1-std::exp(-Slope0*Q2))
1377 + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1)
1378 + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
1379
1380 if (verboseLevel>1)
1381 G4cout<<"Old: Coeff0 Coeff1 Coeff2 "<<Coeff0<<" "
1382 <<Coeff1<<" "<<Coeff2<<" Slope Slope0 Slope1 Slope2 "
1383 <<HadrSlope<<" "<<Slope0<<" "<<Slope1<<" "<<Slope2
1384 <<" Fdistr "<<Fdistr<<G4endl;
1385 return Fdistr;
1386}

Referenced by GetDistrFun(), and GetQ2().

◆ GetHeavyFq2()

G4double G4ElasticHadrNucleusHE::GetHeavyFq2 ( G4int  Z,
G4int  Nucleus,
G4double *  LineFq2 
)

Definition at line 596 of file G4ElasticHadrNucleusHE.cc.

597{
598 G4int ii, jj, aSimp;
599 G4double curQ2, curSec;
600 G4double curSum = 0.0;
601 G4double totSum = 0.0;
602
603 G4double ddQ2 = dQ2/20;
604 G4double Q2l = 0;
605
606 LineF[0] = 0;
607 for(ii = 1; ii<ONQ2; ii++)
608 {
609 curSum = 0;
610 aSimp = 4;
611
612 for(jj = 0; jj<20; jj++)
613 {
614 curQ2 = Q2l+jj*ddQ2;
615
616 curSec = HadrNucDifferCrSec(Z, Nucleus, curQ2);
617 curSum += curSec*aSimp;
618
619 if(aSimp > 3) aSimp = 2;
620 else aSimp = 4;
621
622 if(jj == 0 && verboseLevel>2)
623 G4cout<<" Q2 "<<curQ2<<" AIm "<<aAIm<<" DIm "<<aDIm
624 <<" Diff "<<curSec<<" totSum "<<totSum<<G4endl;
625 }
626
627 Q2l += dQ2;
628 curSum *= ddQ2/2.3; // $$$$$$$$$$$$$$$$$$$$$$$
629 totSum += curSum;
630
631 LineF[ii] = totSum;
632
633 if (verboseLevel>2)
634 G4cout<<" GetHeavy: Q2 dQ2 totSum "<<Q2l<<" "<<dQ2<<" "<<totSum
635 <<" curSec "
636 <<curSec<<" totSum "<< totSum<<" DTot "
637 <<curSum<<G4endl;
638 }
639 return totSum;
640}
G4ElasticHadrNucleusHE::HadrNucDifferCrSec
G4double HadrNucDifferCrSec(G4int Z, G4int Nucleus, G4double Q2)
Definition: G4ElasticHadrNucleusHE.cc:759

Referenced by HadronNucleusQ2_2().

◆ GetKinematics()

void G4ElasticHadrNucleusHE::GetKinematics ( const G4ParticleDefinition *  aHadron,
G4double  MomentumH 
)

Definition at line 1328 of file G4ElasticHadrNucleusHE.cc.

1331{
1332 if (verboseLevel>1)
1333 G4cout<<"1 GetKin.: HadronName MomentumH "
1334 <<aHadron->GetParticleName()<<" "<<MomentumH<<G4endl;
1335
1336 DefineHadronValues(1);
1337
1338 G4double Sh = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV
1339
1340 ConstU = 2*protonM2+2*hMass2-Sh;
1341
1342 G4double MaxT = 4*MomentumCM*MomentumCM;
1343
1344 BoundaryTL[0] = MaxT; //2.0;
1345 BoundaryTL[1] = MaxT;
1346 BoundaryTL[3] = MaxT;
1347 BoundaryTL[4] = MaxT;
1348 BoundaryTL[5] = MaxT;
1349
1350 G4int NumberH=0;
1351
1352 while(iHadrCode!=HadronCode[NumberH]) NumberH++;
1353
1354 NumberH = HadronType1[NumberH];
1355
1356 if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH];
1357 else MaxTR = BoundaryTG[NumberH];
1358
1359 if (verboseLevel>1)
1360 G4cout<<"3 GetKin. : NumberH "<<NumberH
1361 <<" Bound.P[NumberH] "<<BoundaryP[NumberH]
1362 <<" Bound.TL[NumberH] "<<BoundaryTL[NumberH]
1363 <<" Bound.TG[NumberH] "<<BoundaryTG[NumberH]
1364 <<" MaxT MaxTR "<<MaxT<<" "<<MaxTR<<G4endl;
1365
1366// GetParametersHP(aHadron, MomentumH);
1367}
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:115

Referenced by HadronProtonQ2().

◆ GetLightFq2()

G4double G4ElasticHadrNucleusHE::GetLightFq2 ( G4int  Z,
G4int  A,
G4double  Q 
)

Definition at line 646 of file G4ElasticHadrNucleusHE.cc.

648{
649 // Scattering of proton
650 if(Z == 1)
651 {
652 G4double SqrQ2 = std::sqrt(Q2);
653 G4double valueConstU = 2.*(hMass2 + protonM2) - Q2;
654
655 G4double y = (1.-Coeff1-Coeff0)/HadrSlope*(1.-std::exp(-HadrSlope*Q2))
656 + Coeff0*(1.-std::exp(-Slope0*Q2))
657 + Coeff2/Slope2*std::exp(Slope2*valueConstU)*(std::exp(Slope2*Q2)-1.)
658 + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
659
660 return y;
661 }
662
663 // The preparing of probability function
664
665 G4double prec = Nucleus > 208 ? 1.0e-7 : 1.0e-6;
666
667 G4double Stot = HadrTot*MbToGeV2; // Gev^-2
668 G4double Bhad = HadrSlope; // GeV^-2
669 G4double Asq = 1+HadrReIm*HadrReIm;
670 G4double Rho2 = std::sqrt(Asq);
671
672// if(verboseLevel >1)
673 G4cout<<" Fq2 Before for i Tot B Im "<<HadrTot<<" "<<HadrSlope<<" "
674 <<HadrReIm<<G4endl;
675
676 if(verboseLevel > 1) {
677 G4cout << "GetFq2: Stot= " << Stot << " Bhad= " << Bhad
678 <<" Im "<<HadrReIm
679 << " Asq= " << Asq << G4endl;
680 G4cout << "R1= " << R1 << " R2= " << R2 << " Pnucl= " << Pnucl <<G4endl;
681 }
682 G4double R12 = R1*R1;
683 G4double R22 = R2*R2;
684 G4double R12B = R12+2*Bhad;
685 G4double R22B = R22+2*Bhad;
686
687 G4double Norm = (R12*R1-Pnucl*R22*R2); // HP->Aeff;
688
689 G4double R13 = R12*R1/R12B;
690 G4double R23 = Pnucl*R22*R2/R22B;
691 G4double Unucl = Stot/twopi/Norm*R13;
692 G4double UnucRho2 = -Unucl*Rho2;
693
694 G4double FiH = std::asin(HadrReIm/Rho2);
695 G4double NN2 = R23/R13;
696
697 if(verboseLevel > 2)
698 G4cout << "UnucRho2= " << UnucRho2 << " FiH= " << FiH << " NN2= " << NN2
699 << " Norm= " << Norm << G4endl;
700
701 G4double dddd;
702
703 G4double Prod0 = 0;
704 G4double N1 = -1.0;
705 //G4double Tot0 = 0;
706 G4double exp1;
707
708 G4double Prod3 ;
709 G4double exp2 ;
710 G4double N4, N5, N2, Prod1, Prod2;
711 G4int i1, i2, j1, j2;
712
713 for(i1 = 1; i1<= Nucleus; i1++) ////++++++++++ i1
714 {
715 N1 = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2;
716 Prod1 = 0;
717 //Tot0 = 0;
718 N2 = -1;
719
720 for(i2 = 1; i2<=Nucleus; i2++) ////+++++++++ i2
721 {
722 N2 = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2;
723 Prod2 = 0;
724 N5 = -1/NN2;
725 for(j2=0; j2<= i2; j2++) ////+++++++++ j2
726 {
727 Prod3 = 0;
728 exp2 = 1/(j2/R22B+(i2-j2)/R12B);
729 N5 = -N5*NN2;
730 N4 = -1/NN2;
731 for(j1=0; j1<=i1; j1++) ////++++++++ j1
732 {
733 exp1 = 1/(j1/R22B+(i1-j1)/R12B);
734 dddd = exp1+exp2;
735 N4 = -N4*NN2;
736 Prod3 = Prod3+N4*exp1*exp2*
737 (1-std::exp(-Q2*dddd/4))/dddd*4*SetBinom[i1][j1];
738 } // j1
739 Prod2 = Prod2 +Prod3*N5*SetBinom[i2][j2];
740 } // j2
741 Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2));
742
743 if (std::fabs(Prod2*N2/Prod1)<prec) break;
744 } // i2
745 Prod0 = Prod0 + Prod1*N1;
746 if(std::fabs(N1*Prod1/Prod0) < prec) break;
747
748 } // i1
749
750 Prod0 *= 0.25*pi/MbToGeV2; // This is in mb
751
752 if(verboseLevel>1)
753 G4cout << "GetLightFq2 Z= " << Z << " A= " << Nucleus
754 <<" Q2= " << Q2 << " Res= " << Prod0 << G4endl;
755 return Prod0;
756}
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:70

◆ GetQ2()

G4double G4ElasticHadrNucleusHE::GetQ2 ( G4double  Ran)

Definition at line 1388 of file G4ElasticHadrNucleusHE.cc.

1389{
1390 G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta;
1391 G4double Q2=0;
1392
1393 FmaxT = GetFt(MaxTR);
1394 delta = GetDistrFun(DDD0)-Ran;
1395
1396 while(std::fabs(delta) > 0.0001)
1397 {
1398 if(delta>0)
1399 {
1400 DDD2 = DDD0;
1401 DDD0 = (DDD0+DDD1)*0.5;
1402 }
1403 else if(delta<0)
1404 {
1405 DDD1 = DDD0;
1406 DDD0 = (DDD0+DDD2)*0.5;
1407 }
1408 delta = GetDistrFun(DDD0)-Ran;
1409 }
1410
1411 Q2 = DDD0;
1412
1413 return Q2;
1414}
G4ElasticHadrNucleusHE::GetDistrFun
G4double GetDistrFun(G4double Q2)
Definition: G4ElasticHadrNucleusHE.hh:278

Referenced by HadronProtonQ2().

◆ GetQ2_2()

G4double G4ElasticHadrNucleusHE::GetQ2_2 ( G4int  N,
G4double *  Q,
G4double *  F,
G4double  R 
)

Definition at line 542 of file G4ElasticHadrNucleusHE.cc.

544{
545 G4double ranQ2;
546
547 G4double F1 = F[kk-2];
548 G4double F2 = F[kk-1];
549 G4double F3 = F[kk];
550 G4double X1 = Q[kk-2];
551 G4double X2 = Q[kk-1];
552 G4double X3 = Q[kk];
553
554 if(verboseLevel > 2)
555 G4cout << "GetQ2_2 kk= " << kk << " X2= " << X2 << " X3= " << X3
556 << " F2= " << F2 << " F3= " << F3 << " Rndm= " << ranUni << G4endl;
557
558 if(kk == 1 || kk == 0)
559 {
560 F1 = F[0];
561 F2 = F[1];
562 F3 = F[2];
563 X1 = Q[0];
564 X2 = Q[1];
565 X3 = Q[2];
566 }
567
568 G4double F12 = F1*F1;
569 G4double F22 = F2*F2;
570 G4double F32 = F3*F3;
571
572 G4double D0 = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3;
573
574 if(verboseLevel > 2)
575 G4cout << " X1= " << X1 << " F1= " << F1 << " D0= "
576 << D0 << G4endl;
577
578 if(std::abs(D0) < 0.00000001)
579 {
580 ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
581 }
582 else
583 {
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;
589 }
590 return ranQ2; // MeV^2
591}
F22
#define F22
Definition: G4ErrorMatrix.cc:901
F32
#define F32
Definition: G4ErrorMatrix.cc:906
F12
#define F12
Definition: G4ErrorMatrix.cc:896

Referenced by HadronNucleusQ2_2().

◆ HadrNucDifferCrSec()

G4double G4ElasticHadrNucleusHE::HadrNucDifferCrSec ( G4int  Z,
G4int  Nucleus,
G4double  Q2 
)

GeV/GeV;

Definition at line 758 of file G4ElasticHadrNucleusHE.cc.

760{
761// ------ All external kinematical variables are in MeV -------
762// ------ but internal in GeV !!!! ------
763
764 G4int NWeight = int( nistManager->GetAtomicMassAmu(Z) + 0.5 );
765
766 G4double theQ2 = aQ2; ///GeV/GeV;
767
768 // Scattering of proton
769 if(NWeight == 1)
770 {
771 G4double SqrQ2 = std::sqrt(aQ2);
772 G4double valueConstU = hMass2 + protonM2-2*protonM*HadrEnergy - aQ2;
773
774 G4double MaxT = 4*MomentumCM*MomentumCM;
775
776 BoundaryTL[0] = MaxT;
777 BoundaryTL[1] = MaxT;
778 BoundaryTL[3] = MaxT;
779 BoundaryTL[4] = MaxT;
780 BoundaryTL[5] = MaxT;
781
782 G4double dSigPodT;
783
784 dSigPodT = HadrTot*HadrTot*(1+HadrReIm*HadrReIm)*
785 (
786 Coeff1*std::exp(-Slope1*SqrQ2)+
787 Coeff2*std::exp( Slope2*(valueConstU)+aQ2)+
788 (1-Coeff1-Coeff0)*std::exp(-HadrSlope*aQ2)+
789 +Coeff0*std::exp(-Slope0*aQ2)
790// +0.1*(1-std::fabs(CosTh))
791 )/16/3.1416*2.568;
792
793 return dSigPodT;
794 }
795
796 G4double Stot = HadrTot*MbToGeV2;
797 G4double Bhad = HadrSlope;
798 G4double Asq = 1+HadrReIm*HadrReIm;
799 G4double Rho2 = std::sqrt(Asq);
800 G4double Pnuclp = 0.001;
801 Pnuclp = Pnucl;
802 G4double R12 = R1*R1;
803 G4double R22 = R2*R2;
804 G4double R12B = R12+2*Bhad;
805 G4double R22B = R22+2*Bhad;
806 G4double R12Ap = R12+20;
807 G4double R22Ap = R22+20;
808 G4double R13Ap = R12*R1/R12Ap;
809 G4double R23Ap = R22*R2/R22Ap*Pnuclp;
810 G4double R23dR13 = R23Ap/R13Ap;
811 G4double R12Apd = 2/R12Ap;
812 G4double R22Apd = 2/R22Ap;
813 G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);
814
815 G4double DDSec1p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R1/4));
816 G4double DDSec2p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/
817 std::sqrt((R12+R22)/2)/4));
818 G4double DDSec3p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R2/4));
819
820 G4double Norm = (R12*R1-Pnucl*R22*R2)*Aeff;
821 G4double Normp = (R12*R1-Pnuclp*R22*R2)*Aeff;
822 G4double R13 = R12*R1/R12B;
823 G4double R23 = Pnucl*R22*R2/R22B;
824 G4double Unucl = Stot/twopi/Norm*R13;
825 G4double UnuclScr = Stot/twopi/Normp*R13Ap;
826 G4double SinFi = HadrReIm/Rho2;
827 G4double FiH = std::asin(SinFi);
828 G4double N = -1;
829 G4double N2 = R23/R13;
830
831 G4double ImElasticAmpl0 = 0;
832 G4double ReElasticAmpl0 = 0;
833
834 G4double exp1;
835 G4double N4;
836 G4double Prod1, Tot1=0, medTot, DTot1, DmedTot;
837 G4int i;
838
839 for( i=1; i<=NWeight; i++)
840 {
841 N = -N*Unucl*(NWeight-i+1)/i*Rho2;
842 N4 = 1;
843 Prod1 = std::exp(-theQ2/i*R12B/4)/i*R12B;
844 medTot = R12B/i;
845
846 for(G4int l=1; l<=i; l++)
847 {
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;
852 } // end l
853
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;
858 } // i
859
860 ImElasticAmpl0 = ImElasticAmpl0*pi/2.568; // The amplitude in mB
861 ReElasticAmpl0 = ReElasticAmpl0*pi/2.568; // The amplitude in mB
862 Tot1 = Tot1*twopi/2.568;
863
864 G4double C1 = R13Ap*R13Ap/2*DDSec1p;
865 G4double C2 = 2*R23Ap*R13Ap/2*DDSec2p;
866 G4double C3 = R23Ap*R23Ap/2*DDSec3p;
867
868 G4double N1p = 1;
869
870 G4double Din1 = 0.5;
871
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));
875
876 DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2);
877
878 G4double exp1p;
879 G4double exp2p;
880 G4double exp3p;
881 G4double N2p;
882 G4double Din2, BinCoeff;
883
884 BinCoeff = 1;
885
886 for( i = 1; i<= NWeight-2; i++)
887 {
888 N1p = -N1p*UnuclScr*(NWeight-i-1)/i*Rho2;
889 N2p = 1;
890 Din2 = 0;
891 DmedTot = 0;
892 for(G4int l = 0; l<=i; l++)
893 {
894 if(l == 0) BinCoeff = 1;
895 else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l;
896
897 exp1 = l/R22B+(i-l)/R12B;
898 exp1p = exp1+R12Apd;
899 exp2p = exp1+R12ApdR22Ap;
900 exp3p = exp1+R22Apd;
901
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));
906
907 DmedTot = DmedTot + N2p*BinCoeff*
908 (C1/exp1p-C2/exp2p+C3/exp3p);
909
910 N2p = -N2p*R23dR13;
911 } // l
912
913 Din1 = Din1+Din2*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
914 DTot1 = DTot1+DmedTot*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
915
916 if(std::fabs(Din2*N1p/Din1) < 0.000001) break;
917 } // i
918
919 Din1 = -Din1*NWeight*(NWeight-1)
920 /2/pi/Normp/2/pi/Normp*16*pi*pi;
921
922 DTot1 = DTot1*NWeight*(NWeight-1)
923 /2/pi/Normp/2/pi/Normp*16*pi*pi;
924
925 DTot1 *= 5; // $$$$$$$$$$$$$$$$$$$$$$$$
926// Din1 *= 0.2; // %%%%%%%%%%%%%%%%%%%%%%% proton
927// Din1 *= 0.05; // %%%%%%%%%%%%%%%%%%%%%%% pi+
928// ---------------- dSigma/d|-t|, mb/(GeV/c)^-2 -----------------
929
930 G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
931 (ImElasticAmpl0+Din1)*
932 (ImElasticAmpl0+Din1))*2/4/pi;
933
934 Tot1 = Tot1-DTot1;
935 // Tott1 = Tot1*1.0;
936 Dtot11 = DTot1;
937 aAIm = ImElasticAmpl0;
938 aDIm = Din1;
939
940 return DiffCrSec2*1.0; // dSig/d|-t|, mb/(GeV/c)^-2
941} // function
C1
#define C1
C3
#define C3
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:343

Referenced by GetHeavyFq2().

◆ HadronNucleusQ2_2()

G4double G4ElasticHadrNucleusHE::HadronNucleusQ2_2 ( G4ElasticData *  pElD,
G4int  Z,
G4double  plabGeV,
G4double  tmax 
)

Definition at line 441 of file G4ElasticHadrNucleusHE.cc.

444{
445 G4double LineFq2[ONQ2];
446
447 G4double Rand = G4UniformRand();
448
449 G4int iNumbQ2 = 0;
450 G4double Q2 = 0.0;
451
452 G4double ptot2 = plab*plab;
453 G4double ekin = std::sqrt(hMass2 + ptot2) - hMass;
454
455 if(verboseLevel > 1)
456 G4cout<<"Q2_2: ekin plab "<<ekin<<" "<<plab<<" tmax "<<tmax<<G4endl;
457
458 // Find closest energy bin
459 G4int NumbOnE;
460 for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ )
461 {
462 if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break;
463 }
464 G4double* dNumbQ2 = pElD->TableQ2;
465
466 G4int index = NumbOnE*ONQ2;
467
468 // Select kinematics for node energy
469 G4double T = Energy[NumbOnE];
470 hLabMomentum2 = T*(T + 2.*hMass);
471 G4double Q2max = pElD->maxQ2[NumbOnE];
472 G4int length = pElD->dnkE[NumbOnE];
473
474 // Build vector
475 if(length == 0)
476 {
477 R1 = pElD->R1;
478 R2 = pElD->R2;
479 Aeff = pElD->Aeff;
480 Pnucl = pElD->Pnucl;
481 hLabMomentum = std::sqrt(hLabMomentum2);
482
483 DefineHadronValues(Z);
484
485 if(verboseLevel >0)
486 {
487 G4cout<<"1 plab T "<<plab<<" "<<T<<" sigTot B ReIm "
488 <<HadrTot<<" "<<HadrSlope<<" "<<HadrReIm<<G4endl;
489 G4cout<<" R1 R2 Aeff p "<<R1<<" "<<R2<<" "<<Aeff<<" "
490 <<Pnucl<<G4endl;
491 }
492
493 //pElD->CrossSecMaxQ2[NumbOnE] = 1.0;
494
495 if(verboseLevel > 1)
496 G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
497 << " length= " << length
498 << " Q2max= " << Q2max
499 << " ekin= " << ekin <<G4endl;
500
501 pElD->TableCrossSec[index] = 0;
502
503
504 dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0];
505
506 GetHeavyFq2(Z, NumbN, LineFq2); // %%%%%%%%%%%%%%%%%%%%%%%%%
507
508 for(G4int ii=0; ii<ONQ2; ii++)
509 {
510 //if(verboseLevel > 2)
511 // G4cout<<" ii LineFq2 "<<ii<<" "<<LineFq2[ii]/LineFq2[ONQ2-1]
512 // <<" dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl;
513
514 pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1];
515 }
516
517 pElD->dnkE[NumbOnE] = ONQ2;
518 length = ONQ2;
519 }
520
521 G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]);
522
523 for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ )
524 {
525 if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break;
526 }
527 Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
528
529 if(tmax < Q2max) Q2 *= tmax/Q2max;
530
531 if(verboseLevel > 1)
532 G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" iNumbQ2= " << iNumbQ2
533 << " rand= " << Rand << G4endl;
534
535 return Q2;
536}
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4ElasticData::TableQ2
G4double TableQ2[ONQ2]
Definition: G4ElasticHadrNucleusHE.hh:100
G4ElasticData::TableCrossSec
G4double TableCrossSec[NQTABLE]
Definition: G4ElasticHadrNucleusHE.hh:101
G4ElasticData::maxQ2
G4double maxQ2[NENERGY]
Definition: G4ElasticHadrNucleusHE.hh:97
G4ElasticHadrNucleusHE::GetQ2_2
G4double GetQ2_2(G4int N, G4double *Q, G4double *F, G4double R)
Definition: G4ElasticHadrNucleusHE.cc:542
G4ElasticHadrNucleusHE::GetHeavyFq2
G4double GetHeavyFq2(G4int Z, G4int Nucleus, G4double *LineFq2)
Definition: G4ElasticHadrNucleusHE.cc:596

Referenced by SampleInvariantT().

◆ HadronProtonQ2()

G4double G4ElasticHadrNucleusHE::HadronProtonQ2 ( const G4ParticleDefinition *  aHadron,
G4double  inLabMom 
)

Definition at line 1416 of file G4ElasticHadrNucleusHE.cc.

1419{
1420
1421 hMass = p->GetPDGMass()/GeV;
1422 hMass2 = hMass*hMass;
1423 hLabMomentum = inLabMom;
1424 hLabMomentum2 = hLabMomentum*hLabMomentum;
1425 HadrEnergy = sqrt(hLabMomentum2+hMass2);
1426
1427 G4double Rand = G4UniformRand();
1428
1429 GetKinematics(p, inLabMom);
1430
1431 G4double Q2 = GetQ2(Rand);
1432
1433 return Q2;
1434}
G4ElasticHadrNucleusHE::GetQ2
G4double GetQ2(G4double Ran)
Definition: G4ElasticHadrNucleusHE.cc:1388
G4ElasticHadrNucleusHE::GetKinematics
void GetKinematics(const G4ParticleDefinition *aHadron, G4double MomentumH)
Definition: G4ElasticHadrNucleusHE.cc:1329

Referenced by SampleInvariantT().

◆ InterpolateHN()

void G4ElasticHadrNucleusHE::InterpolateHN ( G4int  n,
const G4double  EnP[],
const G4double  C0P[],
const G4double  C1P[],
const G4double  B0P[],
const G4double  B1P[] 
)
inline

Definition at line 247 of file G4ElasticHadrNucleusHE.hh.

250{
251 G4int i;
252
253 for(i=1; i<n; i++) if(hLabMomentum <= EnP[i]) break;
254
255 if(i == n) i = n - 1;
256
257 Coeff0 = LineInterpol(EnP[i], EnP[i-1], C0P[i], C0P[i-1], hLabMomentum);
258 Coeff1 = LineInterpol(EnP[i], EnP[i-1], C1P[i], C1P[i-1], hLabMomentum);
259 Slope0 = LineInterpol(EnP[i], EnP[i-1], B0P[i], B0P[i-1], hLabMomentum);
260 Slope1 = LineInterpol(EnP[i], EnP[i-1], B1P[i], B1P[i-1], hLabMomentum);
261
262// G4cout<<" InterpolHN: n i "<<n<<" "<<i<<" Mom "
263// <<hLabMomentum<<G4endl;
264}
G4ElasticHadrNucleusHE::LineInterpol
G4double LineInterpol(G4double p0, G4double p2, G4double c1, G4double c2, G4double p)
Definition: G4ElasticHadrNucleusHE.hh:233

Referenced by DefineHadronValues().

◆ LineInterpol()

G4double G4ElasticHadrNucleusHE::LineInterpol ( G4double  p0,
G4double  p2,
G4double  c1,
G4double  c2,
G4double  p 
)
inline

Definition at line 233 of file G4ElasticHadrNucleusHE.hh.

236{
237// G4cout<<" LineInterpol: p1 p2 c1 c2 "<<p1<<" "<<p2<<" "
238// <<c1<<" "<<c2<<" c "<<c1+(p-p1)*(c2-c1)/(p2-p1)<<G4endl;
239
240
241 return c1+(p-p1)*(c2-c1)/(p2-p1);
242}

Referenced by InterpolateHN().

◆ operator=()

G4ElasticHadrNucleusHE & G4ElasticHadrNucleusHE::operator= ( const G4ElasticHadrNucleusHE &  right)

◆ SampleInvariantT()

G4double G4ElasticHadrNucleusHE::SampleInvariantT ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)
virtual

Reimplemented from G4HadronElastic.

Definition at line 343 of file G4ElasticHadrNucleusHE.cc.

346{
347 G4double plab = inLabMom/GeV; // (GeV/c)
348 G4double Q2 = 0;
349
350 iHadrCode = p->GetPDGEncoding();
351
352 NumbN = N;
353
354 if(verboseLevel > 1)
355 {
356 G4cout<< " G4ElasticHadrNucleusHE::SampleT: "
357 << " for " << p->GetParticleName()
358 << " at Z= " << Z << " A= " << N
359 << " plab(GeV)= " << plab
360 << G4endl;
361 }
362 G4int idx;
363
364 for( idx = 0 ; idx < NHADRONS; idx++)
365 {
366 if(iHadrCode == HadronCode[idx]) break;
367 }
368
369 // Hadron is not in the list
370
371 if( idx >= NHADRONS ) return Q2;
372
373 iHadron = HadronType[idx];
374 iHadrCode = HadronCode[idx];
375
376 if(Z==1)
377 {
378 hMass = p->GetPDGMass()/GeV;
379 hMass2 = hMass*hMass;
380
381 G4double T = sqrt(plab*plab+hMass2)-hMass;
382
383 if(T > 0.4) Q2 = HadronProtonQ2(p, plab);
384
385 if (verboseLevel>1)
386 G4cout<<" Proton : Q2 "<<Q2<<G4endl;
387 }
388 else
389 {
390 G4ElasticData* ElD1 = SetOfElasticData[idx][Z];
391
392 // Construct elastic data
393 if(!ElD1)
394 {
395 G4double AWeight = nistManager->GetAtomicMassAmu(Z);
396 ElD1 = new G4ElasticData(p, Z, AWeight, Energy);
397 SetOfElasticData[idx][Z] = ElD1;
398
399 if(verboseLevel > 1)
400 {
401 G4cout<< " G4ElasticHadrNucleusHE::SampleT: new record " << idx
402 << " for " << p->GetParticleName() << " Z= " << Z
403 << G4endl;
404 }
405 }
406 hMass = ElD1->massGeV;
407 hMass2 = ElD1->mass2GeV2;
408 G4double M = ElD1->massA;
409 G4double M2 = ElD1->massA2;
410 G4double plab2 = plab*plab;
411 G4double Q2max = 4.*plab2*M2/
412 (hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2));
413
414 // sample scattering
415 G4double T = sqrt(plab2+hMass2)-hMass;
416
417 if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max);
418
419 if(verboseLevel > 1)
420 G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< " t/tmax= " << Q2/Q2max <<G4endl;
421 }
422 return Q2*GeV2;
423}
G4ElasticHadrNucleusHE::HadronNucleusQ2_2
G4double HadronNucleusQ2_2(G4ElasticData *pElD, G4int Z, G4double plabGeV, G4double tmax)
Definition: G4ElasticHadrNucleusHE.cc:442
G4ElasticHadrNucleusHE::HadronProtonQ2
G4double HadronProtonQ2(const G4ParticleDefinition *aHadron, G4double inLabMom)
Definition: G4ElasticHadrNucleusHE.cc:1417

Referenced by SampleT().

◆ SampleT()

G4double G4ElasticHadrNucleusHE::SampleT ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)

Definition at line 430 of file G4ElasticHadrNucleusHE.cc.

433{
434 return SampleInvariantT(p, inLabMom, Z, N);
435}
G4ElasticHadrNucleusHE::SampleInvariantT
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4ElasticHadrNucleusHE.cc:343
The documentation for this class was generated from the following files: