Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
G4ElasticHadrNucleusHE Class Reference

#include <G4ElasticHadrNucleusHE.hh>

Inheritance diagram for G4ElasticHadrNucleusHE:
Inheritance graph
[legend]
Collaboration diagram for G4ElasticHadrNucleusHE:
Collaboration graph
[legend]

Public Member Functions

 G4ElasticHadrNucleusHE (const G4String &name="hElasticGlauber")
 
virtual ~G4ElasticHadrNucleusHE ()
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual void ModelDescription (std::ostream &) 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 &aTrack, G4Nucleus &targetNucleus)
 
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)
 
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
 
virtual const std::pair
< G4double, G4double > 		GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double,
G4double > 		GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void InitialiseModel ()
 

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 239 of file G4ElasticHadrNucleusHE.cc.

240  : G4HadronElastic(name)
241 {
242  //ANDREA->
243  G4AutoLock l(&aMutex);
244  if ( onlyOnceInit ) {
245  for ( int i = 0 ; i< NHADRONS ; ++i) {
246  for (int j = 0 ; j<ZMAX ; ++j ) {
247  SetOfElasticData[i][j]=0;
248  G4MUTEXINIT(eldata_m[i][j]);
249  }
250  }
251  onlyOnceInit = false;
252  }
253  l.unlock();
254  //ANDREA<-
255 
256  dQ2 = hMass = hMass2 = hLabMomentum = hLabMomentum2 = MomentumCM = HadrEnergy
257  = R1 = R2 = Pnucl = Aeff = HadrTot = HadrSlope = HadrReIm = TotP = DDSect2
258  = DDSect3 = ConstU = FmaxT = Slope1 = Slope2 = Coeff1 = Coeff2 = MaxTR
259  = Slope0 = Coeff0 = aAIm = aDIm = Dtot11 = 0.0;
260  NumbN = iHadrCode = iHadron = 0;
261 
262  verboseLevel = 0;
263  plabLowLimit = 20.0*MeV;
264  lowestEnergyLimit = 0.0;
265  //Description();
266 
267  MbToGeV2 = 2.568;
268  sqMbToGeV = 1.602;
269  Fm2ToGeV2 = 25.68;
270  GeV2 = GeV*GeV;
271  protonM = proton_mass_c2/GeV;
272  protonM2 = protonM*protonM;
273 
274  BoundaryP[0]=9.0;BoundaryTG[0]=5.0;BoundaryTL[0]=0.;
275  BoundaryP[1]=20.0;BoundaryTG[1]=1.5;BoundaryTL[1]=0.;
276  BoundaryP[2]=5.0; BoundaryTG[2]=1.0;BoundaryTL[2]=1.5;
277  BoundaryP[3]=8.0; BoundaryTG[3]=3.0;BoundaryTL[3]=0.;
278  BoundaryP[4]=7.0; BoundaryTG[4]=3.0;BoundaryTL[4]=0.;
279  BoundaryP[5]=5.0; BoundaryTG[5]=2.0;BoundaryTL[5]=0.;
280  BoundaryP[6]=5.0; BoundaryTG[6]=1.5;BoundaryTL[6]=3.0;
281 
282  Binom();
283  // energy in GeV
284  Energy[0] = 0.4;
285  Energy[1] = 0.6;
286  Energy[2] = 0.8;
287  LowEdgeEnergy[0] = 0.0;
288  LowEdgeEnergy[1] = 0.5;
289  LowEdgeEnergy[2] = 0.7;
290  G4double e = 1.0;
291  G4double f = G4Exp(G4Log(10.)*0.1);
292  for(G4int i=3; i<NENERGY; i++) {
293  Energy[i] = e;
294  LowEdgeEnergy[i] = e/f;
295  e *= f*f;
296  }
297  nistManager = G4NistManager::Instance();
298 
299  // PDG code for hadrons
300  G4int ic[NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311,
301  3122,3222,3112,3212,3312,3322,3334,
302  -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334};
303  // internal index
304  G4int id[NHADRONS] = {2,3,6,0,4,5,4,4,4,5,
305  0,0,0,0,0,0,0,
306  1,7,1,1,1,1,1,1,1};
307 
308  G4int id1[NHADRONS] = {3,4,1,0,5,6,5,5,5,6,
309  0,0,0,0,0,0,0,
310  2,2,2,2,2,2,2,2,2};
311 
312  for(G4int j=0; j<NHADRONS; j++)
313  {
314  HadronCode[j] = ic[j];
315  HadronType[j] = id[j];
316  HadronType1[j] = id1[j];
317 
318  for(G4int k = 0; k < ZMAX; k++) { SetOfElasticData[j][k] = 0; }
319  }
320 }
G4HadronElastic::G4HadronElastic
G4HadronElastic(const G4String &name="hElasticLHEP")
Definition: G4HadronElastic.cc:48
G4MUTEXINIT
#define G4MUTEXINIT(mutex)
Definition: G4Threading.hh:177
G4int
int G4int
Definition: G4Types.hh:78
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
ZMAX
static const G4int ZMAX
Definition: G4ElasticHadrNucleusHE.hh:62
aMutex
G4Mutex aMutex
Definition: G4DNAMolecularMaterial.cc:42
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
NENERGY
static const G4int NENERGY
Definition: G4ElasticHadrNucleusHE.hh:61
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
NHADRONS
static const G4int NHADRONS
Definition: G4ElasticHadrNucleusHE.hh:58

Here is the call graph for this function:

G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusHE ( )
virtual

Definition at line 338 of file G4ElasticHadrNucleusHE.cc.

339 {
340  //ANDREA->
341  G4AutoLock l(&aMutex);
342  if ( onlyOnceDestroy ) {
343  for(G4int j = 0; j < NHADRONS; j++)
344  {
345  for(G4int k = 0; k < ZMAX; k++)
346  {
347  if ( SetOfElasticData[j][k] ) {
348  delete SetOfElasticData[j][k];
349  SetOfElasticData[j][k]=0;
350  G4MUTEXDESTROY(eldata_m[j][k]);
351  }
352  }
353  }
354  onlyOnceDestroy = false;
355  }
356 }
G4int
int G4int
Definition: G4Types.hh:78
ZMAX
static const G4int ZMAX
Definition: G4ElasticHadrNucleusHE.hh:62
aMutex
G4Mutex aMutex
Definition: G4DNAMolecularMaterial.cc:42
G4MUTEXDESTROY
#define G4MUTEXDESTROY(mutex)
Definition: G4Threading.hh:178
NHADRONS
static const G4int NHADRONS
Definition: G4ElasticHadrNucleusHE.hh:58
G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE ( const G4ElasticHadrNucleusHE &  )

Member Function Documentation

void G4ElasticHadrNucleusHE::DefineHadronValues ( G4int  Z)

Definition at line 975 of file G4ElasticHadrNucleusHE.cc.

