Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4ElasticHadrNucleusHE Class Reference

#include <G4ElasticHadrNucleusHE.hh>

Inheritance diagram for G4ElasticHadrNucleusHE:
G4HadronElastic G4HadronicInteraction

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)
 
void SetLowestEnergyLimit (G4double value)
 
G4double LowestEnergyLimit () const
 
G4double ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
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
 
virtual G4int GetVerboseLevel () const
 
virtual 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::G4ElasticHadrNucleusHE ( const G4String &  name = "hElasticGlauber")

Definition at line 226 of file G4ElasticHadrNucleusHE.cc.

References python.hepunit::GeV, G4NistManager::Instance(), python.hepunit::MeV, python.hepunit::proton_mass_c2, and G4HadronicInteraction::verboseLevel.

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 }
G4HadronElastic::G4HadronElastic
G4HadronElastic(const G4String &name="hElasticLHEP")
Definition: G4HadronElastic.cc:45
G4int
int G4int
Definition: G4Types.hh:78
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4double
double G4double
Definition: G4Types.hh:76
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 }
G4int
int G4int
Definition: G4Types.hh:78
G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE ( const G4ElasticHadrNucleusHE &  )

Member Function Documentation

void G4ElasticHadrNucleusHE::DefineHadronValues ( G4int  Z)

Definition at line 948 of file G4ElasticHadrNucleusHE.cc.

References G4cout, G4endl, InterpolateHN(), G4HadronicInteraction::verboseLevel, and test::x.

Referenced by GetKinematics(), and HadronNucleusQ2_2().

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 }
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
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
G4double
double G4double
Definition: G4Types.hh:76
void G4ElasticHadrNucleusHE::Description ( ) const
virtual

Reimplemented from G4HadronElastic.

Definition at line 296 of file G4ElasticHadrNucleusHE.cc.

References G4HadronicInteraction::GetModelName(), and outFile.

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 }
outFile
std::ofstream outFile
Definition: GammaRayTel.cc:68
G4HadronicInteraction::GetModelName
const G4String & GetModelName() const
Definition: G4HadronicInteraction.hh:119
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 }
G4double G4ElasticHadrNucleusHE::GetDistrFun ( G4double  Q2)
inline

Definition at line 278 of file G4ElasticHadrNucleusHE.hh.

References GetFt().

Referenced by GetQ2().

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

Definition at line 1369 of file G4ElasticHadrNucleusHE.cc.

References G4cout, G4endl, and G4HadronicInteraction::verboseLevel.

Referenced by GetDistrFun(), and GetQ2().

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 }
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4double G4ElasticHadrNucleusHE::GetHeavyFq2 ( G4int  Z,
G4int  Nucleus,
G4double *  LineFq2 
)

Definition at line 596 of file G4ElasticHadrNucleusHE.cc.

References G4cout, G4endl, HadrNucDifferCrSec(), and G4HadronicInteraction::verboseLevel.

Referenced by HadronNucleusQ2_2().

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 }
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
G4ElasticHadrNucleusHE::HadrNucDifferCrSec
G4double HadrNucDifferCrSec(G4int Z, G4int Nucleus, G4double Q2)
Definition: G4ElasticHadrNucleusHE.cc:759
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
void G4ElasticHadrNucleusHE::GetKinematics ( const G4ParticleDefinition *  aHadron,
G4double  MomentumH 
)

Definition at line 1329 of file G4ElasticHadrNucleusHE.cc.

References DefineHadronValues(), G4cout, G4endl, G4ParticleDefinition::GetParticleName(), and G4HadronicInteraction::verboseLevel.

Referenced by HadronProtonQ2().

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 }
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4double G4ElasticHadrNucleusHE::GetLightFq2 ( G4int  Z,
G4int  A,
G4double  Q 
)

Definition at line 646 of file G4ElasticHadrNucleusHE.cc.

References G4cout, G4endl, python.hepunit::pi, python.hepunit::twopi, and G4HadronicInteraction::verboseLevel.

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 }
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4double G4ElasticHadrNucleusHE::GetQ2 ( G4double  Ran)

Definition at line 1388 of file G4ElasticHadrNucleusHE.cc.

References GetDistrFun(), and GetFt().

Referenced by HadronProtonQ2().

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
G4double
double G4double
Definition: G4Types.hh:76
G4double G4ElasticHadrNucleusHE::GetQ2_2 ( G4int  N,
G4double *  Q,
G4double *  F,
G4double  R 
)

Definition at line 542 of file G4ElasticHadrNucleusHE.cc.

References F12, F22, F32, G4cout, G4endl, and G4HadronicInteraction::verboseLevel.

Referenced by HadronNucleusQ2_2().

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
G4cout
G4GLOB_DLL std::ostream G4cout
F32
#define F32
Definition: G4ErrorMatrix.cc:906
G4endl
#define G4endl
Definition: G4ios.hh:61
F12
#define F12
Definition: G4ErrorMatrix.cc:896
G4double
double G4double
Definition: G4Types.hh:76
G4double G4ElasticHadrNucleusHE::HadrNucDifferCrSec ( G4int  Z,
G4int  Nucleus,
G4double  Q2 
)

GeV/GeV;

Definition at line 759 of file G4ElasticHadrNucleusHE.cc.

References C1, C3, G4NistManager::GetAtomicMassAmu(), int(), N, python.hepunit::pi, and python.hepunit::twopi.

