#include <G4QANuANuNuclearCrossSection.hh>
Inheritance diagram for G4QANuANuNuclearCrossSection:
Public Member Functions | |
~G4QANuANuNuclearCrossSection () | |
G4double | ThresholdEnergy (G4int Z, G4int N, G4int PDG=-14) |
virtual G4double | GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=0) |
G4double | CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum) |
G4int | GetExchangePDGCode () |
G4double | GetDirectPart (G4double Q2) |
G4double | GetNPartons (G4double Q2) |
G4double | GetQEL_ExchangeQ2 () |
G4double | GetNQE_ExchangeQ2 () |
G4double | GetLastTOTCS () |
G4double | GetLastQELCS () |
Static Public Member Functions | |
static G4VQCrossSection * | GetPointer () |
Protected Member Functions | |
G4QANuANuNuclearCrossSection () |
Definition at line 49 of file G4QANuANuNuclearCrossSection.hh.
G4QANuANuNuclearCrossSection::G4QANuANuNuclearCrossSection | ( | ) | [inline, protected] |
G4QANuANuNuclearCrossSection::~G4QANuANuNuclearCrossSection | ( | ) |
G4double G4QANuANuNuclearCrossSection::CalculateCrossSection | ( | G4bool | CS, | |
G4int | F, | |||
G4int | I, | |||
G4int | PDG, | |||
G4int | Z, | |||
G4int | N, | |||
G4double | Momentum | |||
) | [virtual] |
Implements G4VQCrossSection.
Definition at line 262 of file G4QANuANuNuclearCrossSection.cc.
References G4cerr, G4cout, and G4endl.
Referenced by GetCrossSection().
00264 { 00265 static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2 00266 static const G4int nE=65; // !! If change this, change it in GetFunctions() (*.hh) !! 00267 static const G4int mL=nE-1; 00268 static const G4double EMi=0.; // Universal threshold of the reaction in GeV 00269 static const G4double EMa=300.; // Maximum tabulated Energy of nu_mu in GeV 00270 // *** Begin of the Associative memory for acceleration of the cross section calculations 00271 static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional) 00272 static G4bool first=true; // Flag of initialization of the energy axis 00273 // *** End of Static Definitions (Associative Memory) *** 00274 //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon 00275 //G4double TotEnergy2=Momentum; 00276 onlyCS=CS; // Flag to calculate only CS (not TX & QE) 00277 lastE=Momentum/GeV; // Kinetic energy of the muon neutrino (in GeV!) 00278 if (lastE<=EMi) // Energy is below the minimum energy in the table 00279 { 00280 lastE=0.; 00281 lastSig=0.; 00282 return 0.; 00283 } 00284 G4int Z=targZ; // New Z, which can change the sign 00285 if(F<=0) // This isotope was not the last used isotop 00286 { 00287 if(F<0) // This isotope was found in DAMDB =-------=> RETRIEVE 00288 { 00289 lastTX =(*TX)[I]; // Pointer to the prepared TX function (same isotope) 00290 lastQE =(*QE)[I]; // Pointer to the prepared QE function (same isotope) 00291 } 00292 else // This isotope wasn't calculated previously => CREATE 00293 { 00294 if(first) 00295 { 00296 lastEN = new G4double[nE]; // This must be done only once! 00297 Z=-Z; // To explain GetFunctions that E-axis must be filled 00298 first=false; // To make it only once 00299 } 00300 lastTX = new G4double[nE]; // Allocate memory for the new TX function 00301 lastQE = new G4double[nE]; // Allocate memory for the new QE function 00302 G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK) 00303 if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl; 00304 // *** The synchronization check *** 00305 G4int sync=TX->size(); 00306 if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl; 00307 TX->push_back(lastTX); 00308 QE->push_back(lastQE); 00309 } // End of creation of the new set of parameters 00310 } // End of parameters udate 00311 // =-----------------= NOW Calculate the Cross Section =---------------= 00312 if (lastE<=EMi) // Check that antiNuEnergy is higher than ThreshE 00313 { 00314 lastE=0.; 00315 lastSig=0.; 00316 return 0.; 00317 } 00318 if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization 00319 { 00320 G4int chk=1; 00321 G4int ran=mL/2; 00322 G4int sep=ran; // as a result = an index of the left edge of the interval 00323 while(ran>=2) 00324 { 00325 G4int newran=ran/2; 00326 if(lastE<=lastEN[sep]) sep-=newran; 00327 else sep+=newran; 00328 ran=newran; 00329 chk=chk+chk; 00330 } 00331 if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl; 00332 G4double lowE=lastEN[sep]; 00333 G4double highE=lastEN[sep+1]; 00334 G4double lowTX=lastTX[sep]; 00335 if(lastE<lowE||sep>=mL||lastE>highE) 00336 G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE 00337 <<", sep="<<sep<<", mL="<<mL<<G4endl; 00338 lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E 00339 if(!onlyCS) // Skip the differential cross-section parameters 00340 { 00341 G4double lowQE=lastQE[sep]; 00342 lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE; 00343 #ifdef pdebug 00344 G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl; 00345 #endif 00346 } 00347 } 00348 else 00349 { 00350 lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking... 00351 lastQEL=lastQE[mL]; 00352 } 00353 if(lastQEL<0.) lastQEL = 0.; 00354 if(lastSig<0.) lastSig = 0.; 00355 // The cross-sections are expected to be in mb 00356 lastSig*=mb38; 00357 if(!onlyCS) lastQEL*=mb38; 00358 return lastSig; 00359 }
G4double G4QANuANuNuclearCrossSection::GetCrossSection | ( | G4bool | fCS, | |
G4double | pMom, | |||
G4int | tgZ, | |||
G4int | tgN, | |||
G4int | pPDG = 0 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 90 of file G4QANuANuNuclearCrossSection.cc.
References CalculateCrossSection(), G4cout, G4endl, ThresholdEnergy(), and G4VQCrossSection::tolerance.
00092 { 00093 static G4int j; // A#0f records found in DB for this projectile 00094 static std::vector <G4int> colPDG;// Vector of the projectile PDG code 00095 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) 00096 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) 00097 static std::vector <G4double> colP; // Vector of last momenta for the reaction 00098 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction 00099 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction 00100 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** 00101 G4double pEn=pMom; 00102 #ifdef debug 00103 G4cout<<"G4QAMNCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN 00104 <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz=" 00105 <<colN.size()<<G4endl; 00106 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00107 #endif 00108 if(pPDG!=-14) 00109 { 00110 #ifdef debug 00111 G4cout<<"G4QAMNCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl; 00112 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00113 #endif 00114 return 0.; // projectile PDG=0 is a mistake (?!) @@ 00115 } 00116 G4bool in=false; // By default the isotope must be found in the AMDB 00117 if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope 00118 { 00119 in = false; // By default the isotope haven't be found in AMDB 00120 lastP = 0.; // New momentum history (nothing to compare with) 00121 lastPDG = pPDG; // The last PDG of the projectile 00122 lastN = tgN; // The last N of the calculated nucleus 00123 lastZ = tgZ; // The last Z of the calculated nucleus 00124 lastI = colN.size(); // Size of the Associative Memory DB in the heap 00125 j = 0; // A#0f records found in DB for this projectile 00126 if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found 00127 { // The nucleus with projPDG is found in AMDB 00128 if(colN[i]==tgN && colZ[i]==tgZ) 00129 { 00130 lastI=i; 00131 lastTH =colTH[i]; // Last THreshold (A-dependent) 00132 #ifdef pdebug 00133 G4cout<<"G4QAMNCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl; 00134 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00135 #endif 00136 if(pEn<=lastTH) 00137 { 00138 #ifdef pdebug 00139 G4cout<<"G4QAMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl; 00140 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00141 #endif 00142 return 0.; // Energy is below the Threshold value 00143 } 00144 lastP =colP [i]; // Last Momentum (A-dependent) 00145 lastCS =colCS[i]; // Last CrossSect (A-dependent) 00146 if(std::fabs(lastP/pMom-1.)<tolerance) 00147 { 00148 #ifdef pdebug 00149 G4cout<<"G4QAMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; 00150 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00151 #endif 00152 return lastCS*millibarn; // Use theLastCS 00153 } 00154 in = true; // This is the case when the isotop is found in DB 00155 // Momentum pMom is in IU ! @@ Units 00156 #ifdef pdebug 00157 G4cout<<"G4QAMNCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl; 00158 #endif 00159 lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update 00160 #ifdef pdebug 00161 G4cout<<"G4QAMNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl; 00162 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00163 #endif 00164 if(lastCS<=0. && pEn>lastTH) // Correct the threshold 00165 { 00166 #ifdef pdebug 00167 G4cout<<"G4QAMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00168 #endif 00169 lastTH=pEn; 00170 } 00171 break; // Go out of the LOOP 00172 } 00173 #ifdef pdebug 00174 G4cout<<"---G4QAMNCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i] 00175 <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl; 00176 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00177 #endif 00178 j++; // Increment a#0f records found in DB for this pPDG 00179 } 00180 if(!in) // This nucleus has not been calculated previously 00181 { 00182 #ifdef pdebug 00183 G4cout<<"G4QAMNCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<lastI<<G4endl; 00184 #endif 00186 lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create 00187 if(lastCS<=0.) 00188 { 00189 lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last 00190 #ifdef pdebug 00191 G4cout<<"G4QAMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl; 00192 #endif 00193 if(pEn>lastTH) 00194 { 00195 #ifdef pdebug 00196 G4cout<<"G4QAMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00197 #endif 00198 lastTH=pEn; 00199 } 00200 } 00201 #ifdef pdebug 00202 G4cout<<"G4QAMNCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl; 00203 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00204 #endif 00205 colN.push_back(tgN); 00206 colZ.push_back(tgZ); 00207 colPDG.push_back(pPDG); 00208 colP.push_back(pMom); 00209 colTH.push_back(lastTH); 00210 colCS.push_back(lastCS); 00211 #ifdef pdebug 00212 G4cout<<"G4QAMNCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl; 00213 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00214 #endif 00215 return lastCS*millibarn; 00216 } // End of creation of the new set of parameters 00217 else 00218 { 00219 #ifdef pdebug 00220 G4cout<<"G4QAMNCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl; 00221 #endif 00222 colP[lastI]=pMom; 00223 colPDG[lastI]=pPDG; 00224 colCS[lastI]=lastCS; 00225 } 00226 } // End of parameters udate 00227 else if(pEn<=lastTH) 00228 { 00229 #ifdef pdebug 00230 G4cout<<"G4QAMNCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl; 00231 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00232 #endif 00233 return 0.; // Momentum is below the Threshold Value -> CS=0 00234 } 00235 else if(std::fabs(lastP/pMom-1.)<tolerance) 00236 { 00237 #ifdef pdebug 00238 G4cout<<"G4QAMNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; 00239 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00240 #endif 00241 return lastCS*millibarn; // Use theLastCS 00242 } 00243 else 00244 { 00245 #ifdef pdebug 00246 G4cout<<"G4QAMNCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl; 00247 #endif 00248 lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB 00249 lastP=pMom; 00250 } 00251 #ifdef pdebug 00252 G4cout<<"G4QAMNCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00253 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00254 #endif 00255 return lastCS*millibarn; 00256 }
Reimplemented from G4VQCrossSection.
Definition at line 725 of file G4QANuANuNuclearCrossSection.cc.
00726 { 00727 G4double f=Q2/4.62; 00728 G4double ff=f*f; 00729 G4double r=ff*ff; 00730 G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3))); 00731 //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par) 00732 return 1.-s_value*(1.-s_value/2); 00733 }
G4int G4QANuANuNuclearCrossSection::GetExchangePDGCode | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 742 of file G4QANuANuNuclearCrossSection.cc.
G4double G4QANuANuNuclearCrossSection::GetLastQELCS | ( | ) | [inline, virtual] |
G4double G4QANuANuNuclearCrossSection::GetLastTOTCS | ( | ) | [inline, virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 736 of file G4QANuANuNuclearCrossSection.cc.
00737 { 00738 return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space 00739 }
G4double G4QANuANuNuclearCrossSection::GetNQE_ExchangeQ2 | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 501 of file G4QANuANuNuclearCrossSection.cc.
References G4UniformRand.
00502 { 00503 static const double mpi=.13957018; // charged pi meson mass in GeV 00504 static const double MN=.931494043; // Nucleon mass (inside nucleus,atomicMassUnit,GeV) 00505 static const double dMN=MN+MN; // 2*M_N in GeV 00506 static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic 00507 static const G4int power=7; // direct power for the magic variable 00508 static const G4double pconv=1./power; // conversion power for the magic variable 00509 static const G4int nX=21; // #Of point in the Xl table (in GeV^2) 00510 static const G4int lX=nX-1; // index of the last in the Xl table 00511 static const G4int bX=lX-1; // @@ index of the before last in the Xl table 00512 static const G4int nE=20; // #Of point in the El table (in GeV^2) 00513 static const G4int bE=nE-1; // index of the last in the El table 00514 static const G4int pE=bE-1; // index of the before last in the El table 00515 // Reversed table 00516 static const G4double X0[nX]={5.21412e-05, 00517 .437860, .681908, .891529, 1.08434, 1.26751, 1.44494, 1.61915, 1.79198, 1.96493, 2.13937, 00518 2.31664, 2.49816, 2.68559, 2.88097, 3.08705, 3.30774, 3.54917, 3.82233, 4.15131, 4.62182}; 00519 static const G4double X1[nX]={.00102591, 00520 1.00443, 1.55828, 2.03126, 2.46406, 2.87311, 3.26723, 3.65199, 4.03134, 4.40835, 4.78561, 00521 5.16549, 5.55031, 5.94252, 6.34484, 6.76049, 7.19349, 7.64917, 8.13502, 8.66246, 9.25086}; 00522 static const G4double X2[nX]={.0120304, 00523 2.59903, 3.98637, 5.15131, 6.20159, 7.18024, 8.10986, 9.00426, 9.87265, 10.7217, 11.5564, 00524 12.3808, 13.1983, 14.0116, 14.8234, 15.6359, 16.4515, 17.2723, 18.1006, 18.9386, 19.7892}; 00525 static const G4double X3[nX]={.060124, 00526 5.73857, 8.62595, 10.9849, 13.0644, 14.9636, 16.7340, 18.4066, 20.0019, 21.5342, 23.0142, 00527 24.4497, 25.8471, 27.2114, 28.5467, 29.8564, 31.1434, 32.4102, 33.6589, 34.8912, 36.1095}; 00528 static const G4double X4[nX]={.0992363, 00529 8.23746, 12.1036, 15.1740, 17.8231, 20.1992, 22.3792, 24.4092, 26.3198, 28.1320, 29.8615, 00530 31.5200, 33.1169, 34.6594, 36.1536, 37.6044, 39.0160, 40.3920, 41.7353, 43.0485, 44.3354}; 00531 static const G4double X5[nX]={.0561127, 00532 7.33661, 10.5694, 13.0778, 15.2061, 17.0893, 18.7973, 20.3717, 21.8400, 23.2211, 24.5291, 00533 25.7745, 26.9655, 28.1087, 29.2094, 30.2721, 31.3003, 32.2972, 33.2656, 34.2076, 35.1265}; 00534 static const G4double X6[nX]={.0145859, 00535 4.81774, 6.83565, 8.37399, 9.66291, 10.7920, 11.8075, 12.7366, 13.5975, 14.4025, 15.1608, 00536 15.8791, 16.5628, 17.2162, 17.8427, 18.4451, 19.0259, 19.5869, 20.1300, 20.6566, 21.1706}; 00537 static const G4double X7[nX]={.00241155, 00538 2.87095, 4.02492, 4.89243, 5.61207, 6.23747, 6.79613, 7.30433, 7.77270, 8.20858, 8.61732, 00539 9.00296, 9.36863, 9.71682, 10.0495, 10.3684, 10.6749, 10.9701, 11.2550, 11.5306, 11.7982}; 00540 static const G4double X8[nX]={.000316863, 00541 1.76189, 2.44632, 2.95477, 3.37292, 3.73378, 4.05420, 4.34415, 4.61009, 4.85651, 5.08666, 00542 5.30299, 5.50738, 5.70134, 5.88609, 6.06262, 6.23178, 6.39425, 6.55065, 6.70149, 6.84742}; 00543 static const G4double X9[nX]={3.73544e-05, 00544 1.17106, 1.61289, 1.93763, 2.20259, 2.42976, 2.63034, 2.81094, 2.97582, 3.12796, 3.26949, 00545 3.40202, 3.52680, 3.64482, 3.75687, 3.86360, 3.96557, 4.06323, 4.15697, 4.24713, 4.33413}; 00546 static const G4double XA[nX]={4.19131e-06, 00547 .849573, 1.16208, 1.38955, 1.57379, 1.73079, 1.86867, 1.99221, 2.10451, 2.20770, 2.30332, 00548 2.39252, 2.47622, 2.55511, 2.62977, 2.70066, 2.76818, 2.83265, 2.89437, 2.95355, 3.01051}; 00549 static const G4double XB[nX]={4.59981e-07, 00550 .666131, .905836, 1.07880, 1.21796, 1.33587, 1.43890, 1.53080, 1.61399, 1.69011, 1.76040, 00551 1.82573, 1.88682, 1.94421, 1.99834, 2.04959, 2.09824, 2.14457, 2.18878, 2.23107, 2.27162}; 00552 static const G4double XC[nX]={4.99861e-08, 00553 .556280, .752730, .893387, 1.00587, 1.10070, 1.18317, 1.25643, 1.32247, 1.38269, 1.43809, 00554 1.48941, 1.53724, 1.58203, 1.62416, 1.66391, 1.70155, 1.73728, 1.77128, 1.80371, 1.83473}; 00555 static const G4double XD[nX]={5.40832e-09, 00556 .488069, .657650, .778236, .874148, .954621, 1.02432, 1.08599, 1.14138, 1.19172, 1.23787, 00557 1.28049, 1.32008, 1.35705, 1.39172, 1.42434, 1.45514, 1.48429, 1.51197, 1.53829, 1.56339}; 00558 static const G4double XE[nX]={5.84029e-10, 00559 .445057, .597434, .705099, .790298, .861468, .922865, .976982, 1.02542, 1.06930, 1.10939, 00560 1.14630, 1.18050, 1.21233, 1.24208, 1.27001, 1.29630, 1.32113, 1.34462, 1.36691, 1.38812}; 00561 static const G4double XF[nX]={6.30137e-11, 00562 .418735, .560003, .659168, .737230, .802138, .857898, .906854, .950515, .989915, 1.02580, 00563 1.05873, 1.08913, 1.11734, 1.14364, 1.16824, 1.19133, 1.21306, 1.23358, 1.25298, 1.27139}; 00564 static const G4double XG[nX]={6.79627e-12, 00565 .405286, .539651, .633227, .706417, .766929, .818642, .863824, .903931, .939963, .972639, 00566 1.00250, 1.02995, 1.05532, 1.07887, 1.10082, 1.12134, 1.14058, 1.15867, 1.17572, 1.19183}; 00567 static const G4double XH[nX]={7.32882e-13, 00568 .404391, .535199, .625259, .695036, .752243, .800752, .842823, .879906, .912994, .942802, 00569 .969862, .994583, 1.01729, 1.03823, 1.05763, 1.07566, 1.09246, 1.10816, 1.12286, 1.13667}; 00570 static const G4double XI[nX]={7.90251e-14, 00571 .418084, .548382, .636489, .703728, .758106, .803630, .842633, .876608, .906576, .933269, 00572 .957233, .978886, .998556, 1.01651, 1.03295, 1.04807, 1.06201, 1.07489, 1.08683, 1.09792}; 00573 static const G4double XJ[nX]={8.52083e-15, 00574 .447299, .579635, .666780, .731788, .783268, .825512, .861013, .891356, .917626, .940597, 00575 .960842, .978802, .994820, 1.00917, 1.02208, 1.03373, 1.04427, 1.05383, 1.06253, 1.07046}; 00576 // Direct table 00577 static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0], 00578 X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]}; 00579 static const G4double dX[nE]={ 00580 (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX, 00581 (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX, 00582 (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX, 00583 (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX, 00584 (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX}; 00585 static const G4double* Xl[nE]= 00586 {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ}; 00587 static const G4double I0[nX]={0, 00588 .354631, 1.08972, 2.05138, 3.16564, 4.38343, 5.66828, 6.99127, 8.32858, 9.65998, 10.9680, 00589 12.2371, 13.4536, 14.6050, 15.6802, 16.6686, 17.5609, 18.3482, 19.0221, 19.5752, 20.0000}; 00590 static const G4double I1[nX]={0, 00591 .281625, .877354, 1.67084, 2.60566, 3.64420, 4.75838, 5.92589, 7.12829, 8.34989, 9.57708, 00592 10.7978, 12.0014, 13.1781, 14.3190, 15.4162, 16.4620, 17.4496, 18.3724, 19.2245, 20.0000}; 00593 static const G4double I2[nX]={0, 00594 .201909, .642991, 1.24946, 1.98463, 2.82370, 3.74802, 4.74263, 5.79509, 6.89474, 8.03228, 00595 9.19947, 10.3889, 11.5938, 12.8082, 14.0262, 15.2427, 16.4527, 17.6518, 18.8356, 20.0000}; 00596 static const G4double I3[nX]={0, 00597 .140937, .461189, .920216, 1.49706, 2.17728, 2.94985, 3.80580, 4.73758, 5.73867, 6.80331, 00598 7.92637, 9.10316, 10.3294, 11.6013, 12.9150, 14.2672, 15.6548, 17.0746, 18.5239, 20.0000}; 00599 static const G4double I4[nX]={0, 00600 .099161, .337358, .694560, 1.16037, 1.72761, 2.39078, 3.14540, 3.98768, 4.91433, 5.92245, 00601 7.00942, 8.17287, 9.41060, 10.7206, 12.1010, 13.5500, 15.0659, 16.6472, 18.2924, 20.0000}; 00602 static const G4double I5[nX]={0, 00603 .071131, .255084, .543312, .932025, 1.41892, 2.00243, 2.68144, 3.45512, 4.32283, 5.28411, 00604 6.33859, 7.48602, 8.72621, 10.0590, 11.4844, 13.0023, 14.6128, 16.3158, 18.1115, 20.0000}; 00605 static const G4double I6[nX]={0, 00606 .053692, .202354, .443946, .778765, 1.20774, 1.73208, 2.35319, 3.07256, 3.89177, 4.81249, 00607 5.83641, 6.96528, 8.20092, 9.54516, 10.9999, 12.5670, 14.2486, 16.0466, 17.9630, 20.0000}; 00608 static const G4double I7[nX]={0, 00609 .043065, .168099, .376879, .672273, 1.05738, 1.53543, 2.10973, 2.78364, 3.56065, 4.44429, 00610 5.43819, 6.54610, 7.77186, 9.11940, 10.5928, 12.1963, 13.9342, 15.8110, 17.8313, 20.0000}; 00611 static const G4double I8[nX]={0, 00612 .036051, .143997, .327877, .592202, .941572, 1.38068, 1.91433, 2.54746, 3.28517, 4.13277, 00613 5.09574, 6.17984, 7.39106, 8.73568, 10.2203, 11.8519, 13.6377, 15.5854, 17.7033, 20.0000}; 00614 static const G4double I9[nX]={0, 00615 .030977, .125727, .289605, .528146, .846967, 1.25183, 1.74871, 2.34384, 3.04376, 3.85535, 00616 4.78594, 5.84329, 7.03567, 8.37194, 9.86163, 11.5150, 13.3430, 15.3576, 17.5719, 20.0000}; 00617 static const G4double IA[nX]={0, 00618 .027129, .111420, .258935, .475812, .768320, 1.14297, 1.60661, 2.16648, 2.83034, 3.60650, 00619 4.50394, 5.53238, 6.70244, 8.02569, 9.51488, 11.1841, 13.0488, 15.1264, 17.4362, 20.0000}; 00620 static const G4double IB[nX]={0, 00621 .024170, .100153, .234345, .433198, .703363, 1.05184, 1.48607, 2.01409, 2.64459, 3.38708, 00622 4.25198, 5.25084, 6.39647, 7.70319, 9.18708, 10.8663, 12.7617, 14.8968, 17.2990, 20.0000}; 00623 static const G4double IC[nX]={0, 00624 .021877, .091263, .214670, .398677, .650133, .976322, 1.38510, 1.88504, 2.48555, 3.19709, 00625 4.03129, 5.00127, 6.12184, 7.40989, 8.88482, 10.5690, 12.4888, 14.6748, 17.1638, 20.0000}; 00626 static const G4double ID[nX]={0, 00627 .020062, .084127, .198702, .370384, .606100, .913288, 1.30006, 1.77535, 2.34912, 3.03253, 00628 3.83822, 4.78063, 5.87634, 7.14459, 8.60791, 10.2929, 12.2315, 14.4621, 17.0320, 20.0000}; 00629 static const G4double IE[nX]={0, 00630 .018547, .078104, .185102, .346090, .567998, .858331, 1.22535, 1.67824, 2.22735, 2.88443, 00631 3.66294, 4.57845, 5.64911, 6.89637, 8.34578, 10.0282, 11.9812, 14.2519, 16.8993, 20.0000}; 00632 static const G4double IF[nX]={0, 00633 .017143, .072466, .172271, .323007, .531545, .805393, 1.15288, 1.58338, 2.10754, 2.73758, 00634 3.48769, 4.37450, 5.41770, 6.64092, 8.07288, 9.74894, 11.7135, 14.0232, 16.7522, 20.0000}; 00635 static const G4double IG[nX]={0, 00636 .015618, .066285, .158094, .297316, .490692, .745653, 1.07053, 1.47479, 1.96931, 2.56677, 00637 3.28205, 4.13289, 5.14068, 6.33158, 7.73808, 9.40133, 11.3745, 13.7279, 16.5577, 20.0000}; 00638 static const G4double IH[nX]={0, 00639 .013702, .058434, .139923, .264115, .437466, .667179, .961433, 1.32965, 1.78283, 2.33399, 00640 2.99871, 3.79596, 4.74916, 5.88771, 7.24937, 8.88367, 10.8576, 13.2646, 16.2417, 20.0000}; 00641 static const G4double II[nX]={0, 00642 .011264, .048311, .116235, .220381, .366634, .561656, .813132, 1.13008, 1.52322, 2.00554, 00643 2.59296, 3.30542, 4.16834, 5.21490, 6.48964, 8.05434, 9.99835, 12.4580, 15.6567, 20.0000}; 00644 static const G4double IJ[nX]={0, 00645 .008628, .037206, .089928, .171242, .286114, .440251, .640343, .894382, 1.21208, 1.60544, 00646 2.08962, 2.68414, 3.41486, 4.31700, 5.44048, 6.85936, 8.69067, 11.1358, 14.5885, 20.0000}; 00647 static const G4double* Il[nE]= 00648 {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ}; 00649 static const G4double lE[nE]={ 00650 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292, 00651 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218}; 00652 static const G4double lEmi=lE[0]; 00653 static const G4double lEma=lE[nE-1]; 00654 static const G4double dlE=(lEma-lEmi)/bE; 00655 //*************************************************************************************** 00656 G4double Enu=lastE; // Get energy of the last calculated cross-section 00657 G4double lEn=std::log(Enu); // log(E) for interpolation 00658 G4double rE=(lEn-lEmi)/dlE; // Position of the energy 00659 G4int fE=static_cast<int>(rE); // Left bin for interpolation 00660 if(fE<0) fE=0; 00661 if(fE>pE)fE=pE; 00662 G4int sE=fE+1; // Right bin for interpolation 00663 G4double dE=rE-fE; // relative log shift from the left bin 00664 G4double dEnu=Enu+Enu; // doubled energy of nu/anu 00665 G4double Enu2=Enu*Enu; // squared energy of nu/anu 00666 G4double Emu=Enu; // Free Energy of neutrino/anti-neutrino 00667 G4double ME=Enu*MN; // M*E 00668 G4double dME=ME+ME; // 2*M*E 00669 G4double dEMN=(dEnu+MN)*ME; 00670 G4double sqE=Enu*ME; 00671 G4double E2M=MN*Enu2; 00672 G4double ymax=(E2M+sqE)/dEMN; 00673 G4double Q2mi=0.; // Q2_min(E_nu) 00674 G4double Q2ma=dME*ymax; // Q2_max(E_nu) 00675 G4double Q2nq=Emu*dMN-mcV; 00676 if(Q2ma>Q2nq) Q2ma=Q2nq; // Correction for Non Quasi Elastic 00677 // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma --- 00678 G4double Rmi=Q2mi/Q2ma; 00679 G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78); 00680 // --- E-interpolation must be done in a log scale --- 00681 G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu) 00682 G4double Xma=std::pow((shift-1.),power); // X_max(E_nu) 00683 // Find the integral values integ(Xmi) & integ(Xma) using the direct table 00684 G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step 00685 G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum 00686 G4double rXi=(Xmi-iXmi)/idX; 00687 G4int iXi=static_cast<int>(rXi); 00688 if(iXi<0) iXi=0; 00689 if(iXi>bX) iXi=bX; 00690 G4double dXi=rXi-iXi; 00691 G4double bntil=Il[fE][iXi]; 00692 G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil); 00693 G4double bntir=Il[sE][iXi]; 00694 G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir); 00695 G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral 00696 // 00697 G4double rXa=(Xma-iXmi)/idX; 00698 G4int iXa=static_cast<int>(rXa); 00699 if(iXa<0) iXa=0; 00700 if(iXa>bX) iXa=bX; 00701 G4double dXa=rXa-iXa; 00702 G4double bntal=Il[fE][iXa]; 00703 G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal); 00704 G4double bntar=Il[sE][iXa]; 00705 G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar); 00706 G4double inta=intal+dE*(intar-intal);// interpolated end of the integral 00707 // 00708 // *** Find X using the reversed table *** 00709 G4double intx=inti+(inta-inti)*G4UniformRand(); 00710 G4int intc=static_cast<int>(intx); 00711 if(intc<0) intc=0; 00712 if(intc>bX) intc=bX; 00713 G4double dint=intx-intc; 00714 G4double mXl=Xl[fE][intc]; 00715 G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl); 00716 G4double mXr=Xl[sE][intc]; 00717 G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr); 00718 G4double X=Xlb+dE*(Xrb-Xlb); // interpolated X value 00719 G4double R=shift-std::pow(X,pconv); 00720 G4double Q2=R*Q2ma; 00721 return Q2*GeV*GeV; 00722 }
G4VQCrossSection * G4QANuANuNuclearCrossSection::GetPointer | ( | ) | [static] |
Definition at line 72 of file G4QANuANuNuclearCrossSection.cc.
Referenced by G4QInelastic::GetMeanFreePath(), and G4QInelastic::PostStepDoIt().
00073 { 00074 static G4QANuANuNuclearCrossSection theCrossSection; //*Static body of the Cross Section* 00075 return &theCrossSection; 00076 }
G4double G4QANuANuNuclearCrossSection::GetQEL_ExchangeQ2 | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 424 of file G4QANuANuNuclearCrossSection.cc.
References G4UniformRand.
00425 { 00426 static const double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV) 00427 static const G4double power=-3.5; // direct power for the magic variable 00428 static const G4double pconv=1./power;// conversion power for the magic variable 00429 static const G4int nQ2=101; // #Of point in the Q2l table (in GeV^2) 00430 static const G4int lQ2=nQ2-1; // index of the last in the Q2l table 00431 static const G4int bQ2=lQ2-1; // index of the before last in the Q2 ltable 00432 // Reversed table 00433 static const G4double Xl[nQ2]={5.20224e-16, 00434 .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898, 00435 .093878, .103332, .112844, .122410, .132023, .141680, .151376, .161109, .170875, .180672, 00436 .190499, .200352, .210230, .220131, .230055, .239999, .249963, .259945, .269944, .279960, 00437 .289992, .300039, .310099, .320173, .330260, .340359, .350470, .360592, .370724, .380867, 00438 .391019, .401181, .411352, .421531, .431719, .441915, .452118, .462329, .472547, .482771, 00439 .493003, .503240, .513484, .523734, .533989, .544250, .554517, .564788, .575065, .585346, 00440 .595632, .605923, .616218, .626517, .636820, .647127, .657438, .667753, .678072, .688394, 00441 .698719, .709048, .719380, .729715, .740053, .750394, .760738, .771085, .781434, .791786, 00442 .802140, .812497, .822857, .833219, .843582, .853949, .864317, .874687, .885060, .895434, 00443 .905810, .916188, .926568, .936950, .947333, .957719, .968105, .978493, .988883, .999275}; 00444 // Direct table 00445 static const G4double Xmax=Xl[lQ2]; 00446 static const G4double Xmin=Xl[0]; 00447 static const G4double dX=(Xmax-Xmin)/lQ2; // step in X(Q2, GeV^2) 00448 static const G4double inl[nQ2]={0, 00449 1.52225, 2.77846, 3.96651, 5.11612, 6.23990, 7.34467, 8.43466, 9.51272, 10.5809, 11.6406, 00450 12.6932, 13.7394, 14.7801, 15.8158, 16.8471, 17.8743, 18.8979, 19.9181, 20.9353, 21.9496, 00451 22.9614, 23.9707, 24.9777, 25.9826, 26.9855, 27.9866, 28.9860, 29.9837, 30.9798, 31.9745, 00452 32.9678, 33.9598, 34.9505, 35.9400, 36.9284, 37.9158, 38.9021, 39.8874, 40.8718, 41.8553, 00453 42.8379, 43.8197, 44.8007, 45.7810, 46.7605, 47.7393, 48.7174, 49.6950, 50.6718, 51.6481, 00454 52.6238, 53.5990, 54.5736, 55.5476, 56.5212, 57.4943, 58.4670, 59.4391, 60.4109, 61.3822, 00455 62.3531, 63.3236, 64.2937, 65.2635, 66.2329, 67.2019, 68.1707, 69.1390, 70.1071, 71.0748, 00456 72.0423, 73.0095, 73.9763, 74.9429, 75.9093, 76.8754, 77.8412, 78.8068, 79.7721, 80.7373, 00457 81.7022, 82.6668, 83.6313, 84.5956, 85.5596, 86.5235, 87.4872, 88.4507, 89.4140, 90.3771, 00458 91.3401, 92.3029, 93.2656, 94.2281, 95.1904, 96.1526, 97.1147, 98.0766, 99.0384, 100.000}; 00459 G4double Enu=lastE; // Get energy of the last calculated cross-section 00460 G4double dEnu=Enu+Enu; // doubled energy of nu/anu 00461 G4double Enu2=Enu*Enu; // squared energy of nu/anu 00462 G4double ME=Enu*MN; // M*E 00463 G4double dME=ME+ME; // 2*M*E 00464 G4double dEMN=(dEnu+MN)*ME; 00465 G4double sqE=Enu*ME; 00466 G4double E2M=MN*Enu2; 00467 G4double ymax=(E2M+sqE)/dEMN; 00468 G4double Q2mi=0.; // Q2_min(E_nu) 00469 G4double Q2ma=dME*ymax; // Q2_max(E_nu) 00470 G4double Xma=std::pow((1.+Q2mi),power); // X_max(E_nu) 00471 G4double Xmi=std::pow((1.+Q2ma),power); // X_min(E_nu) 00472 // Find the integral values integ(Xmi) & integ(Xma) using the direct table 00473 G4double rXi=(Xmi-Xmin)/dX; 00474 G4int iXi=static_cast<int>(rXi); 00475 if(iXi<0) iXi=0; 00476 if(iXi>bQ2) iXi=bQ2; 00477 G4double dXi=rXi-iXi; 00478 G4double bnti=inl[iXi]; 00479 G4double inti=bnti+dXi*(inl[iXi+1]-bnti); 00480 // 00481 G4double rXa=(Xma-Xmin)/dX; 00482 G4int iXa=static_cast<int>(rXa); 00483 if(iXa<0) iXa=0; 00484 if(iXa>bQ2) iXa=bQ2; 00485 G4double dXa=rXa-iXa; 00486 G4double bnta=inl[iXa]; 00487 G4double inta=bnta+dXa*(inl[iXa+1]-bnta); 00488 // *** Find X using the reversed table *** 00489 G4double intx=inti+(inta-inti)*G4UniformRand(); 00490 G4int intc=static_cast<int>(intx); 00491 if(intc<0) intc=0; 00492 if(intc>bQ2) intc=bQ2; // If it is more than max, then the BAD extrapolation 00493 G4double dint=intx-intc; 00494 G4double mX=Xl[intc]; 00495 G4double X=mX+dint*(Xl[intc+1]-mX); 00496 G4double Q2=std::pow(X,pconv)-1.; 00497 return Q2*GeV*GeV; 00498 }
G4double G4QANuANuNuclearCrossSection::ThresholdEnergy | ( | G4int | Z, | |
G4int | N, | |||
G4int | PDG = -14 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 259 of file G4QANuANuNuclearCrossSection.cc.
Referenced by GetCrossSection().