976 {
977  HadrEnergy = std::sqrt(hMass2 + hLabMomentum2);
978 
979  G4double sHadr = 2.*HadrEnergy*protonM+protonM2+hMass2;
980  G4double sqrS = std::sqrt(sHadr);
981  G4double Ecm = 0.5*(sHadr-hMass2+protonM2)/sqrS;
982  MomentumCM = std::sqrt(Ecm*Ecm-protonM2);
983 
984  if(verboseLevel>2)
985  G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron
986  << " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS
987  << " plab= " << hLabMomentum
988  <<" E - m "<<HadrEnergy - hMass<< G4endl;
989 
990  G4double TotN = 0.0;
991  G4double logE = G4Log(HadrEnergy);
992  G4double logS = G4Log(sHadr);
993  TotP = 0.0;
994 
995  switch (iHadron)
996  {
997  case 0: // proton, neutron
998  case 6:
999 
1000  if(hLabMomentum > 10)
1001  TotP = TotN = 7.5*logE - 40.12525 + 103*G4Exp(-G4Log(sHadr)*0.165);// mb
1002 
1003  else
1004  {
1005 // ================== neutron ================
1006 
1008 
1009 
1010  if( hLabMomentum > 1.4 )
1011  TotN = 33.3+15.2*(hLabMomentum2-1.35)/
1012  (G4Exp(G4Log(hLabMomentum)*2.37)+0.95);
1013 
1014  else if(hLabMomentum > 0.8)
1015  {
1016  G4double A0 = logE + 0.0513;
1017  TotN = 33.0 + 25.5*A0*A0;
1018  }
1019  else
1020  {
1021  G4double A0 = logE - 0.2634; // log(1.3)
1022  TotN = 33.0 + 30.*A0*A0*A0*A0;
1023  }
1024 // ================= proton ===============
1025 // else if(iHadrCode == 2212)
1026  {
1027  if(hLabMomentum >= 1.05)
1028  {
1029  TotP = 39.0+75.*(hLabMomentum-1.2)/
1030  (hLabMomentum2*hLabMomentum+0.15);
1031  }
1032 
1033  else if(hLabMomentum >= 0.7)
1034  {
1035  G4double A0 = logE + 0.3147;
1036  TotP = 23.0 + 40.*A0*A0;
1037  }
1038  else
1039  {
1040  TotP = 23.+50.*G4Exp(G4Log(G4Log(0.73/hLabMomentum))*3.5);
1041  }
1042  }
1043  }
1044 
1045 // HadrTot = 0.5*(82*TotP+126*TotN)/104; // $$$$$$$$$$$$$$$$$$
1046  HadrTot = 0.5*(TotP+TotN);
1047 // ...................................................
1048  // Proton slope
1049  if(hLabMomentum >= 2.) HadrSlope = 5.44 + 0.88*logS;
1050 
1051  else if(hLabMomentum >= 0.5) HadrSlope = 3.73*hLabMomentum-0.37;
1052 
1053  else HadrSlope = 1.5;
1054 
1055 // ...................................................
1056  if(hLabMomentum >= 1.2)
1057  HadrReIm = 0.13*(logS - 5.8579332)*G4Exp(-G4Log(sHadr)*0.18);
1058 
1059  else if(hLabMomentum >= 0.6)
1060  HadrReIm = -75.5*(G4Exp(G4Log(hLabMomentum)*0.25)-0.95)/
1061  (G4Exp(G4Log(3*hLabMomentum)*2.2)+1);
1062 
1063  else
1064  HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2);
1065 // ...................................................
1066  DDSect2 = 2.2; //mb*GeV-2
1067  DDSect3 = 0.6; //mb*GeV-2
1068  // ================== lambda ==================
1069  if( iHadrCode == 3122)
1070  {
1071  HadrTot *= 0.88;
1072  HadrSlope *=0.85;
1073  }
1074  // ================== sigma + ==================
1075  else if( iHadrCode == 3222)
1076  {
1077  HadrTot *=0.81;
1078  HadrSlope *=0.85;
1079  }
1080  // ================== sigma 0,- ==================
1081  else if(iHadrCode == 3112 || iHadrCode == 3212 )
1082  {
1083  HadrTot *=0.88;
1084  HadrSlope *=0.85;
1085  }
1086  // =================== xi =================
1087  else if( iHadrCode == 3312 || iHadrCode == 3322 )
1088  {
1089  HadrTot *=0.77;
1090  HadrSlope *=0.75;
1091  }
1092  // ================= omega =================
1093  else if( iHadrCode == 3334)
1094  {
1095  HadrTot *=0.78;
1096  HadrSlope *=0.7;
1097  }
1098 
1099  break;
1100 // ===========================================================
1101  case 1: // antiproton
1102  case 7: // antineutron
1103 
1104  HadrTot = 5.2+5.2*logE + 123.2/sqrS; // mb
1105  HadrSlope = 8.32+0.57*logS; //(GeV/c)^-2
1106 
1107  if( HadrEnergy < 1000 )
1108  HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*G4Exp(-G4Log(sHadr)*0.8);
1109  else
1110  HadrReIm = 0.6*(logS - 5.8579332)*G4Exp(-G4Log(sHadr)*0.25);
1111 
1112  DDSect2 = 11; //mb*(GeV/c)^-2
1113  DDSect3 = 3; //mb*(GeV/c)^-2
1114  // ================== lambda ==================
1115  if( iHadrCode == -3122)
1116  {
1117  HadrTot *= 0.88;
1118  HadrSlope *=0.85;
1119  }
1120  // ================== sigma + ==================
1121  else if( iHadrCode == -3222)
1122  {
1123  HadrTot *=0.81;
1124  HadrSlope *=0.85;
1125  }
1126  // ================== sigma 0,- ==================
1127  else if(iHadrCode == -3112 || iHadrCode == -3212 )
1128  {
1129  HadrTot *=0.88;
1130  HadrSlope *=0.85;
1131  }
1132  // =================== xi =================
1133  else if( iHadrCode == -3312 || iHadrCode == -3322 )
1134  {
1135  HadrTot *=0.77;
1136  HadrSlope *=0.75;
1137  }
1138  // ================= omega =================
1139  else if( iHadrCode == -3334)
1140  {
1141  HadrTot *=0.78;
1142  HadrSlope *=0.7;
1143  }
1144 
1145  break;
1146 // -------------------------------------------
1147  case 2: // pi plus, pi minus
1148  case 3:
1149 
1150  if(hLabMomentum >= 3.5)
1151  TotP = 10.6+2.*logE + 25.*G4Exp(-G4Log(HadrEnergy)*0.43); // mb
1152 // =========================================
1153  else if(hLabMomentum >= 1.15)
1154  {
1155  G4double x = (hLabMomentum - 2.55)/0.55;
1156  G4double y = (hLabMomentum - 1.47)/0.225;
1157  TotP = 3.2*G4Exp(-x*x) + 12.*G4Exp(-y*y) + 27.5;
1158  }
1159 // =========================================
1160  else if(hLabMomentum >= 0.4)
1161  {
1162  TotP = 88*(logE+0.2877)*(logE+0.2877)+14.0;
1163  }
1164 // =========================================
1165  else
1166  {
1167  G4double x = (hLabMomentum - 0.29)/0.085;
1168  TotP = 20. + 180.*G4Exp(-x*x);
1169  }
1170 // -------------------------------------------
1171 
1172  if(hLabMomentum >= 3.0 )
1173  TotN = 10.6 + 2.*logE + 30.*G4Exp(-G4Log(HadrEnergy)*0.43); // mb
1174 
1175  else if(hLabMomentum >= 1.3)
1176  {
1177  G4double x = (hLabMomentum - 2.1)/0.4;
1178  G4double y = (hLabMomentum - 1.4)/0.12;
1179  TotN = 36.1+0.079 - 4.313*logE + 3.*G4Exp(-x*x) +
1180  1.5*G4Exp(-y*y);
1181  }
1182  else if(hLabMomentum >= 0.65)
1183  {
1184  G4double x = (hLabMomentum - 0.72)/0.06;
1185  G4double y = (hLabMomentum - 1.015)/0.075;
1186  TotN = 36.1 + 10.*G4Exp(-x*x) + 24*G4Exp(-y*y);
1187  }
1188  else if(hLabMomentum >= 0.37)
1189  {
1190  G4double x = G4Log(hLabMomentum/0.48);
1191  TotN = 26. + 110.*x*x;
1192  }
1193  else
1194  {
1195  G4double x = (hLabMomentum - 0.29)/0.07;
1196  TotN = 28.0 + 40.*G4Exp(-x*x);
1197  }
1198  HadrTot = (TotP+TotN)/2;
1199 // ........................................
1200  HadrSlope = 7.28+0.245*logS; // GeV-2
1201  HadrReIm = 0.2*(logS - 4.6051702)*G4Exp(-G4Log(sHadr)*0.15);
1202 
1203  DDSect2 = 0.7; //mb*GeV-2
1204  DDSect3 = 0.27; //mb*GeV-2
1205 
1206  break;
1207 // ==========================================================
1208  case 4: // K plus
1209 
1210  HadrTot = 10.6+1.8*logE + 9.0*G4Exp(-G4Log(HadrEnergy)*0.55); // mb
1211  if(HadrEnergy>100) HadrSlope = 15.0;
1212  else HadrSlope = 1.0+1.76*logS - 2.84/sqrS; // GeV-2
1213 
1214  HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*G4Exp(-G4Log(sHadr+50)*2.1);
1215  DDSect2 = 0.7; //mb*GeV-2
1216  DDSect3 = 0.21; //mb*GeV-2
1217  break;
1218 // =========================================================
1219  case 5: // K minus
1220 
1221  HadrTot = 10+1.8*logE + 25./sqrS; // mb
1222  HadrSlope = 6.98+0.127*logS; // GeV-2
1223  HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*G4Exp(-G4Log(sHadr+50)*2.1);
1224  DDSect2 = 0.7; //mb*GeV-2
1225  DDSect3 = 0.27; //mb*GeV-2
1226  break;
1227  }
1228 // =========================================================
1229  if(verboseLevel>2)
1230  G4cout << "HadrTot= " << HadrTot << " HadrSlope= " << HadrSlope
1231  << " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2
1232  << " DDSect3= " << DDSect3 << G4endl;
1233 
1234  if(Z != 1) return;
1235 
1236  // Scattering of protons
1237 
1238  Coeff0 = Coeff1 = Coeff2 = 0.0;
1239  Slope0 = Slope1 = 1.0;
1240  Slope2 = 5.0;
1241 
1242  // data for iHadron=0
1243  static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0};
1244  static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002};
1245  static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0};
1246  static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25};
1247  static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8};
1248 
1249  // data for iHadron=6,7
1250  static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0};
1251  static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01};
1252  static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003};
1253  static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5};
1254  static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8};
1255 
1256  // data for iHadron=1
1257  static const G4double EnP[2]={1.5,4.0};
1258  static const G4double C0P[2]={0.001,0.0005};
1259  static const G4double C1P[2]={0.003,0.001};
1260  static const G4double B0P[2]={2.5,4.5};
1261  static const G4double B1P[2]={1.0,4.0};
1262 
1263  // data for iHadron=2
1264  static const G4double EnPP[4]={1.0,2.0,3.0,4.0};
1265  static const G4double C0PP[4]={0.0,0.0,0.0,0.0};
1266  static const G4double C1PP[4]={0.15,0.08,0.02,0.01};
1267  static const G4double B0PP[4]={1.5,2.8,3.8,3.8};
1268  static const G4double B1PP[4]={0.8,1.6,3.6,4.6};
1269 
1270  // data for iHadron=3
1271  static const G4double EnPPN[4]={1.0,2.0,3.0,4.0};
1272  static const G4double C0PPN[4]={0.0,0.0,0.0,0.0};
1273  static const G4double C1PPN[4]={0.0,0.0,0.0,0.0};
1274  static const G4double B0PPN[4]={1.5,2.8,3.8,3.8};
1275  static const G4double B1PPN[4]={0.8,1.6,3.6,4.6};
1276 
1277  // data for iHadron=4
1278  static const G4double EnK[4]={1.4,2.33,3.0,5.0};
1279  static const G4double C0K[4]={0.0,0.0,0.0,0.0};
1280  static const G4double C1K[4]={0.01,0.007,0.005,0.003};
1281  static const G4double B0K[4]={1.5,2.0,3.8,3.8};
1282  static const G4double B1K[4]={1.6,1.6,1.6,1.6};
1283 
1284  // data for iHadron=5
1285  static const G4double EnKM[2]={1.4,4.0};
1286  static const G4double C0KM[2]={0.006,0.002};
1287  static const G4double C1KM[2]={0.00,0.00};
1288  static const G4double B0KM[2]={2.5,3.5};
1289  static const G4double B1KM[2]={1.6,1.6};
1290 
1291  switch(iHadron)
1292  {
1293  case 0 :
1294 
1295  if(hLabMomentum <BoundaryP[0])
1296  InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0);
1297 
1298  Coeff2 = 0.8/hLabMomentum/hLabMomentum;
1299  break;
1300 
1301  case 6 :
1302 
1303  if(hLabMomentum < BoundaryP[1])
1304  InterpolateHN(5,EnN,C0N,C1N,B0N,B1N);
1305 
1306  Coeff2 = 0.8/hLabMomentum/hLabMomentum;
1307  break;
1308 
1309  case 1 :
1310  case 7 :
1311  if(hLabMomentum < BoundaryP[2])
1312  InterpolateHN(2,EnP,C0P,C1P,B0P,B1P);
1313  break;
1314 
1315  case 2 :
1316 
1317  if(hLabMomentum < BoundaryP[3])
1318  InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1PP);
1319 
1320  Coeff2 = 0.02/hLabMomentum;
1321  break;
1322 
1323  case 3 :
1324 
1325  if(hLabMomentum < BoundaryP[4])
1326  InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN,B1PPN);
1327 
1328  Coeff2 = 0.02/hLabMomentum;
1329  break;
1330 
1331  case 4 :
1332 
1333  if(hLabMomentum < BoundaryP[5])
1334  InterpolateHN(4,EnK,C0K,C1K,B0K,B1K);
1335 
1336  if(hLabMomentum < 1) Coeff2 = 0.34;
1337  else Coeff2 = 0.34/hLabMomentum2/hLabMomentum;
1338  break;
1339 
1340  case 5 :
1341  if(hLabMomentum < BoundaryP[6])
1342  InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1KM);
1343 
1344  if(hLabMomentum < 1) Coeff2 = 0.01;
1345  else Coeff2 = 0.01/hLabMomentum2/hLabMomentum;
1346  break;
1347  }
1348 
1349  if(verboseLevel > 2)
1350  G4cout<<" HadrVal : Plasb "<<hLabMomentum
1351  <<" iHadron "<<iHadron<<" HadrTot "<<HadrTot<<G4endl;
1352 }
test.x
tuple x
Definition: test.py:50
G4cout
G4GLOB_DLL std::ostream G4cout
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
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:248
G4double
double G4double
Definition: G4Types.hh:76
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 270 of file G4ElasticHadrNucleusHE.hh.

