#include <G4ElasticHadrNucleusHE.hh>
Inheritance diagram for G4ElasticHadrNucleusHE:
Definition at line 108 of file G4ElasticHadrNucleusHE.hh.
| G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE | ( | const G4String & | name = "hElasticGlauber" |
) |
Definition at line 226 of file G4ElasticHadrNucleusHE.cc.
References G4NistManager::Instance(), and G4HadronicInteraction::verboseLevel.
00227 : G4HadronElastic(name) 00228 { 00229 dQ2 = hMass = hMass2 = hLabMomentum = hLabMomentum2 = MomentumCM = HadrEnergy 00230 = R1 = R2 = Pnucl = Aeff = HadrTot = HadrSlope = HadrReIm = TotP = DDSect2 00231 = DDSect3 = ConstU = FmaxT = Slope1 = Slope2 = Coeff1 = Coeff2 = MaxTR 00232 = Slope0 = Coeff0 = aAIm = aDIm = Dtot11 = 0.0; 00233 NumbN = iHadrCode = iHadron = 0; 00234 00235 verboseLevel = 0; 00236 plabLowLimit = 20.0*MeV; 00237 lowestEnergyLimit = 0.0; 00238 //Description(); 00239 00240 MbToGeV2 = 2.568; 00241 sqMbToGeV = 1.602; 00242 Fm2ToGeV2 = 25.68; 00243 GeV2 = GeV*GeV; 00244 protonM = proton_mass_c2/GeV; 00245 protonM2 = protonM*protonM; 00246 00247 BoundaryP[0]=9.0;BoundaryTG[0]=5.0;BoundaryTL[0]=0.; 00248 BoundaryP[1]=20.0;BoundaryTG[1]=1.5;BoundaryTL[1]=0.; 00249 BoundaryP[2]=5.0; BoundaryTG[2]=1.0;BoundaryTL[2]=1.5; 00250 BoundaryP[3]=8.0; BoundaryTG[3]=3.0;BoundaryTL[3]=0.; 00251 BoundaryP[4]=7.0; BoundaryTG[4]=3.0;BoundaryTL[4]=0.; 00252 BoundaryP[5]=5.0; BoundaryTG[5]=2.0;BoundaryTL[5]=0.; 00253 BoundaryP[6]=5.0; BoundaryTG[6]=1.5;BoundaryTL[6]=3.0; 00254 00255 Binom(); 00256 // energy in GeV 00257 Energy[0] = 0.4; 00258 Energy[1] = 0.6; 00259 Energy[2] = 0.8; 00260 LowEdgeEnergy[0] = 0.0; 00261 LowEdgeEnergy[1] = 0.5; 00262 LowEdgeEnergy[2] = 0.7; 00263 G4double e = 1.0; 00264 G4double f = std::pow(10.,0.1); 00265 for(G4int i=3; i<NENERGY; i++) { 00266 Energy[i] = e; 00267 LowEdgeEnergy[i] = e/f; 00268 e *= f*f; 00269 } 00270 nistManager = G4NistManager::Instance(); 00271 00272 // PDG code for hadrons 00273 G4int ic[NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311, 00274 3122,3222,3112,3212,3312,3322,3334, 00275 -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334}; 00276 // internal index 00277 G4int id[NHADRONS] = {2,3,6,0,4,5,4,4,4,5, 00278 0,0,0,0,0,0,0, 00279 1,7,1,1,1,1,1,1,1}; 00280 00281 G4int id1[NHADRONS] = {3,4,1,0,5,6,5,5,5,6, 00282 0,0,0,0,0,0,0, 00283 2,2,2,2,2,2,2,2,2}; 00284 00285 for(G4int j=0; j<NHADRONS; j++) 00286 { 00287 HadronCode[j] = ic[j]; 00288 HadronType[j] = id[j]; 00289 HadronType1[j] = id1[j]; 00290 00291 for(G4int k = 0; k < 93; k++) { SetOfElasticData[j][k] = 0; } 00292 } 00293 }
| G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusHE | ( | ) | [virtual] |
Definition at line 327 of file G4ElasticHadrNucleusHE.cc.
00328 { 00329 for(G4int j = 0; j < NHADRONS; j++) 00330 { 00331 for(G4int k = 0; k < 93; k++) 00332 { 00333 if(SetOfElasticData[j][k]) delete SetOfElasticData[j][k]; 00334 } 00335 } 00336 }
| G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE | ( | const G4ElasticHadrNucleusHE & | ) |
| void G4ElasticHadrNucleusHE::DefineHadronValues | ( | G4int | Z | ) |
Definition at line 948 of file G4ElasticHadrNucleusHE.cc.
References G4cout, G4endl, InterpolateHN(), and G4HadronicInteraction::verboseLevel.
Referenced by GetKinematics(), and HadronNucleusQ2_2().
00949 { 00950 HadrEnergy = std::sqrt(hMass2 + hLabMomentum2); 00951 00952 G4double sHadr = 2.*HadrEnergy*protonM+protonM2+hMass2; 00953 G4double sqrS = std::sqrt(sHadr); 00954 G4double Ecm = 0.5*(sHadr-hMass2+protonM2)/sqrS; 00955 MomentumCM = std::sqrt(Ecm*Ecm-protonM2); 00956 00957 if(verboseLevel>2) 00958 G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron 00959 << " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS 00960 << " plab= " << hLabMomentum 00961 <<" E - m "<<HadrEnergy - hMass<< G4endl; 00962 00963 G4double TotN = 0.0; 00964 G4double logE = std::log(HadrEnergy); 00965 G4double logS = std::log(sHadr); 00966 TotP = 0.0; 00967 00968 switch (iHadron) 00969 { 00970 case 0: // proton, neutron 00971 case 6: 00972 00973 if(hLabMomentum > 10) 00974 TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165); // mb 00975 00976 else 00977 { 00978 // ================== neutron ================ 00979 00981 00982 00983 if( hLabMomentum > 1.4 ) 00984 TotN = 33.3+15.2*(hLabMomentum2-1.35)/ 00985 (std::pow(hLabMomentum,2.37)+0.95); 00986 00987 else if(hLabMomentum > 0.8) 00988 { 00989 G4double A0 = logE + 0.0513; 00990 TotN = 33.0 + 25.5*A0*A0; 00991 } 00992 else 00993 { 00994 G4double A0 = logE - 0.2634; // log(1.3) 00995 TotN = 33.0 + 30.*A0*A0*A0*A0; 00996 } 00997 // ================= proton =============== 00998 // else if(iHadrCode == 2212) 00999 { 01000 if(hLabMomentum >= 1.05) 01001 { 01002 TotP = 39.0+75.*(hLabMomentum-1.2)/ 01003 (hLabMomentum2*hLabMomentum+0.15); 01004 } 01005 01006 else if(hLabMomentum >= 0.7) 01007 { 01008 G4double A0 = logE + 0.3147; 01009 TotP = 23.0 + 40.*A0*A0; 01010 } 01011 else 01012 { 01013 TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5); 01014 } 01015 } 01016 } 01017 01018 // HadrTot = 0.5*(82*TotP+126*TotN)/104; // $$$$$$$$$$$$$$$$$$ 01019 HadrTot = 0.5*(TotP+TotN); 01020 // ................................................... 01021 // Proton slope 01022 if(hLabMomentum >= 2.) HadrSlope = 5.44 + 0.88*logS; 01023 01024 else if(hLabMomentum >= 0.5) HadrSlope = 3.73*hLabMomentum-0.37; 01025 01026 else HadrSlope = 1.5; 01027 01028 // ................................................... 01029 if(hLabMomentum >= 1.2) 01030 HadrReIm = 0.13*(logS - 5.8579332)*std::pow(sHadr,-0.18); 01031 01032 else if(hLabMomentum >= 0.6) 01033 HadrReIm = -75.5*(std::pow(hLabMomentum,0.25)-0.95)/ 01034 (std::pow(3*hLabMomentum,2.2)+1); 01035 01036 else 01037 HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2); 01038 // ................................................... 01039 DDSect2 = 2.2; //mb*GeV-2 01040 DDSect3 = 0.6; //mb*GeV-2 01041 // ================== lambda ================== 01042 if( iHadrCode == 3122) 01043 { 01044 HadrTot *= 0.88; 01045 HadrSlope *=0.85; 01046 } 01047 // ================== sigma + ================== 01048 else if( iHadrCode == 3222) 01049 { 01050 HadrTot *=0.81; 01051 HadrSlope *=0.85; 01052 } 01053 // ================== sigma 0,- ================== 01054 else if(iHadrCode == 3112 || iHadrCode == 3212 ) 01055 { 01056 HadrTot *=0.88; 01057 HadrSlope *=0.85; 01058 } 01059 // =================== xi ================= 01060 else if( iHadrCode == 3312 || iHadrCode == 3322 ) 01061 { 01062 HadrTot *=0.77; 01063 HadrSlope *=0.75; 01064 } 01065 // ================= omega ================= 01066 else if( iHadrCode == 3334) 01067 { 01068 HadrTot *=0.78; 01069 HadrSlope *=0.7; 01070 } 01071 01072 break; 01073 // =========================================================== 01074 case 1: // antiproton 01075 case 7: // antineutron 01076 01077 HadrTot = 5.2+5.2*logE + 123.2/sqrS; // mb 01078 HadrSlope = 8.32+0.57*logS; //(GeV/c)^-2 01079 01080 if( HadrEnergy < 1000 ) 01081 HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8); 01082 else 01083 HadrReIm = 0.6*(logS - 5.8579332)*std::pow(sHadr,-0.25); 01084 01085 DDSect2 = 11; //mb*(GeV/c)^-2 01086 DDSect3 = 3; //mb*(GeV/c)^-2 01087 // ================== lambda ================== 01088 if( iHadrCode == -3122) 01089 { 01090 HadrTot *= 0.88; 01091 HadrSlope *=0.85; 01092 } 01093 // ================== sigma + ================== 01094 else if( iHadrCode == -3222) 01095 { 01096 HadrTot *=0.81; 01097 HadrSlope *=0.85; 01098 } 01099 // ================== sigma 0,- ================== 01100 else if(iHadrCode == -3112 || iHadrCode == -3212 ) 01101 { 01102 HadrTot *=0.88; 01103 HadrSlope *=0.85; 01104 } 01105 // =================== xi ================= 01106 else if( iHadrCode == -3312 || iHadrCode == -3322 ) 01107 { 01108 HadrTot *=0.77; 01109 HadrSlope *=0.75; 01110 } 01111 // ================= omega ================= 01112 else if( iHadrCode == -3334) 01113 { 01114 HadrTot *=0.78; 01115 HadrSlope *=0.7; 01116 } 01117 01118 break; 01119 // ------------------------------------------- 01120 case 2: // pi plus, pi minus 01121 case 3: 01122 01123 if(hLabMomentum >= 3.5) 01124 TotP = 10.6+2.*logE + 25.*std::pow(HadrEnergy,-0.43); // mb 01125 // ========================================= 01126 else if(hLabMomentum >= 1.15) 01127 { 01128 G4double x = (hLabMomentum - 2.55)/0.55; 01129 G4double y = (hLabMomentum - 1.47)/0.225; 01130 TotP = 3.2*std::exp(-x*x) + 12.*std::exp(-y*y) + 27.5; 01131 } 01132 // ========================================= 01133 else if(hLabMomentum >= 0.4) 01134 { 01135 TotP = 88*(logE+0.2877)*(logE+0.2877)+14.0; 01136 } 01137 // ========================================= 01138 else 01139 { 01140 G4double x = (hLabMomentum - 0.29)/0.085; 01141 TotP = 20. + 180.*std::exp(-x*x); 01142 } 01143 // ------------------------------------------- 01144 01145 if(hLabMomentum >= 3.0 ) 01146 TotN = 10.6 + 2.*logE + 30.*std::pow(HadrEnergy,-0.43); // mb 01147 01148 else if(hLabMomentum >= 1.3) 01149 { 01150 G4double x = (hLabMomentum - 2.1)/0.4; 01151 G4double y = (hLabMomentum - 1.4)/0.12; 01152 TotN = 36.1+0.079 - 4.313*logE + 3.*std::exp(-x*x) + 01153 1.5*std::exp(-y*y); 01154 } 01155 else if(hLabMomentum >= 0.65) 01156 { 01157 G4double x = (hLabMomentum - 0.72)/0.06; 01158 G4double y = (hLabMomentum - 1.015)/0.075; 01159 TotN = 36.1 + 10.*std::exp(-x*x) + 24*std::exp(-y*y); 01160 } 01161 else if(hLabMomentum >= 0.37) 01162 { 01163 G4double x = std::log(hLabMomentum/0.48); 01164 TotN = 26. + 110.*x*x; 01165 } 01166 else 01167 { 01168 G4double x = (hLabMomentum - 0.29)/0.07; 01169 TotN = 28.0 + 40.*std::exp(-x*x); 01170 } 01171 HadrTot = (TotP+TotN)/2; 01172 // ........................................ 01173 HadrSlope = 7.28+0.245*logS; // GeV-2 01174 HadrReIm = 0.2*(logS - 4.6051702)*std::pow(sHadr,-0.15); 01175 01176 DDSect2 = 0.7; //mb*GeV-2 01177 DDSect3 = 0.27; //mb*GeV-2 01178 01179 break; 01180 // ========================================================== 01181 case 4: // K plus 01182 01183 HadrTot = 10.6+1.8*logE + 9.0*std::pow(HadrEnergy,-0.55); // mb 01184 if(HadrEnergy>100) HadrSlope = 15.0; 01185 else HadrSlope = 1.0+1.76*logS - 2.84/sqrS; // GeV-2 01186 01187 HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1); 01188 DDSect2 = 0.7; //mb*GeV-2 01189 DDSect3 = 0.21; //mb*GeV-2 01190 break; 01191 // ========================================================= 01192 case 5: // K minus 01193 01194 HadrTot = 10+1.8*logE + 25./sqrS; // mb 01195 HadrSlope = 6.98+0.127*logS; // GeV-2 01196 HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1); 01197 DDSect2 = 0.7; //mb*GeV-2 01198 DDSect3 = 0.27; //mb*GeV-2 01199 break; 01200 } 01201 // ========================================================= 01202 if(verboseLevel>2) 01203 G4cout << "HadrTot= " << HadrTot << " HadrSlope= " << HadrSlope 01204 << " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2 01205 << " DDSect3= " << DDSect3 << G4endl; 01206 01207 if(Z != 1) return; 01208 01209 // Scattering of protons 01210 01211 Coeff0 = Coeff1 = Coeff2 = 0.0; 01212 Slope0 = Slope1 = 1.0; 01213 Slope2 = 5.0; 01214 01215 // data for iHadron=0 01216 static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0}; 01217 static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002}; 01218 static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0}; 01219 static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25}; 01220 static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8}; 01221 01222 // data for iHadron=6,7 01223 static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0}; 01224 static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01}; 01225 static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003}; 01226 static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5}; 01227 static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8}; 01228 01229 // data for iHadron=1 01230 static const G4double EnP[2]={1.5,4.0}; 01231 static const G4double C0P[2]={0.001,0.0005}; 01232 static const G4double C1P[2]={0.003,0.001}; 01233 static const G4double B0P[2]={2.5,4.5}; 01234 static const G4double B1P[2]={1.0,4.0}; 01235 01236 // data for iHadron=2 01237 static const G4double EnPP[4]={1.0,2.0,3.0,4.0}; 01238 static const G4double C0PP[4]={0.0,0.0,0.0,0.0}; 01239 static const G4double C1PP[4]={0.15,0.08,0.02,0.01}; 01240 static const G4double B0PP[4]={1.5,2.8,3.8,3.8}; 01241 static const G4double B1PP[4]={0.8,1.6,3.6,4.6}; 01242 01243 // data for iHadron=3 01244 static const G4double EnPPN[4]={1.0,2.0,3.0,4.0}; 01245 static const G4double C0PPN[4]={0.0,0.0,0.0,0.0}; 01246 static const G4double C1PPN[4]={0.0,0.0,0.0,0.0}; 01247 static const G4double B0PPN[4]={1.5,2.8,3.8,3.8}; 01248 static const G4double B1PPN[4]={0.8,1.6,3.6,4.6}; 01249 01250 // data for iHadron=4 01251 static const G4double EnK[4]={1.4,2.33,3.0,5.0}; 01252 static const G4double C0K[4]={0.0,0.0,0.0,0.0}; 01253 static const G4double C1K[4]={0.01,0.007,0.005,0.003}; 01254 static const G4double B0K[4]={1.5,2.0,3.8,3.8}; 01255 static const G4double B1K[4]={1.6,1.6,1.6,1.6}; 01256 01257 // data for iHadron=5 01258 static const G4double EnKM[2]={1.4,4.0}; 01259 static const G4double C0KM[2]={0.006,0.002}; 01260 static const G4double C1KM[2]={0.00,0.00}; 01261 static const G4double B0KM[2]={2.5,3.5}; 01262 static const G4double B1KM[2]={1.6,1.6}; 01263 01264 switch(iHadron) 01265 { 01266 case 0 : 01267 01268 if(hLabMomentum <BoundaryP[0]) 01269 InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0); 01270 01271 Coeff2 = 0.8/hLabMomentum/hLabMomentum; 01272 break; 01273 01274 case 6 : 01275 01276 if(hLabMomentum < BoundaryP[1]) 01277 InterpolateHN(5,EnN,C0N,C1N,B0N,B1N); 01278 01279 Coeff2 = 0.8/hLabMomentum/hLabMomentum; 01280 break; 01281 01282 case 1 : 01283 case 7 : 01284 if(hLabMomentum < BoundaryP[2]) 01285 InterpolateHN(2,EnP,C0P,C1P,B0P,B1P); 01286 break; 01287 01288 case 2 : 01289 01290 if(hLabMomentum < BoundaryP[3]) 01291 InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1PP); 01292 01293 Coeff2 = 0.02/hLabMomentum; 01294 break; 01295 01296 case 3 : 01297 01298 if(hLabMomentum < BoundaryP[4]) 01299 InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN,B1PPN); 01300 01301 Coeff2 = 0.02/hLabMomentum; 01302 break; 01303 01304 case 4 : 01305 01306 if(hLabMomentum < BoundaryP[5]) 01307 InterpolateHN(4,EnK,C0K,C1K,B0K,B1K); 01308 01309 if(hLabMomentum < 1) Coeff2 = 0.34; 01310 else Coeff2 = 0.34/hLabMomentum2/hLabMomentum; 01311 break; 01312 01313 case 5 : 01314 if(hLabMomentum < BoundaryP[6]) 01315 InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1KM); 01316 01317 if(hLabMomentum < 1) Coeff2 = 0.01; 01318 else Coeff2 = 0.01/hLabMomentum2/hLabMomentum; 01319 break; 01320 } 01321 01322 if(verboseLevel > 2) 01323 G4cout<<" HadrVal : Plasb "<<hLabMomentum 01324 <<" iHadron "<<iHadron<<" HadrTot "<<HadrTot<<G4endl; 01325 }
| void G4ElasticHadrNucleusHE::Description | ( | ) | const [virtual] |
Reimplemented from G4HadronElastic.
Definition at line 296 of file G4ElasticHadrNucleusHE.cc.
References G4HadronicInteraction::GetModelName().
00297 { 00298 char* dirName = getenv("G4PhysListDocDir"); 00299 if (dirName) { 00300 std::ofstream outFile; 00301 G4String outFileName = GetModelName() + ".html"; 00302 G4String pathName = G4String(dirName) + "/" + outFileName; 00303 outFile.open(pathName); 00304 outFile << "<html>\n"; 00305 outFile << "<head>\n"; 00306 00307 outFile << "<title>Description of G4ElasticHadrNucleusHE Model</title>\n"; 00308 outFile << "</head>\n"; 00309 outFile << "<body>\n"; 00310 00311 outFile << "G4ElasticHadrNucleusHE is a hadron-nucleus elastic scattering\n" 00312 << "model developed by N. Starkov which uses a Glauber model\n" 00313 << "parameterization to calculate the final state. It is valid\n" 00314 << "for all hadrons with incident energies above 1 GeV.\n"; 00315 00316 outFile << "</body>\n"; 00317 outFile << "</html>\n"; 00318 outFile.close(); 00319 } 00320 }
Definition at line 269 of file G4ElasticHadrNucleusHE.hh.
00270 { 00271 if ( numN >= numM && numN <= 240) return SetBinom[numN][numM]; 00272 else return 0.; 00273 }
Definition at line 278 of file G4ElasticHadrNucleusHE.hh.
References GetFt().
Referenced by GetQ2().
00279 { 00280 return GetFt(Q2)/FmaxT; 00281 }
Definition at line 1369 of file G4ElasticHadrNucleusHE.cc.
References G4cout, G4endl, and G4HadronicInteraction::verboseLevel.
Referenced by GetDistrFun(), and GetQ2().
01370 { 01371 G4double Fdistr=0; 01372 G4double SqrQ2 = std::sqrt(Q2); 01373 01374 Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*std::exp(ConstU)) 01375 /HadrSlope*(1-std::exp(-HadrSlope*Q2)) 01376 + Coeff0*(1-std::exp(-Slope0*Q2)) 01377 + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1) 01378 + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2)); 01379 01380 if (verboseLevel>1) 01381 G4cout<<"Old: Coeff0 Coeff1 Coeff2 "<<Coeff0<<" " 01382 <<Coeff1<<" "<<Coeff2<<" Slope Slope0 Slope1 Slope2 " 01383 <<HadrSlope<<" "<<Slope0<<" "<<Slope1<<" "<<Slope2 01384 <<" Fdistr "<<Fdistr<<G4endl; 01385 return Fdistr; 01386 }
Definition at line 596 of file G4ElasticHadrNucleusHE.cc.
References G4cout, G4endl, HadrNucDifferCrSec(), and G4HadronicInteraction::verboseLevel.
Referenced by HadronNucleusQ2_2().
00597 { 00598 G4int ii, jj, aSimp; 00599 G4double curQ2, curSec; 00600 G4double curSum = 0.0; 00601 G4double totSum = 0.0; 00602 00603 G4double ddQ2 = dQ2/20; 00604 G4double Q2l = 0; 00605 00606 LineF[0] = 0; 00607 for(ii = 1; ii<ONQ2; ii++) 00608 { 00609 curSum = 0; 00610 aSimp = 4; 00611 00612 for(jj = 0; jj<20; jj++) 00613 { 00614 curQ2 = Q2l+jj*ddQ2; 00615 00616 curSec = HadrNucDifferCrSec(Z, Nucleus, curQ2); 00617 curSum += curSec*aSimp; 00618 00619 if(aSimp > 3) aSimp = 2; 00620 else aSimp = 4; 00621 00622 if(jj == 0 && verboseLevel>2) 00623 G4cout<<" Q2 "<<curQ2<<" AIm "<<aAIm<<" DIm "<<aDIm 00624 <<" Diff "<<curSec<<" totSum "<<totSum<<G4endl; 00625 } 00626 00627 Q2l += dQ2; 00628 curSum *= ddQ2/2.3; // $$$$$$$$$$$$$$$$$$$$$$$ 00629 totSum += curSum; 00630 00631 LineF[ii] = totSum; 00632 00633 if (verboseLevel>2) 00634 G4cout<<" GetHeavy: Q2 dQ2 totSum "<<Q2l<<" "<<dQ2<<" "<<totSum 00635 <<" curSec " 00636 <<curSec<<" totSum "<< totSum<<" DTot " 00637 <<curSum<<G4endl; 00638 } 00639 return totSum; 00640 }
| 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().
01331 { 01332 if (verboseLevel>1) 01333 G4cout<<"1 GetKin.: HadronName MomentumH " 01334 <<aHadron->GetParticleName()<<" "<<MomentumH<<G4endl; 01335 01336 DefineHadronValues(1); 01337 01338 G4double Sh = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV 01339 01340 ConstU = 2*protonM2+2*hMass2-Sh; 01341 01342 G4double MaxT = 4*MomentumCM*MomentumCM; 01343 01344 BoundaryTL[0] = MaxT; //2.0; 01345 BoundaryTL[1] = MaxT; 01346 BoundaryTL[3] = MaxT; 01347 BoundaryTL[4] = MaxT; 01348 BoundaryTL[5] = MaxT; 01349 01350 G4int NumberH=0; 01351 01352 while(iHadrCode!=HadronCode[NumberH]) NumberH++; 01353 01354 NumberH = HadronType1[NumberH]; 01355 01356 if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH]; 01357 else MaxTR = BoundaryTG[NumberH]; 01358 01359 if (verboseLevel>1) 01360 G4cout<<"3 GetKin. : NumberH "<<NumberH 01361 <<" Bound.P[NumberH] "<<BoundaryP[NumberH] 01362 <<" Bound.TL[NumberH] "<<BoundaryTL[NumberH] 01363 <<" Bound.TG[NumberH] "<<BoundaryTG[NumberH] 01364 <<" MaxT MaxTR "<<MaxT<<" "<<MaxTR<<G4endl; 01365 01366 // GetParametersHP(aHadron, MomentumH); 01367 }
Definition at line 646 of file G4ElasticHadrNucleusHE.cc.
References G4cout, G4endl, G4INCL::Math::pi, and G4HadronicInteraction::verboseLevel.
00648 { 00649 // Scattering of proton 00650 if(Z == 1) 00651 { 00652 G4double SqrQ2 = std::sqrt(Q2); 00653 G4double valueConstU = 2.*(hMass2 + protonM2) - Q2; 00654 00655 G4double y = (1.-Coeff1-Coeff0)/HadrSlope*(1.-std::exp(-HadrSlope*Q2)) 00656 + Coeff0*(1.-std::exp(-Slope0*Q2)) 00657 + Coeff2/Slope2*std::exp(Slope2*valueConstU)*(std::exp(Slope2*Q2)-1.) 00658 + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2)); 00659 00660 return y; 00661 } 00662 00663 // The preparing of probability function 00664 00665 G4double prec = Nucleus > 208 ? 1.0e-7 : 1.0e-6; 00666 00667 G4double Stot = HadrTot*MbToGeV2; // Gev^-2 00668 G4double Bhad = HadrSlope; // GeV^-2 00669 G4double Asq = 1+HadrReIm*HadrReIm; 00670 G4double Rho2 = std::sqrt(Asq); 00671 00672 // if(verboseLevel >1) 00673 G4cout<<" Fq2 Before for i Tot B Im "<<HadrTot<<" "<<HadrSlope<<" " 00674 <<HadrReIm<<G4endl; 00675 00676 if(verboseLevel > 1) { 00677 G4cout << "GetFq2: Stot= " << Stot << " Bhad= " << Bhad 00678 <<" Im "<<HadrReIm 00679 << " Asq= " << Asq << G4endl; 00680 G4cout << "R1= " << R1 << " R2= " << R2 << " Pnucl= " << Pnucl <<G4endl; 00681 } 00682 G4double R12 = R1*R1; 00683 G4double R22 = R2*R2; 00684 G4double R12B = R12+2*Bhad; 00685 G4double R22B = R22+2*Bhad; 00686 00687 G4double Norm = (R12*R1-Pnucl*R22*R2); // HP->Aeff; 00688 00689 G4double R13 = R12*R1/R12B; 00690 G4double R23 = Pnucl*R22*R2/R22B; 00691 G4double Unucl = Stot/twopi/Norm*R13; 00692 G4double UnucRho2 = -Unucl*Rho2; 00693 00694 G4double FiH = std::asin(HadrReIm/Rho2); 00695 G4double NN2 = R23/R13; 00696 00697 if(verboseLevel > 2) 00698 G4cout << "UnucRho2= " << UnucRho2 << " FiH= " << FiH << " NN2= " << NN2 00699 << " Norm= " << Norm << G4endl; 00700 00701 G4double dddd; 00702 00703 G4double Prod0 = 0; 00704 G4double N1 = -1.0; 00705 //G4double Tot0 = 0; 00706 G4double exp1; 00707 00708 G4double Prod3 ; 00709 G4double exp2 ; 00710 G4double N4, N5, N2, Prod1, Prod2; 00711 G4int i1, i2, j1, j2; 00712 00713 for(i1 = 1; i1<= Nucleus; i1++) 00714 { 00715 N1 = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2; 00716 Prod1 = 0; 00717 //Tot0 = 0; 00718 N2 = -1; 00719 00720 for(i2 = 1; i2<=Nucleus; i2++) 00721 { 00722 N2 = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2; 00723 Prod2 = 0; 00724 N5 = -1/NN2; 00725 for(j2=0; j2<= i2; j2++) 00726 { 00727 Prod3 = 0; 00728 exp2 = 1/(j2/R22B+(i2-j2)/R12B); 00729 N5 = -N5*NN2; 00730 N4 = -1/NN2; 00731 for(j1=0; j1<=i1; j1++) 00732 { 00733 exp1 = 1/(j1/R22B+(i1-j1)/R12B); 00734 dddd = exp1+exp2; 00735 N4 = -N4*NN2; 00736 Prod3 = Prod3+N4*exp1*exp2* 00737 (1-std::exp(-Q2*dddd/4))/dddd*4*SetBinom[i1][j1]; 00738 } // j1 00739 Prod2 = Prod2 +Prod3*N5*SetBinom[i2][j2]; 00740 } // j2 00741 Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2)); 00742 00743 if (std::fabs(Prod2*N2/Prod1)<prec) break; 00744 } // i2 00745 Prod0 = Prod0 + Prod1*N1; 00746 if(std::fabs(N1*Prod1/Prod0) < prec) break; 00747 00748 } // i1 00749 00750 Prod0 *= 0.25*pi/MbToGeV2; // This is in mb 00751 00752 if(verboseLevel>1) 00753 G4cout << "GetLightFq2 Z= " << Z << " A= " << Nucleus 00754 <<" Q2= " << Q2 << " Res= " << Prod0 << G4endl; 00755 return Prod0; 00756 }
Definition at line 1388 of file G4ElasticHadrNucleusHE.cc.
References GetDistrFun(), and GetFt().
Referenced by HadronProtonQ2().
01389 { 01390 G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta; 01391 G4double Q2=0; 01392 01393 FmaxT = GetFt(MaxTR); 01394 delta = GetDistrFun(DDD0)-Ran; 01395 01396 while(std::fabs(delta) > 0.0001) 01397 { 01398 if(delta>0) 01399 { 01400 DDD2 = DDD0; 01401 DDD0 = (DDD0+DDD1)*0.5; 01402 } 01403 else if(delta<0) 01404 { 01405 DDD1 = DDD0; 01406 DDD0 = (DDD0+DDD2)*0.5; 01407 } 01408 delta = GetDistrFun(DDD0)-Ran; 01409 } 01410 01411 Q2 = DDD0; 01412 01413 return Q2; 01414 }
Definition at line 542 of file G4ElasticHadrNucleusHE.cc.
References F12, F22, F32, G4cout, G4endl, and G4HadronicInteraction::verboseLevel.
Referenced by HadronNucleusQ2_2().
00544 { 00545 G4double ranQ2; 00546 00547 G4double F1 = F[kk-2]; 00548 G4double F2 = F[kk-1]; 00549 G4double F3 = F[kk]; 00550 G4double X1 = Q[kk-2]; 00551 G4double X2 = Q[kk-1]; 00552 G4double X3 = Q[kk]; 00553 00554 if(verboseLevel > 2) 00555 G4cout << "GetQ2_2 kk= " << kk << " X2= " << X2 << " X3= " << X3 00556 << " F2= " << F2 << " F3= " << F3 << " Rndm= " << ranUni << G4endl; 00557 00558 if(kk == 1 || kk == 0) 00559 { 00560 F1 = F[0]; 00561 F2 = F[1]; 00562 F3 = F[2]; 00563 X1 = Q[0]; 00564 X2 = Q[1]; 00565 X3 = Q[2]; 00566 } 00567 00568 G4double F12 = F1*F1; 00569 G4double F22 = F2*F2; 00570 G4double F32 = F3*F3; 00571 00572 G4double D0 = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3; 00573 00574 if(verboseLevel > 2) 00575 G4cout << " X1= " << X1 << " F1= " << F1 << " D0= " 00576 << D0 << G4endl; 00577 00578 if(std::abs(D0) < 0.00000001) 00579 { 00580 ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2); 00581 } 00582 else 00583 { 00584 G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1; 00585 G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*F22; 00586 G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22 00587 -X1*F2*F32-X2*F3*F12-X3*F1*F22; 00588 ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0; 00589 } 00590 return ranQ2; // MeV^2 00591 }
Definition at line 759 of file G4ElasticHadrNucleusHE.cc.
References C1, C3, G4NistManager::GetAtomicMassAmu(), and G4INCL::Math::pi.
Referenced by GetHeavyFq2().
00760 { 00761 // ------ All external kinematical variables are in MeV ------- 00762 // ------ but internal in GeV !!!! ------ 00763 00764 G4int NWeight = int( nistManager->GetAtomicMassAmu(Z) + 0.5 ); 00765 00766 G4double theQ2 = aQ2; 00767 00768 // Scattering of proton 00769 if(NWeight == 1) 00770 { 00771 G4double SqrQ2 = std::sqrt(aQ2); 00772 G4double valueConstU = hMass2 + protonM2-2*protonM*HadrEnergy - aQ2; 00773 00774 G4double MaxT = 4*MomentumCM*MomentumCM; 00775 00776 BoundaryTL[0] = MaxT; 00777 BoundaryTL[1] = MaxT; 00778 BoundaryTL[3] = MaxT; 00779 BoundaryTL[4] = MaxT; 00780 BoundaryTL[5] = MaxT; 00781 00782 G4double dSigPodT; 00783 00784 dSigPodT = HadrTot*HadrTot*(1+HadrReIm*HadrReIm)* 00785 ( 00786 Coeff1*std::exp(-Slope1*SqrQ2)+ 00787 Coeff2*std::exp( Slope2*(valueConstU)+aQ2)+ 00788 (1-Coeff1-Coeff0)*std::exp(-HadrSlope*aQ2)+ 00789 +Coeff0*std::exp(-Slope0*aQ2) 00790 // +0.1*(1-std::fabs(CosTh)) 00791 )/16/3.1416*2.568; 00792 00793 return dSigPodT; 00794 } 00795 00796 G4double Stot = HadrTot*MbToGeV2; 00797 G4double Bhad = HadrSlope; 00798 G4double Asq = 1+HadrReIm*HadrReIm; 00799 G4double Rho2 = std::sqrt(Asq); 00800 G4double Pnuclp = 0.001; 00801 Pnuclp = Pnucl; 00802 G4double R12 = R1*R1; 00803 G4double R22 = R2*R2; 00804 G4double R12B = R12+2*Bhad; 00805 G4double R22B = R22+2*Bhad; 00806 G4double R12Ap = R12+20; 00807 G4double R22Ap = R22+20; 00808 G4double R13Ap = R12*R1/R12Ap; 00809 G4double R23Ap = R22*R2/R22Ap*Pnuclp; 00810 G4double R23dR13 = R23Ap/R13Ap; 00811 G4double R12Apd = 2/R12Ap; 00812 G4double R22Apd = 2/R22Ap; 00813 G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd); 00814 00815 G4double DDSec1p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R1/4)); 00816 G4double DDSec2p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/ 00817 std::sqrt((R12+R22)/2)/4)); 00818 G4double DDSec3p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R2/4)); 00819 00820 G4double Norm = (R12*R1-Pnucl*R22*R2)*Aeff; 00821 G4double Normp = (R12*R1-Pnuclp*R22*R2)*Aeff; 00822 G4double R13 = R12*R1/R12B; 00823 G4double R23 = Pnucl*R22*R2/R22B; 00824 G4double Unucl = Stot/twopi/Norm*R13; 00825 G4double UnuclScr = Stot/twopi/Normp*R13Ap; 00826 G4double SinFi = HadrReIm/Rho2; 00827 G4double FiH = std::asin(SinFi); 00828 G4double N = -1; 00829 G4double N2 = R23/R13; 00830 00831 G4double ImElasticAmpl0 = 0; 00832 G4double ReElasticAmpl0 = 0; 00833 00834 G4double exp1; 00835 G4double N4; 00836 G4double Prod1, Tot1=0, medTot, DTot1, DmedTot; 00837 G4int i; 00838 00839 for( i=1; i<=NWeight; i++) 00840 { 00841 N = -N*Unucl*(NWeight-i+1)/i*Rho2; 00842 N4 = 1; 00843 Prod1 = std::exp(-theQ2/i*R12B/4)/i*R12B; 00844 medTot = R12B/i; 00845 00846 for(G4int l=1; l<=i; l++) 00847 { 00848 exp1 = l/R22B+(i-l)/R12B; 00849 N4 = -N4*(i-l+1)/l*N2; 00850 Prod1 = Prod1+N4/exp1*std::exp(-theQ2/exp1/4); 00851 medTot = medTot+N4/exp1; 00852 } // end l 00853 00854 ReElasticAmpl0 = ReElasticAmpl0+Prod1*N*std::sin(FiH*i); 00855 ImElasticAmpl0 = ImElasticAmpl0+Prod1*N*std::cos(FiH*i); 00856 Tot1 = Tot1+medTot*N*std::cos(FiH*i); 00857 if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001) break; 00858 } // i 00859 00860 ImElasticAmpl0 = ImElasticAmpl0*pi/2.568; // The amplitude in mB 00861 ReElasticAmpl0 = ReElasticAmpl0*pi/2.568; // The amplitude in mB 00862 Tot1 = Tot1*twopi/2.568; 00863 00864 G4double C1 = R13Ap*R13Ap/2*DDSec1p; 00865 G4double C2 = 2*R23Ap*R13Ap/2*DDSec2p; 00866 G4double C3 = R23Ap*R23Ap/2*DDSec3p; 00867 00868 G4double N1p = 1; 00869 00870 G4double Din1 = 0.5; 00871 00872 Din1 = 0.5*(C1*std::exp(-theQ2/8*R12Ap)/2*R12Ap- 00873 C2/R12ApdR22Ap*std::exp(-theQ2/4/R12ApdR22Ap)+ 00874 C3*R22Ap/2*std::exp(-theQ2/8*R22Ap)); 00875 00876 DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2); 00877 00878 G4double exp1p; 00879 G4double exp2p; 00880 G4double exp3p; 00881 G4double N2p; 00882 G4double Din2, BinCoeff; 00883 00884 BinCoeff = 1; 00885 00886 for( i = 1; i<= NWeight-2; i++) 00887 { 00888 N1p = -N1p*UnuclScr*(NWeight-i-1)/i*Rho2; 00889 N2p = 1; 00890 Din2 = 0; 00891 DmedTot = 0; 00892 for(G4int l = 0; l<=i; l++) 00893 { 00894 if(l == 0) BinCoeff = 1; 00895 else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l; 00896 00897 exp1 = l/R22B+(i-l)/R12B; 00898 exp1p = exp1+R12Apd; 00899 exp2p = exp1+R12ApdR22Ap; 00900 exp3p = exp1+R22Apd; 00901 00902 Din2 = Din2 + N2p*BinCoeff* 00903 (C1/exp1p*std::exp(-theQ2/4/exp1p)- 00904 C2/exp2p*std::exp(-theQ2/4/exp2p)+ 00905 C3/exp3p*std::exp(-theQ2/4/exp3p)); 00906 00907 DmedTot = DmedTot + N2p*BinCoeff* 00908 (C1/exp1p-C2/exp2p+C3/exp3p); 00909 00910 N2p = -N2p*R23dR13; 00911 } // l 00912 00913 Din1 = Din1+Din2*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i); 00914 DTot1 = DTot1+DmedTot*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i); 00915 00916 if(std::fabs(Din2*N1p/Din1) < 0.000001) break; 00917 } // i 00918 00919 Din1 = -Din1*NWeight*(NWeight-1) 00920 /2/pi/Normp/2/pi/Normp*16*pi*pi; 00921 00922 DTot1 = DTot1*NWeight*(NWeight-1) 00923 /2/pi/Normp/2/pi/Normp*16*pi*pi; 00924 00925 DTot1 *= 5; // $$$$$$$$$$$$$$$$$$$$$$$$ 00926 // Din1 *= 0.2; // %%%%%%%%%%%%%%%%%%%%%%% proton 00927 // Din1 *= 0.05; // %%%%%%%%%%%%%%%%%%%%%%% pi+ 00928 // ---------------- dSigma/d|-t|, mb/(GeV/c)^-2 ----------------- 00929 00930 G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+ 00931 (ImElasticAmpl0+Din1)* 00932 (ImElasticAmpl0+Din1))*2/4/pi; 00933 00934 Tot1 = Tot1-DTot1; 00935 // Tott1 = Tot1*1.0; 00936 Dtot11 = DTot1; 00937 aAIm = ImElasticAmpl0; 00938 aDIm = Din1; 00939 00940 return DiffCrSec2*1.0; // dSig/d|-t|, mb/(GeV/c)^-2 00941 } // function
| 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().
00444 { 00445 G4double LineFq2[ONQ2]; 00446 00447 G4double Rand = G4UniformRand(); 00448 00449 G4int iNumbQ2 = 0; 00450 G4double Q2 = 0.0; 00451 00452 G4double ptot2 = plab*plab; 00453 G4double ekin = std::sqrt(hMass2 + ptot2) - hMass; 00454 00455 if(verboseLevel > 1) 00456 G4cout<<"Q2_2: ekin plab "<<ekin<<" "<<plab<<" tmax "<<tmax<<G4endl; 00457 00458 // Find closest energy bin 00459 G4int NumbOnE; 00460 for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ ) 00461 { 00462 if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break; 00463 } 00464 G4double* dNumbQ2 = pElD->TableQ2; 00465 00466 G4int index = NumbOnE*ONQ2; 00467 00468 // Select kinematics for node energy 00469 G4double T = Energy[NumbOnE]; 00470 hLabMomentum2 = T*(T + 2.*hMass); 00471 G4double Q2max = pElD->maxQ2[NumbOnE]; 00472 G4int length = pElD->dnkE[NumbOnE]; 00473 00474 // Build vector 00475 if(length == 0) 00476 { 00477 R1 = pElD->R1; 00478 R2 = pElD->R2; 00479 Aeff = pElD->Aeff; 00480 Pnucl = pElD->Pnucl; 00481 hLabMomentum = std::sqrt(hLabMomentum2); 00482 00483 DefineHadronValues(Z); 00484 00485 if(verboseLevel >0) 00486 { 00487 G4cout<<"1 plab T "<<plab<<" "<<T<<" sigTot B ReIm " 00488 <<HadrTot<<" "<<HadrSlope<<" "<<HadrReIm<<G4endl; 00489 G4cout<<" R1 R2 Aeff p "<<R1<<" "<<R2<<" "<<Aeff<<" " 00490 <<Pnucl<<G4endl; 00491 } 00492 00493 //pElD->CrossSecMaxQ2[NumbOnE] = 1.0; 00494 00495 if(verboseLevel > 1) 00496 G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE 00497 << " length= " << length 00498 << " Q2max= " << Q2max 00499 << " ekin= " << ekin <<G4endl; 00500 00501 pElD->TableCrossSec[index] = 0; 00502 00503 00504 dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0]; 00505 00506 GetHeavyFq2(Z, NumbN, LineFq2); // %%%%%%%%%%%%%%%%%%%%%%%%% 00507 00508 for(G4int ii=0; ii<ONQ2; ii++) 00509 { 00510 //if(verboseLevel > 2) 00511 // G4cout<<" ii LineFq2 "<<ii<<" "<<LineFq2[ii]/LineFq2[ONQ2-1] 00512 // <<" dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl; 00513 00514 pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1]; 00515 } 00516 00517 pElD->dnkE[NumbOnE] = ONQ2; 00518 length = ONQ2; 00519 } 00520 00521 G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]); 00522 00523 for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ ) 00524 { 00525 if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break; 00526 } 00527 Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand); 00528 00529 if(tmax < Q2max) Q2 *= tmax/Q2max; 00530 00531 if(verboseLevel > 1) 00532 G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" iNumbQ2= " << iNumbQ2 00533 << " rand= " << Rand << G4endl; 00534 00535 return Q2; 00536 }
| G4double G4ElasticHadrNucleusHE::HadronProtonQ2 | ( | const G4ParticleDefinition * | aHadron, | |
| G4double | inLabMom | |||
| ) |
Definition at line 1417 of file G4ElasticHadrNucleusHE.cc.
References G4UniformRand, GetKinematics(), G4ParticleDefinition::GetPDGMass(), and GetQ2().
Referenced by SampleInvariantT().
01419 { 01420 01421 hMass = p->GetPDGMass()/GeV; 01422 hMass2 = hMass*hMass; 01423 hLabMomentum = inLabMom; 01424 hLabMomentum2 = hLabMomentum*hLabMomentum; 01425 HadrEnergy = sqrt(hLabMomentum2+hMass2); 01426 01427 G4double Rand = G4UniformRand(); 01428 01429 GetKinematics(p, inLabMom); 01430 01431 G4double Q2 = GetQ2(Rand); 01432 01433 return Q2; 01434 }
| 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().
Referenced by DefineHadronValues().
00250 { 00251 G4int i; 00252 00253 for(i=1; i<n; i++) if(hLabMomentum <= EnP[i]) break; 00254 00255 if(i == n) i = n - 1; 00256 00257 Coeff0 = LineInterpol(EnP[i], EnP[i-1], C0P[i], C0P[i-1], hLabMomentum); 00258 Coeff1 = LineInterpol(EnP[i], EnP[i-1], C1P[i], C1P[i-1], hLabMomentum); 00259 Slope0 = LineInterpol(EnP[i], EnP[i-1], B0P[i], B0P[i-1], hLabMomentum); 00260 Slope1 = LineInterpol(EnP[i], EnP[i-1], B1P[i], B1P[i-1], hLabMomentum); 00261 00262 // G4cout<<" InterpolHN: n i "<<n<<" "<<i<<" Mom " 00263 // <<hLabMomentum<<G4endl; 00264 }
| 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().
00236 { 00237 // G4cout<<" LineInterpol: p1 p2 c1 c2 "<<p1<<" "<<p2<<" " 00238 // <<c1<<" "<<c2<<" c "<<c1+(p-p1)*(c2-c1)/(p2-p1)<<G4endl; 00239 00240 00241 return c1+(p-p1)*(c2-c1)/(p2-p1); 00242 }
| 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(), HadronNucleusQ2_2(), HadronProtonQ2(), G4ElasticData::mass2GeV2, G4ElasticData::massA, G4ElasticData::massA2, G4ElasticData::massGeV, and G4HadronicInteraction::verboseLevel.
Referenced by SampleT().
00346 { 00347 G4double plab = inLabMom/GeV; // (GeV/c) 00348 G4double Q2 = 0; 00349 00350 iHadrCode = p->GetPDGEncoding(); 00351 00352 NumbN = N; 00353 00354 if(verboseLevel > 1) 00355 { 00356 G4cout<< " G4ElasticHadrNucleusHE::SampleT: " 00357 << " for " << p->GetParticleName() 00358 << " at Z= " << Z << " A= " << N 00359 << " plab(GeV)= " << plab 00360 << G4endl; 00361 } 00362 G4int idx; 00363 00364 for( idx = 0 ; idx < NHADRONS; idx++) 00365 { 00366 if(iHadrCode == HadronCode[idx]) break; 00367 } 00368 00369 // Hadron is not in the list 00370 00371 if( idx >= NHADRONS ) return Q2; 00372 00373 iHadron = HadronType[idx]; 00374 iHadrCode = HadronCode[idx]; 00375 00376 if(Z==1) 00377 { 00378 hMass = p->GetPDGMass()/GeV; 00379 hMass2 = hMass*hMass; 00380 00381 G4double T = sqrt(plab*plab+hMass2)-hMass; 00382 00383 if(T > 0.4) Q2 = HadronProtonQ2(p, plab); 00384 00385 if (verboseLevel>1) 00386 G4cout<<" Proton : Q2 "<<Q2<<G4endl; 00387 } 00388 else 00389 { 00390 G4ElasticData* ElD1 = SetOfElasticData[idx][Z]; 00391 00392 // Construct elastic data 00393 if(!ElD1) 00394 { 00395 G4double AWeight = nistManager->GetAtomicMassAmu(Z); 00396 ElD1 = new G4ElasticData(p, Z, AWeight, Energy); 00397 SetOfElasticData[idx][Z] = ElD1; 00398 00399 if(verboseLevel > 1) 00400 { 00401 G4cout<< " G4ElasticHadrNucleusHE::SampleT: new record " << idx 00402 << " for " << p->GetParticleName() << " Z= " << Z 00403 << G4endl; 00404 } 00405 } 00406 hMass = ElD1->massGeV; 00407 hMass2 = ElD1->mass2GeV2; 00408 G4double M = ElD1->massA; 00409 G4double M2 = ElD1->massA2; 00410 G4double plab2 = plab*plab; 00411 G4double Q2max = 4.*plab2*M2/ 00412 (hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2)); 00413 00414 // sample scattering 00415 G4double T = sqrt(plab2+hMass2)-hMass; 00416 00417 if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max); 00418 00419 if(verboseLevel > 1) 00420 G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< " t/tmax= " << Q2/Q2max <<G4endl; 00421 } 00422 return Q2*GeV2; 00423 }
| G4double G4ElasticHadrNucleusHE::SampleT | ( | const G4ParticleDefinition * | p, | |
| G4double | plab, | |||
| G4int | Z, | |||
| G4int | A | |||
| ) |
Definition at line 430 of file G4ElasticHadrNucleusHE.cc.
References SampleInvariantT().
00433 { 00434 return SampleInvariantT(p, inLabMom, Z, N); 00435 }
1.4.7