Referenced by GetHeavyFq2().

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
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4int
int G4int
Definition: G4Types.hh:78
C3
#define C3
C1
#define C1
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:344
N
**D E S C R I P T I O N
Definition: HepMCEx01/src/HEPEvtcom.cc:77
G4double
double G4double
Definition: G4Types.hh:76
G4double G4ElasticHadrNucleusHE::HadronNucleusQ2_2 ( G4ElasticData *  pElD,
G4int  Z,
G4double  plabGeV,
G4double  tmax 
)

Definition at line 442 of file G4ElasticHadrNucleusHE.cc.

References G4ElasticData::Aeff, DefineHadronValues(), G4ElasticData::dnkE, G4cout, G4endl, G4UniformRand, GetHeavyFq2(), GetQ2_2(), G4ElasticData::maxQ2, G4ElasticData::Pnucl, G4ElasticData::R1, G4ElasticData::R2, G4ElasticData::TableCrossSec, G4ElasticData::TableQ2, and G4HadronicInteraction::verboseLevel.

Referenced by SampleInvariantT().

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 }
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4ElasticData::TableCrossSec
G4double TableCrossSec[NQTABLE]
Definition: G4ElasticHadrNucleusHE.hh:101
G4ElasticData::maxQ2
G4double maxQ2[NENERGY]
Definition: G4ElasticHadrNucleusHE.hh:97
G4ElasticHadrNucleusHE::GetHeavyFq2
G4double GetHeavyFq2(G4int Z, G4int Nucleus, G4double *LineFq2)
Definition: G4ElasticHadrNucleusHE.cc:596
G4ElasticData::TableQ2
G4double TableQ2[ONQ2]
Definition: G4ElasticHadrNucleusHE.hh:100
G4endl
#define G4endl
Definition: G4ios.hh:61
G4ElasticHadrNucleusHE::GetQ2_2
G4double GetQ2_2(G4int N, G4double *Q, G4double *F, G4double R)
Definition: G4ElasticHadrNucleusHE.cc:542
G4double
double G4double
Definition: G4Types.hh:76
G4double G4ElasticHadrNucleusHE::HadronProtonQ2 ( const G4ParticleDefinition *  aHadron,
G4double  inLabMom 
)

Definition at line 1417 of file G4ElasticHadrNucleusHE.cc.

References G4UniformRand, GetKinematics(), G4ParticleDefinition::GetPDGMass(), GetQ2(), and python.hepunit::GeV.

Referenced by SampleInvariantT().

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
p
const char * p
Definition: xmltok.h:285
G4ElasticHadrNucleusHE::GetKinematics
void GetKinematics(const G4ParticleDefinition *aHadron, G4double MomentumH)
Definition: G4ElasticHadrNucleusHE.cc:1329
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4double
double G4double
Definition: G4Types.hh:76
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.

References LineInterpol(), and n.

Referenced by DefineHadronValues().

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
G4int
int G4int
Definition: G4Types.hh:78
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4double G4ElasticHadrNucleusHE::LineInterpol ( G4double  p0,
G4double  p2,
G4double  c1,
G4double  c2,
G4double  p 
)
inline

Definition at line 233 of file G4ElasticHadrNucleusHE.hh.

Referenced by InterpolateHN().

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 }
p
const char * p
Definition: xmltok.h:285
plottest35.c1
tuple c1
Definition: plottest35.py:14
G4ElasticHadrNucleusHE& G4ElasticHadrNucleusHE::operator= ( const G4ElasticHadrNucleusHE &  right)
G4double G4ElasticHadrNucleusHE::SampleInvariantT ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)
virtual

Reimplemented from G4HadronElastic.

Definition at line 343 of file G4ElasticHadrNucleusHE.cc.

References G4cout, G4endl, G4NistManager::GetAtomicMassAmu(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), python.hepunit::GeV, HadronNucleusQ2_2(), HadronProtonQ2(), G4ElasticData::mass2GeV2, G4ElasticData::massA, G4ElasticData::massA2, G4ElasticData::massGeV, N, and G4HadronicInteraction::verboseLevel.

Referenced by SampleT().

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::HadronProtonQ2
G4double HadronProtonQ2(const G4ParticleDefinition *aHadron, G4double inLabMom)
Definition: G4ElasticHadrNucleusHE.cc:1417
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4cout
G4GLOB_DLL std::ostream G4cout
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:344
G4endl
#define G4endl
Definition: G4ios.hh:61
N
**D E S C R I P T I O N
Definition: HepMCEx01/src/HEPEvtcom.cc:77
G4ElasticHadrNucleusHE::HadronNucleusQ2_2
G4double HadronNucleusQ2_2(G4ElasticData *pElD, G4int Z, G4double plabGeV, G4double tmax)
Definition: G4ElasticHadrNucleusHE.cc:442
G4double
double G4double
Definition: G4Types.hh:76
G4double G4ElasticHadrNucleusHE::SampleT ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)

Definition at line 430 of file G4ElasticHadrNucleusHE.cc.

References SampleInvariantT().

433 {
434  return SampleInvariantT(p, inLabMom, Z, N);
435 }
N
**D E S C R I P T I O N
Definition: HepMCEx01/src/HEPEvtcom.cc:77
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:
Generated on Wed Apr 30 2014 15:56:21 for Geant4.10 by   doxygen 1.8.7