271 {
272  if ( numN >= numM && numN <= 240) return SetBinom[numN][numM];
273  else return 0.;
274 }
G4double G4ElasticHadrNucleusHE::GetDistrFun ( G4double  Q2)
inline

Definition at line 279 of file G4ElasticHadrNucleusHE.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ElasticHadrNucleusHE::GetFt ( G4double  Q2)

Definition at line 1396 of file G4ElasticHadrNucleusHE.cc.

1397 {
1398  G4double Fdistr=0;
1399  G4double SqrQ2 = std::sqrt(Q2);
1400 
1401  Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*G4Exp(ConstU))
1402  /HadrSlope*(1-G4Exp(-HadrSlope*Q2))
1403  + Coeff0*(1-G4Exp(-Slope0*Q2))
1404  + Coeff2/Slope2*G4Exp(Slope2*ConstU)*(G4Exp(Slope2*Q2)-1)
1405  + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*G4Exp(-Slope1*SqrQ2));
1406 
1407  if (verboseLevel>1)
1408  G4cout<<"Old: Coeff0 Coeff1 Coeff2 "<<Coeff0<<" "
1409  <<Coeff1<<" "<<Coeff2<<" Slope Slope0 Slope1 Slope2 "
1410  <<HadrSlope<<" "<<Slope0<<" "<<Slope1<<" "<<Slope2
1411  <<" Fdistr "<<Fdistr<<G4endl;
1412  return Fdistr;
1413 }
G4cout
G4GLOB_DLL std::ostream G4cout
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 625 of file G4ElasticHadrNucleusHE.cc.

626 {
627  G4int ii, jj, aSimp;
628  G4double curQ2, curSec;
629  G4double curSum = 0.0;
630  G4double totSum = 0.0;
631 
632  G4double ddQ2 = dQ2/20;
633  G4double Q2l = 0;
634 
635  LineF[0] = 0;
636  for(ii = 1; ii<ONQ2; ii++)
637  {
638  curSum = 0;
639  aSimp = 4;
640 
641  for(jj = 0; jj<20; jj++)
642  {
643  curQ2 = Q2l+jj*ddQ2;
644 
645  curSec = HadrNucDifferCrSec(Z, Nucleus, curQ2);
646  curSum += curSec*aSimp;
647 
648  if(aSimp > 3) aSimp = 2;
649  else aSimp = 4;
650 
651  if(jj == 0 && verboseLevel>2)
652  G4cout<<" Q2 "<<curQ2<<" AIm "<<aAIm<<" DIm "<<aDIm
653  <<" Diff "<<curSec<<" totSum "<<totSum<<G4endl;
654  }
655 
656  Q2l += dQ2;
657  curSum *= ddQ2/2.3; // $$$$$$$$$$$$$$$$$$$$$$$
658  totSum += curSum;
659 
660  LineF[ii] = totSum;
661 
662  if (verboseLevel>2)
663  G4cout<<" GetHeavy: Q2 dQ2 totSum "<<Q2l<<" "<<dQ2<<" "<<totSum
664  <<" curSec "
665  <<curSec<<" totSum "<< totSum<<" DTot "
666  <<curSum<<G4endl;
667  }
668  return totSum;
669 }
ONQ2
static const G4int ONQ2
Definition: G4ElasticHadrNucleusHE.hh:60
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:788
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 1356 of file G4ElasticHadrNucleusHE.cc.

1358 {
1359  if (verboseLevel>1)
1360  G4cout<<"1 GetKin.: HadronName MomentumH "
1361  <<aHadron->GetParticleName()<<" "<<MomentumH<<G4endl;
1362 
1363  DefineHadronValues(1);
1364 
1365  G4double Sh = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV
1366 
1367  ConstU = 2*protonM2+2*hMass2-Sh;
1368 
1369  G4double MaxT = 4*MomentumCM*MomentumCM;
1370 
1371  BoundaryTL[0] = MaxT; //2.0;
1372  BoundaryTL[1] = MaxT;
1373  BoundaryTL[3] = MaxT;
1374  BoundaryTL[4] = MaxT;
1375  BoundaryTL[5] = MaxT;
1376 
1377  G4int NumberH=0;
1378 
1379  while(iHadrCode!=HadronCode[NumberH]) NumberH++; /* Loop checking, 10.08.2015, A.Ribon */
1380 
1381  NumberH = HadronType1[NumberH];
1382 
1383  if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH];
1384  else MaxTR = BoundaryTG[NumberH];
1385 
1386  if (verboseLevel>1)
1387  G4cout<<"3 GetKin. : NumberH "<<NumberH
1388  <<" Bound.P[NumberH] "<<BoundaryP[NumberH]
1389  <<" Bound.TL[NumberH] "<<BoundaryTL[NumberH]
1390  <<" Bound.TG[NumberH] "<<BoundaryTG[NumberH]
1391  <<" MaxT MaxTR "<<MaxT<<" "<<MaxTR<<G4endl;
1392 
1393 // GetParametersHP(aHadron, MomentumH);
1394 }
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 675 of file G4ElasticHadrNucleusHE.cc.

677 {
678  // Scattering of proton
679  if(Z == 1)
680  {
681  G4double SqrQ2 = std::sqrt(Q2);
682  G4double valueConstU = 2.*(hMass2 + protonM2) - Q2;
683 
684  G4double y = (1.-Coeff1-Coeff0)/HadrSlope*(1.-G4Exp(-HadrSlope*Q2))
685  + Coeff0*(1.-G4Exp(-Slope0*Q2))
686  + Coeff2/Slope2*G4Exp(Slope2*valueConstU)*(G4Exp(Slope2*Q2)-1.)
687  + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1+SqrQ2)*G4Exp(-Slope1*SqrQ2));
688 
689  return y;
690  }
691 
692  // The preparing of probability function
693 
694  G4double prec = Nucleus > 208 ? 1.0e-7 : 1.0e-6;
695 
696  G4double Stot = HadrTot*MbToGeV2; // Gev^-2
697  G4double Bhad = HadrSlope; // GeV^-2
698  G4double Asq = 1+HadrReIm*HadrReIm;
699  G4double Rho2 = std::sqrt(Asq);
700 
701 // if(verboseLevel >1)
702  G4cout<<" Fq2 Before for i Tot B Im "<<HadrTot<<" "<<HadrSlope<<" "
703  <<HadrReIm<<G4endl;
704 
705  if(verboseLevel > 1) {
706  G4cout << "GetFq2: Stot= " << Stot << " Bhad= " << Bhad
707  <<" Im "<<HadrReIm
708  << " Asq= " << Asq << G4endl;
709  G4cout << "R1= " << R1 << " R2= " << R2 << " Pnucl= " << Pnucl <<G4endl;
710  }
711  G4double R12 = R1*R1;
712  G4double R22 = R2*R2;
713  G4double R12B = R12+2*Bhad;
714  G4double R22B = R22+2*Bhad;
715 
716  G4double Norm = (R12*R1-Pnucl*R22*R2); // HP->Aeff;
717 
718  G4double R13 = R12*R1/R12B;
719  G4double R23 = Pnucl*R22*R2/R22B;
720  G4double Unucl = Stot/twopi/Norm*R13;
721  G4double UnucRho2 = -Unucl*Rho2;
722 
723  G4double FiH = std::asin(HadrReIm/Rho2);
724  G4double NN2 = R23/R13;
725 
726  if(verboseLevel > 2)
727  G4cout << "UnucRho2= " << UnucRho2 << " FiH= " << FiH << " NN2= " << NN2
728  << " Norm= " << Norm << G4endl;
729 
730  G4double dddd;
731 
732  G4double Prod0 = 0;
733  G4double N1 = -1.0;
734  //G4double Tot0 = 0;
735  G4double exp1;
736 
737  G4double Prod3 ;
738  G4double exp2 ;
739  G4double N4, N5, N2, Prod1, Prod2;
740  G4int i1, i2, j1, j2;
741 
742  for(i1 = 1; i1<= Nucleus; i1++)
743  {
744  N1 = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2;
745  Prod1 = 0;
746  //Tot0 = 0;
747  N2 = -1;
748 
749  for(i2 = 1; i2<=Nucleus; i2++)
750  {
751  N2 = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2;
752  Prod2 = 0;
753  N5 = -1/NN2;
754  for(j2=0; j2<= i2; j2++)
755  {
756  Prod3 = 0;
757  exp2 = 1/(j2/R22B+(i2-j2)/R12B);
758  N5 = -N5*NN2;
759  N4 = -1/NN2;
760  for(j1=0; j1<=i1; j1++)
761  {
762  exp1 = 1/(j1/R22B+(i1-j1)/R12B);
763  dddd = exp1+exp2;
764  N4 = -N4*NN2;
765  Prod3 = Prod3+N4*exp1*exp2*
766  (1-G4Exp(-Q2*dddd/4))/dddd*4*SetBinom[i1][j1];
767  } // j1
768  Prod2 = Prod2 +Prod3*N5*SetBinom[i2][j2];
769  } // j2
770  Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2));
771 
772  if (std::fabs(Prod2*N2/Prod1)<prec) break;
773  } // i2
774  Prod0 = Prod0 + Prod1*N1;
775  if(std::fabs(N1*Prod1/Prod0) < prec) break;
776 
777  } // i1
778 
779  Prod0 *= 0.25*pi/MbToGeV2; // This is in mb
780 
781  if(verboseLevel>1)
782  G4cout << "GetLightFq2 Z= " << Z << " A= " << Nucleus
783  <<" Q2= " << Q2 << " Res= " << Prod0 << G4endl;
784  return Prod0;
785 }
G4int
int G4int
Definition: G4Types.hh:78
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
CLHEP::prec
static const double prec
Definition: RanecuEngine.cc:58
G4cout
G4GLOB_DLL std::ostream G4cout
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4endl
#define G4endl
Definition: G4ios.hh:61
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4double
double G4double
Definition: G4Types.hh:76
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

G4double G4ElasticHadrNucleusHE::GetQ2 ( G4double  Ran)

Definition at line 1415 of file G4ElasticHadrNucleusHE.cc.

1416 {
1417  G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta;
1418  G4double Q2=0;
1419 
1420  FmaxT = GetFt(MaxTR);
1421  delta = GetDistrFun(DDD0)-Ran;
1422 
1423  const G4int maxNumberOfLoops = 10000;
1424  G4int loopCounter = -1;
1425  while ( (std::fabs(delta) > 0.0001) &&
1426  ++loopCounter < maxNumberOfLoops ) /* Loop checking, 10.08.2015, A.Ribon */
1427  {
1428  if(delta>0)
1429  {
1430  DDD2 = DDD0;
1431  DDD0 = (DDD0+DDD1)*0.5;
1432  }
1433  else if(delta<0)
1434  {
1435  DDD1 = DDD0;
1436  DDD0 = (DDD0+DDD2)*0.5;
1437  }
1438  delta = GetDistrFun(DDD0)-Ran;
1439  }
1440  if ( loopCounter >= maxNumberOfLoops ) {
1441  return 0.0;
1442  }
1443 
1444  Q2 = DDD0;
1445 
1446  return Q2;
1447 }
G4ElasticHadrNucleusHE::GetDistrFun
G4double GetDistrFun(G4double Q2)
Definition: G4ElasticHadrNucleusHE.hh:279
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 571 of file G4ElasticHadrNucleusHE.cc.

573 {
574  G4double ranQ2;
575 
576  G4double F1 = F[kk-2];
577  G4double F2 = F[kk-1];
578  G4double F3 = F[kk];
579  G4double X1 = Q[kk-2];
580  G4double X2 = Q[kk-1];
581  G4double X3 = Q[kk];
582 
583  if(verboseLevel > 2)
584  G4cout << "GetQ2_2 kk= " << kk << " X2= " << X2 << " X3= " << X3
585  << " F2= " << F2 << " F3= " << F3 << " Rndm= " << ranUni << G4endl;
586 
587  if(kk == 1 || kk == 0)
588  {
589  F1 = F[0];
590  F2 = F[1];
591  F3 = F[2];
592  X1 = Q[0];
593  X2 = Q[1];
594  X3 = Q[2];
595  }
596 
597  G4double F12 = F1*F1;
598  G4double F22 = F2*F2;
599  G4double F32 = F3*F3;
600 
601  G4double D0 = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3;
602 
603  if(verboseLevel > 2)
604  G4cout << " X1= " << X1 << " F1= " << F1 << " D0= "
605  << D0 << G4endl;
606 
607  if(std::abs(D0) < 0.00000001)
608  {
609  ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
610  }
611  else
612  {
613  G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1;
614  G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*F22;
615  G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22
616  -X1*F2*F32-X2*F3*F12-X3*F1*F22;
617  ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0;
618  }
619  return ranQ2; // MeV^2
620 }
F22
#define F22
Definition: G4ErrorMatrix.cc:902
Q
static double Q[]
Definition: nf_gammaFunctions.cc:91
G4cout
G4GLOB_DLL std::ostream G4cout
F32
#define F32
Definition: G4ErrorMatrix.cc:907
G4endl
#define G4endl
Definition: G4ios.hh:61
F12
#define F12
Definition: G4ErrorMatrix.cc:897
G4double
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

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

GeV/GeV;

Definition at line 788 of file G4ElasticHadrNucleusHE.cc.

789 {
790 // ------ All external kinematical variables are in MeV -------
791 // ------ but internal in GeV !!!! ------
792 
793  G4int NWeight = G4lrint(nistManager->GetAtomicMassAmu(Z));
794 
795  G4double theQ2 = aQ2;
796 
797  // Scattering of proton
798  if(NWeight == 1)
799  {
800  G4double SqrQ2 = std::sqrt(aQ2);
801  G4double valueConstU = hMass2 + protonM2-2*protonM*HadrEnergy - aQ2;
802 
803  G4double MaxT = 4*MomentumCM*MomentumCM;
804 
805  BoundaryTL[0] = MaxT;
806  BoundaryTL[1] = MaxT;
807  BoundaryTL[3] = MaxT;
808  BoundaryTL[4] = MaxT;
809  BoundaryTL[5] = MaxT;
810 
811  G4double dSigPodT;
812 
813  dSigPodT = HadrTot*HadrTot*(1+HadrReIm*HadrReIm)*
814  (
815  Coeff1*G4Exp(-Slope1*SqrQ2)+
816  Coeff2*G4Exp( Slope2*(valueConstU)+aQ2)+
817  (1-Coeff1-Coeff0)*G4Exp(-HadrSlope*aQ2)+
818  +Coeff0*G4Exp(-Slope0*aQ2)
819 // +0.1*(1-std::fabs(CosTh))
820  )*2.568/(16*pi);
821 
822  return dSigPodT;
823  }
824 
825  G4double Stot = HadrTot*MbToGeV2;
826  G4double Bhad = HadrSlope;
827  G4double Asq = 1+HadrReIm*HadrReIm;
828  G4double Rho2 = std::sqrt(Asq);
829  G4double Pnuclp = 0.001;
830  Pnuclp = Pnucl;
831  G4double R12 = R1*R1;
832  G4double R22 = R2*R2;
833  G4double R12B = R12+2*Bhad;
834  G4double R22B = R22+2*Bhad;
835  G4double R12Ap = R12+20;
836  G4double R22Ap = R22+20;
837  G4double R13Ap = R12*R1/R12Ap;
838  G4double R23Ap = R22*R2/R22Ap*Pnuclp;
839  G4double R23dR13 = R23Ap/R13Ap;
840  G4double R12Apd = 2/R12Ap;
841  G4double R22Apd = 2/R22Ap;
842  G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);
843 
844  G4double DDSec1p = (DDSect2+DDSect3*G4Log(0.53*HadrEnergy/R1));
845  G4double DDSec2p = (DDSect2+DDSect3*G4Log(0.53*HadrEnergy/
846  std::sqrt((R12+R22)/2)));
847  G4double DDSec3p = (DDSect2+DDSect3*G4Log(0.53*HadrEnergy/R2));
848 
849  G4double Norm = (R12*R1-Pnucl*R22*R2)*Aeff;
850  G4double Normp = (R12*R1-Pnuclp*R22*R2)*Aeff;
851  G4double R13 = R12*R1/R12B;
852  G4double R23 = Pnucl*R22*R2/R22B;
853  G4double Unucl = Stot/(twopi*Norm)*R13;
854  G4double UnuclScr = Stot/(twopi*Normp)*R13Ap;
855  G4double SinFi = HadrReIm/Rho2;
856  G4double FiH = std::asin(SinFi);
857  G4double N = -1;
858  G4double N2 = R23/R13;
859 
860  G4double ImElasticAmpl0 = 0;
861  G4double ReElasticAmpl0 = 0;
862 
863  G4double exp1;
864  G4double N4;
865  G4double Prod1, Tot1=0, medTot, DTot1, DmedTot;
866  G4int i;
867 
868  for( i=1; i<=NWeight; i++)
869  {
870  N = -N*Unucl*(NWeight-i+1)/i*Rho2;
871  N4 = 1;
872  Prod1 = G4Exp(-theQ2/i*R12B/4)/i*R12B;
873  medTot = R12B/i;
874 
875  for(G4int l=1; l<=i; l++)
876  {
877  exp1 = l/R22B+(i-l)/R12B;
878  N4 = -N4*(i-l+1)/l*N2;
879  Prod1 = Prod1+N4/exp1*G4Exp(-theQ2/(exp1*4));
880  medTot = medTot+N4/exp1;
881  } // end l
882 
883  ReElasticAmpl0 = ReElasticAmpl0+Prod1*N*std::sin(FiH*i);
884  ImElasticAmpl0 = ImElasticAmpl0+Prod1*N*std::cos(FiH*i);
885  Tot1 = Tot1+medTot*N*std::cos(FiH*i);
886  if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001) break;
887  } // i
888 
889  ImElasticAmpl0 = ImElasticAmpl0*pi/2.568; // The amplitude in mB
890  ReElasticAmpl0 = ReElasticAmpl0*pi/2.568; // The amplitude in mB
891  Tot1 = Tot1*twopi/2.568;
892 
893  G4double C1 = R13Ap*R13Ap*0.5*DDSec1p;
894  G4double C2 = 2*R23Ap*R13Ap*0.5*DDSec2p;
895  G4double C3 = R23Ap*R23Ap*0.5*DDSec3p;
896 
897  G4double N1p = 1;
898 
899  G4double Din1 = 0.5;
900 
901  Din1 = 0.5*(C1*G4Exp(-theQ2/8*R12Ap)/2*R12Ap-
902  C2/R12ApdR22Ap*G4Exp(-theQ2/(4*R12ApdR22Ap))+
903  C3*R22Ap/2*G4Exp(-theQ2/8*R22Ap));
904 
905  DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2);
906 
907  G4double exp1p;
908  G4double exp2p;
909  G4double exp3p;
910  G4double N2p;
911  G4double Din2, BinCoeff;
912 
913  BinCoeff = 1;
914 
915  for( i = 1; i<= NWeight-2; i++)
916  {
917  N1p = -N1p*UnuclScr*(NWeight-i-1)/i*Rho2;
918  N2p = 1;
919  Din2 = 0;
920  DmedTot = 0;
921  for(G4int l = 0; l<=i; l++)
922  {
923  if(l == 0) BinCoeff = 1;
924  else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l;
925 
926  exp1 = l/R22B+(i-l)/R12B;
927  exp1p = exp1+R12Apd;
928  exp2p = exp1+R12ApdR22Ap;
929  exp3p = exp1+R22Apd;
930 
931  Din2 = Din2 + N2p*BinCoeff*
932  (C1/exp1p*G4Exp(-theQ2/(4*exp1p))-
933  C2/exp2p*G4Exp(-theQ2/(4*exp2p))+
934  C3/exp3p*G4Exp(-theQ2/(4*exp3p)));
935 
936  DmedTot = DmedTot + N2p*BinCoeff*
937  (C1/exp1p-C2/exp2p+C3/exp3p);
938 
939  N2p = -N2p*R23dR13;
940  } // l
941 
942  Din1 = Din1+Din2*N1p/*Mnoj[i]*//((i+2)*(i+1))*std::cos(FiH*i);
943  DTot1 = DTot1+DmedTot*N1p/*Mnoj[i]*//((i+2)*(i+1))*std::cos(FiH*i);
944 
945  if(std::fabs(Din2*N1p/Din1) < 0.000001) break;
946  } // i
947 
948  Din1 = -Din1*NWeight*(NWeight-1)*4/(Normp*Normp);
949 
950  DTot1 = DTot1*NWeight*(NWeight-1)*4/(Normp*Normp);
951 
952  DTot1 *= 5; // $$$$$$$$$$$$$$$$$$$$$$$$
953 // Din1 *= 0.2; // %%%%%%%%%%%%%%%%%%%%%%% proton
954 // Din1 *= 0.05; // %%%%%%%%%%%%%%%%%%%%%%% pi+
955 // ---------------- dSigma/d|-t|, mb/(GeV/c)^-2 -----------------
956 
957  G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
958  (ImElasticAmpl0+Din1)*
959  (ImElasticAmpl0+Din1))/twopi;
960 
961  Tot1 = Tot1-DTot1;
962  // Tott1 = Tot1*1.0;
963  Dtot11 = DTot1;
964  aAIm = ImElasticAmpl0;
965  aDIm = Din1;
966 
967  return DiffCrSec2; // dSig/d|-t|, mb/(GeV/c)^-2
968 } // function
C2
const double C2
Definition: MCGIDI_KalbachMann.cc:15
C1
const double C1
Definition: MCGIDI_KalbachMann.cc:15
G4int
int G4int
Definition: G4Types.hh:78
C3
#define C3
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:359
N
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4double
double G4double
Definition: G4Types.hh:76
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 465 of file G4ElasticHadrNucleusHE.cc.

467 {
468  //ANDREA: Important notice on this function
469  // For MT we are sharing among threads the G4ElasticData classes
470  // The *call* to this function is proteced
471  // with a mutex in the calling function
472  G4double LineFq2[ONQ2];
473 
474  G4double Rand = G4UniformRand();
475 
476  G4int iNumbQ2 = 0;
477  G4double Q2 = 0.0;
478 
479  G4double ptot2 = plab*plab;
480  G4double ekin = std::sqrt(hMass2 + ptot2) - hMass;
481 
482  if(verboseLevel > 1)
483  G4cout<<"Q2_2: ekin plab "<<ekin<<" "<<plab<<" tmax "<<tmax<<G4endl;
484 
485  // Find closest energy bin
486  G4int NumbOnE;
487  for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ )
488  {
489  if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break;
490  }
491  G4double* dNumbQ2 = pElD->TableQ2;
492 
493  if(NumbOnE >= NENERGY-1) { NumbOnE = NENERGY-2; }
494  G4int index = NumbOnE*ONQ2;
495 
496  // Select kinematics for node energy
497  G4double T = Energy[NumbOnE];
498  hLabMomentum2 = T*(T + 2.*hMass);
499  G4double Q2max = pElD->maxQ2[NumbOnE];
500  G4int length = pElD->dnkE[NumbOnE];
501 
502  // Build vector
503  if(length == 0)
504  {
505  R1 = pElD->R1;
506  R2 = pElD->R2;
507  Aeff = pElD->Aeff;
508  Pnucl = pElD->Pnucl;
509  hLabMomentum = std::sqrt(hLabMomentum2);
510 
511  DefineHadronValues(Z);
512 
513  if(verboseLevel >0)
514  {
515  G4cout<<"1 plab T "<<plab<<" "<<T<<" sigTot B ReIm "
516  <<HadrTot<<" "<<HadrSlope<<" "<<HadrReIm<<G4endl;
517  G4cout<<" R1 R2 Aeff p "<<R1<<" "<<R2<<" "<<Aeff<<" "
518  <<Pnucl<<G4endl;
519  }
520 
521  //pElD->CrossSecMaxQ2[NumbOnE] = 1.0;
522 
523  if(verboseLevel > 1)
524  G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
525  << " length= " << length
526  << " Q2max= " << Q2max
527  << " ekin= " << ekin <<G4endl;
528 
529  pElD->TableCrossSec[index] = 0;
530 
531 
532  dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0];
533 
534  GetHeavyFq2(Z, NumbN, LineFq2); // %%%%%%%%%%%%%%%%%%%%%%%%%
535 
536  for(G4int ii=0; ii<ONQ2; ++ii)
537  {
538  //if(verboseLevel > 2)
539  // G4cout<<" ii LineFq2 "<<ii<<" "<<LineFq2[ii]/LineFq2[ONQ2-1]
540  // <<" dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl;
541 
542  pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1];
543  }
544 
545  pElD->dnkE[NumbOnE] = ONQ2;
546  length = ONQ2;
547  }
548 
549  G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]);
550 
551  for( iNumbQ2 = 1; iNumbQ2<length; ++iNumbQ2)
552  {
553  if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break;
554  }
555  if(iNumbQ2 >= ONQ2) { iNumbQ2 = ONQ2 - 1; }
556  Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
557 
558  if(tmax < Q2max) Q2 *= tmax/Q2max;
559 
560  if(verboseLevel > 1)
561  G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" iNumbQ2= " << iNumbQ2
562  << " rand= " << Rand << G4endl;
563 
564  return Q2;
565 }
ONQ2
static const G4int ONQ2
Definition: G4ElasticHadrNucleusHE.hh:60
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4ElasticData::TableCrossSec
G4double TableCrossSec[NQTABLE]
Definition: G4ElasticHadrNucleusHE.hh:101
G4ElasticData::maxQ2
G4double maxQ2[NENERGY]
Definition: G4ElasticHadrNucleusHE.hh:98
NENERGY
static const G4int NENERGY
Definition: G4ElasticHadrNucleusHE.hh:61
G4ElasticHadrNucleusHE::GetHeavyFq2
G4double GetHeavyFq2(G4int Z, G4int Nucleus, G4double *LineFq2)
Definition: G4ElasticHadrNucleusHE.cc:625
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:571
G4double
double G4double
Definition: G4Types.hh:76
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 1450 of file G4ElasticHadrNucleusHE.cc.

1452 {
1453 
1454  hMass = p->GetPDGMass()/GeV;
1455  hMass2 = hMass*hMass;
1456  hLabMomentum = inLabMom;
1457  hLabMomentum2 = hLabMomentum*hLabMomentum;
1458  HadrEnergy = sqrt(hLabMomentum2+hMass2);
1459 
1460  G4double Rand = G4UniformRand();
1461 
1462  GetKinematics(p, inLabMom);
1463 
1464  G4double Q2 = GetQ2(Rand);
1465 
1466  return Q2;
1467 }
G4ElasticHadrNucleusHE::GetQ2
G4double GetQ2(G4double Ran)
Definition: G4ElasticHadrNucleusHE.cc:1415
p
const char * p
Definition: xmltok.h:285
G4ElasticHadrNucleusHE::GetKinematics
void GetKinematics(const G4ParticleDefinition *aHadron, G4double MomentumH)
Definition: G4ElasticHadrNucleusHE.cc:1356
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 248 of file G4ElasticHadrNucleusHE.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 234 of file G4ElasticHadrNucleusHE.hh.

237 {
238 // G4cout<<" LineInterpol: p1 p2 c1 c2 "<<p1<<" "<<p2<<" "
239 // <<c1<<" "<<c2<<" c "<<c1+(p-p1)*(c2-c1)/(p2-p1)<<G4endl;
240 
241 
242  return c1+(p-p1)*(c2-c1)/(p2-p1);
243 }
c2
static c2_factory< G4double > c2
Definition: G4ScreenedNuclearRecoil.cc:126
p
const char * p
Definition: xmltok.h:285
plottest35.c1
tuple c1
Definition: plottest35.py:14

Here is the caller graph for this function:

void G4ElasticHadrNucleusHE::ModelDescription ( std::ostream &  outFile) const
virtual

Reimplemented from G4HadronElastic.

Definition at line 323 of file G4ElasticHadrNucleusHE.cc.

324 {
325 
326  outFile << "G4ElasticHadrNucleusHE is a hadron-nucleus elastic scattering\n"
327  << "model developed by N. Starkov which uses a Glauber model\n"
328  << "parameterization to calculate the final state. It is valid\n"
329  << "for all hadrons with incident energies above 1 GeV.\n";
330 
331 }
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 363 of file G4ElasticHadrNucleusHE.cc.

366 {
367  G4int Z = iZ;
368  if(Z >= ZMAX) { Z = ZMAX-1; }
369  G4double plab = inLabMom/GeV; // (GeV/c)
370  G4double Q2 = 0;
371 
372  iHadrCode = p->GetPDGEncoding();
373 
374  NumbN = N;
375 
376  if(verboseLevel > 1)
377  {
378  G4cout<< " G4ElasticHadrNucleusHE::SampleT: "
379  << " for " << p->GetParticleName()
380  << " at Z= " << Z << " A= " << N
381  << " plab(GeV)= " << plab
382  << G4endl;
383  }
384  G4int idx;
385 
386  for( idx = 0 ; idx < NHADRONS; idx++)
387  {
388  if(iHadrCode == HadronCode[idx]) break;
389  }
390 
391  // Hadron is not in the list
392 
393  if( idx >= NHADRONS ) return Q2;
394 
395  iHadron = HadronType[idx];
396  iHadrCode = HadronCode[idx];
397 
398  if(Z==1)
399  {
400  hMass = p->GetPDGMass()/GeV;
401  hMass2 = hMass*hMass;
402 
403  G4double T = sqrt(plab*plab+hMass2)-hMass;
404 
405  if(T > 0.4) Q2 = HadronProtonQ2(p, plab);
406 
407  if (verboseLevel>1)
408  G4cout<<" Proton : Q2 "<<Q2<<G4endl;
409  }
410  else
411  {
412  G4AutoLock l(&(eldata_m[idx][Z]));//ANDREA
413  G4ElasticData* ElD1 = SetOfElasticData[idx][Z];
414 
415  // Construct elastic data
416  if(!ElD1)
417  {
418  G4double AWeight = nistManager->GetAtomicMassAmu(Z);
419  ElD1 = new G4ElasticData(p, Z, AWeight, Energy);
420  SetOfElasticData[idx][Z] = ElD1;
421 
422  if(verboseLevel > 1)
423  {
424  G4cout<< " G4ElasticHadrNucleusHE::SampleT: new record " << idx
425  << " for " << p->GetParticleName() << " Z= " << Z
426  << G4endl;
427  }
428  }
429  hMass = ElD1->massGeV;
430  hMass2 = ElD1->mass2GeV2;
431  G4double M = ElD1->massA;
432  G4double M2 = ElD1->massA2;
433  G4double plab2 = plab*plab;
434  G4double Q2max = 4.*plab2*M2/
435  (hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2));
436 
437  // sample scattering
438  G4double T = sqrt(plab2+hMass2)-hMass;
439 
440  if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max);
441 
442  if(verboseLevel > 1)
443  G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< " t/tmax= " << Q2/Q2max <<G4endl;
444  }
445  return Q2*GeV2;
446 }
G4ElasticHadrNucleusHE::HadronProtonQ2
G4double HadronProtonQ2(const G4ParticleDefinition *aHadron, G4double inLabMom)
Definition: G4ElasticHadrNucleusHE.cc:1450
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
ZMAX
static const G4int ZMAX
Definition: G4ElasticHadrNucleusHE.hh:62
G4cout
G4GLOB_DLL std::ostream G4cout
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:359
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
N
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
G4ElasticHadrNucleusHE::HadronNucleusQ2_2
G4double HadronNucleusQ2_2(G4ElasticData *pElD, G4int Z, G4double plabGeV, G4double tmax)
Definition: G4ElasticHadrNucleusHE.cc:465
G4double
double G4double
Definition: G4Types.hh:76
NHADRONS
static const G4int NHADRONS
Definition: G4ElasticHadrNucleusHE.hh:58
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 453 of file G4ElasticHadrNucleusHE.cc.

456 {
457  return SampleInvariantT(p, inLabMom, Z, N);
458 }
N
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
G4ElasticHadrNucleusHE::SampleInvariantT
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4ElasticHadrNucleusHE.cc:363
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

The documentation for this class was generated from the following files: