#include <G4QNuENuclearCrossSection.hh>
Inheritance diagram for G4QNuENuclearCrossSection:
Public Member Functions | |
~G4QNuENuclearCrossSection () | |
G4double | ThresholdEnergy (G4int Z, G4int N, G4int PDG=12) |
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 | |
G4QNuENuclearCrossSection () |
Definition at line 49 of file G4QNuENuclearCrossSection.hh.
G4QNuENuclearCrossSection::G4QNuENuclearCrossSection | ( | ) | [inline, protected] |
G4QNuENuclearCrossSection::~G4QNuENuclearCrossSection | ( | ) |
G4double G4QNuENuclearCrossSection::CalculateCrossSection | ( | G4bool | CS, | |
G4int | F, | |||
G4int | I, | |||
G4int | PDG, | |||
G4int | Z, | |||
G4int | N, | |||
G4double | Momentum | |||
) | [virtual] |
Implements G4VQCrossSection.
Definition at line 278 of file G4QNuENuclearCrossSection.cc.
References G4cerr, G4cout, and G4endl.
Referenced by GetCrossSection().
00280 { 00281 static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2 00282 static const G4int nE=65; // !! If change this, change it in GetFunctions()33/65!! 00283 static const G4int mL=nE-1; 00284 static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV) 00285 static const G4double dmN=mN+mN; // Doubled nucleon mass (2*AtomicMassUnit, GeV) 00286 static const G4double me=.00051099892;// electron mass in GeV 00287 static const G4double me2=me*me; // Squared mass of an electron in GeV^2 00288 static const G4double EMi=me+me2/dmN; // Universal threshold of the reaction in GeV 00289 static const G4double EMa=300.; // Maximum tabulated Energy of nu_e in GeV 00290 // *** Begin of the Associative memory for acceleration of the cross section calculations 00291 static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional) 00292 static G4bool first=true; // Flag of initialization of the energy axis 00293 // *** End of Static Definitions (Associative Memory) *** 00294 //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Electron 00295 //G4double TotEnergy2=Momentum; 00296 onlyCS=CS; // Flag to calculate only CS (not TX & QE) 00297 lastE=Momentum/GeV; // Kinetic energy of the electron neutrino (in GeV) 00298 if (lastE<=EMi) // Energy is below the minimum energy in the table 00299 { 00300 lastE=0.; 00301 lastSig=0.; 00302 return 0.; 00303 } 00304 G4int Z=targZ; // New Z, which can change the sign 00305 if(F<=0) // This isotope was not the last used isotop 00306 { 00307 if(F<0) // This isotope was found in DAMDB =-------=> RETRIEVE 00308 { 00309 lastTX =(*TX)[I]; // Pointer to the prepared TX function (same isotope) 00310 lastQE =(*QE)[I]; // Pointer to the prepared QE function (same isotope) 00311 } 00312 else // This isotope wasn't calculated previously => CREATE 00313 { 00314 if(first) 00315 { 00316 lastEN = new G4double[nE]; // This must be done only once! 00317 Z=-Z; // To explain GetFunctions that E-axis must be filled 00318 first=false; // To make it only once 00319 } 00320 lastTX = new G4double[nE]; // Allocate memory for the new TX function 00321 lastQE = new G4double[nE]; // Allocate memory for the new QE function 00322 G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK) 00323 if(res<0) G4cerr<<"*W*G4NuENuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl; 00324 // *** The synchronization check *** 00325 G4int sync=TX->size(); 00326 if(sync!=I) G4cerr<<"***G4NuENuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl; 00327 TX->push_back(lastTX); 00328 QE->push_back(lastQE); 00329 } // End of creation of the new set of parameters 00330 } // End of parameters udate 00331 // =----------------= NOW Calculate the Cross Section =-----------------= 00332 if (lastE<=EMi) // Check that the neutrinoEnergy is higher than ThreshE 00333 { 00334 lastE=0.; 00335 lastSig=0.; 00336 return 0.; 00337 } 00338 if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization 00339 { 00340 G4int chk=1; 00341 G4int ran=mL/2; 00342 G4int sep=ran; // as a result = an index of the left edge of the interval 00343 while(ran>=2) 00344 { 00345 G4int newran=ran/2; 00346 if(lastE<=lastEN[sep]) sep-=newran; 00347 else sep+=newran; 00348 ran=newran; 00349 chk=chk+chk; 00350 } 00351 if(chk+chk!=mL) G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Table! mL="<<mL<<G4endl; 00352 G4double lowE=lastEN[sep]; 00353 G4double highE=lastEN[sep+1]; 00354 G4double lowTX=lastTX[sep]; 00355 if(lastE<lowE||sep>=mL||lastE>highE) 00356 G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE 00357 <<", sep="<<sep<<", mL="<<mL<<G4endl; 00358 lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E 00359 if(!onlyCS) // Skip the differential cross-section parameters 00360 { 00361 G4double lowQE=lastQE[sep]; 00362 lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE; 00363 #ifdef pdebug 00364 G4cout<<"G4NuENuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl; 00365 #endif 00366 } 00367 } 00368 else 00369 { 00370 lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking... 00371 lastQEL=lastQE[mL]; 00372 } 00373 if(lastQEL<0.) lastQEL = 0.; 00374 if(lastSig<0.) lastSig = 0.; 00375 // The cross-sections are expected to be in mb 00376 lastSig*=mb38; 00377 if(!onlyCS) lastQEL*=mb38; 00378 return lastSig; 00379 }
G4double G4QNuENuclearCrossSection::GetCrossSection | ( | G4bool | fCS, | |
G4double | pMom, | |||
G4int | tgZ, | |||
G4int | tgN, | |||
G4int | pPDG = 0 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 90 of file G4QNuENuclearCrossSection.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<<"G4QNENCS::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!=12) 00109 { 00110 #ifdef debug 00111 G4cout<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"---G4QNENCrossSec::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<<"G4QNENCS::GetCrSec: CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<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<<"G4QNENCrossSection::GetCrossSect:NewThresh="<<lastTH<<",T="<<pEn<<G4endl; 00192 #endif 00193 if(pEn>lastTH) 00194 { 00195 #ifdef pdebug 00196 G4cout<<"G4QNENCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00197 #endif 00198 lastTH=pEn; 00199 } 00200 } 00201 #ifdef pdebug 00202 G4cout<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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 767 of file G4QNuENuclearCrossSection.cc.
00768 { 00769 G4double f=Q2/4.62; 00770 G4double ff=f*f; 00771 G4double r=ff*ff; 00772 G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3))); 00773 //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par) 00774 return 1.-s_value*(1.-s_value/2); 00775 }
G4int G4QNuENuclearCrossSection::GetExchangePDGCode | ( | ) | [virtual] |
G4double G4QNuENuclearCrossSection::GetLastQELCS | ( | ) | [inline, virtual] |
G4double G4QNuENuclearCrossSection::GetLastTOTCS | ( | ) | [inline, virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 778 of file G4QNuENuclearCrossSection.cc.
00779 { 00780 return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space 00781 }
G4double G4QNuENuclearCrossSection::GetNQE_ExchangeQ2 | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 535 of file G4QNuENuclearCrossSection.cc.
References G4UniformRand.
00536 { 00537 static const double mpi=.13957018; // charged pi meson mass in GeV 00538 static const G4double me=.00051099892;// electron mass in GeV 00539 static const G4double me2=me*me; // Squared mass of an electron in GeV^2 00540 static const G4double hme2=me2/2; // .5*m_e^2 in GeV^2 00541 static const double MN=.931494043; // Nucleon mass (inside nucleus,atomicMassUnit,GeV) 00542 static const double MN2=MN*MN; // M_N^2 in GeV^2 00543 static const double dMN=MN+MN; // 2*M_N in GeV 00544 static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic 00545 static const G4int power=7; // direct power for the magic variable 00546 static const G4double pconv=1./power; // conversion power for the magic variable 00547 static const G4int nX=21; // #Of point in the Xl table (in GeV^2) 00548 static const G4int lX=nX-1; // index of the last in the Xl table 00549 static const G4int bX=lX-1; // @@ index of the before last in the Xl table 00550 static const G4int nE=20; // #Of point in the El table (in GeV^2) 00551 static const G4int bE=nE-1; // index of the last in the El table 00552 static const G4int pE=bE-1; // index of the before last in the El table 00553 // Reversed table 00554 static const G4double X0[nX]={6.14081e-05, 00555 .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244, 00556 2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419}; 00557 static const G4double X1[nX]={.00125268, 00558 .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049, 00559 4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330}; 00560 static const G4double X2[nX]={.015694, 00561 1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728, 00562 10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276}; 00563 static const G4double X3[nX]={.0866877, 00564 4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302, 00565 21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115}; 00566 static const G4double X4[nX]={.160483, 00567 5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127, 00568 29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265}; 00569 static const G4double X5[nX]={.0999307, 00570 5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398, 00571 26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649}; 00572 static const G4double X6[nX]={.0276367, 00573 3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519, 00574 16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148}; 00575 static const G4double X7[nX]={.00472383, 00576 2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473, 00577 9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057}; 00578 static const G4double X8[nX]={.000630783, 00579 1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480, 00580 5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037}; 00581 static const G4double X9[nX]={7.49179e-05, 00582 .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823, 00583 3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028}; 00584 static const G4double XA[nX]={8.43437e-06, 00585 .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704, 00586 2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823}; 00587 static const G4double XB[nX]={9.27028e-07, 00588 .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059, 00589 1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341}; 00590 static const G4double XC[nX]={1.00807e-07, 00591 .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668, 00592 1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014}; 00593 static const G4double XD[nX]={1.09102e-08, 00594 .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878, 00595 1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617}; 00596 static const G4double XE[nX]={1.17831e-09, 00597 .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374, 00598 .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571}; 00599 static const G4double XF[nX]={1.27141e-10, 00600 .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046, 00601 .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338}; 00602 static const G4double XG[nX]={1.3713e-11, 00603 .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870, 00604 .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379}; 00605 static const G4double XH[nX]={1.47877e-12, 00606 .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421, 00607 .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197}; 00608 static const G4double XI[nX]={1.59454e-13, 00609 .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806, 00610 .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872}; 00611 static const G4double XJ[nX]={1.71931e-14, 00612 .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147, 00613 .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813}; 00614 // Direct table 00615 static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0], 00616 X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]}; 00617 static const G4double dX[nE]={ 00618 (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX, 00619 (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX, 00620 (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX, 00621 (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX, 00622 (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX}; 00623 static const G4double* Xl[nE]= 00624 {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ}; 00625 static const G4double I0[nX]={0, 00626 .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555, 00627 13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000}; 00628 static const G4double I1[nX]={0, 00629 .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365, 00630 13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000}; 00631 static const G4double I2[nX]={0, 00632 .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586, 00633 12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000}; 00634 static const G4double I3[nX]={0, 00635 .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461, 00636 12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000}; 00637 static const G4double I4[nX]={0, 00638 .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515, 00639 11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000}; 00640 static const G4double I5[nX]={0, 00641 .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040, 00642 10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000}; 00643 static const G4double I6[nX]={0, 00644 .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243, 00645 9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000}; 00646 static const G4double I7[nX]={0, 00647 .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717, 00648 8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000}; 00649 static const G4double I8[nX]={0, 00650 .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415, 00651 8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000}; 00652 static const G4double I9[nX]={0, 00653 .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672, 00654 8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000}; 00655 static const G4double IA[nX]={0, 00656 .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265, 00657 7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000}; 00658 static const G4double IB[nX]={0, 00659 .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298, 00660 7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000}; 00661 static const G4double IC[nX]={0, 00662 .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696, 00663 7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000}; 00664 static const G4double ID[nX]={0, 00665 .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212, 00666 7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000}; 00667 static const G4double IE[nX]={0, 00668 .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949, 00669 7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000}; 00670 static const G4double IF[nX]={0, 00671 .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164, 00672 7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000}; 00673 static const G4double IG[nX]={0, 00674 .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249, 00675 7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000}; 00676 static const G4double IH[nX]={0, 00677 .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804, 00678 7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000}; 00679 static const G4double II[nX]={0, 00680 .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722, 00681 7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000}; 00682 static const G4double IJ[nX]={0, 00683 .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196, 00684 6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000}; 00685 static const G4double* Il[nE]= 00686 {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ}; 00687 static const G4double lE[nE]={ 00688 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292, 00689 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218}; 00690 static const G4double lEmi=lE[0]; 00691 static const G4double lEma=lE[nE-1]; 00692 static const G4double dlE=(lEma-lEmi)/bE; 00693 //*************************************************************************************** 00694 G4double Enu=lastE; // Get energy of the last calculated cross-section 00695 G4double lEn=std::log(Enu); // log(E) for interpolation 00696 G4double rE=(lEn-lEmi)/dlE; // Position of the energy 00697 G4int fE=static_cast<int>(rE); // Left bin for interpolation 00698 if(fE<0) fE=0; 00699 if(fE>pE)fE=pE; 00700 G4int sE=fE+1; // Right bin for interpolation 00701 G4double dE=rE-fE; // relative log shift from the left bin 00702 G4double dEnu=Enu+Enu; // doubled energy of nu/anu 00703 G4double Enu2=Enu*Enu; // squared energy of nu/anu 00704 G4double Ee=Enu-me; // Free Energy of neutrino/anti-neutrino 00705 G4double ME=Enu*MN; // M*E 00706 G4double dME=ME+ME; // 2*M*E 00707 G4double dEMN=(dEnu+MN)*ME; 00708 G4double MEm=ME-hme2; 00709 G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2); 00710 G4double E2M=MN*Enu2-(Enu+MN)*hme2; 00711 G4double ymax=(E2M+sqE)/dEMN; 00712 G4double ymin=(E2M-sqE)/dEMN; 00713 G4double rmin=1.-ymin; 00714 G4double rhm2E=hme2/Enu2; 00715 G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu) 00716 G4double Q2ma=dME*ymax; // Q2_max(E_nu) 00717 G4double Q2nq=Ee*dMN-mcV; 00718 if(Q2ma>Q2nq) Q2ma=Q2nq; // Correction for Non Quasi Elastic 00719 // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma --- 00720 G4double Rmi=Q2mi/Q2ma; 00721 G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu 00722 // --- E-interpolation must be done in a log scale --- 00723 G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu) 00724 G4double Xma=std::pow((shift-1.),power); // X_max(E_nu) 00725 // Find the integral values integ(Xmi) & integ(Xma) using the direct table 00726 G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step 00727 G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum 00728 G4double rXi=(Xmi-iXmi)/idX; 00729 G4int iXi=static_cast<int>(rXi); 00730 if(iXi<0) iXi=0; 00731 if(iXi>bX) iXi=bX; 00732 G4double dXi=rXi-iXi; 00733 G4double bntil=Il[fE][iXi]; 00734 G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil); 00735 G4double bntir=Il[sE][iXi]; 00736 G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir); 00737 G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral 00738 // 00739 G4double rXa=(Xma-iXmi)/idX; 00740 G4int iXa=static_cast<int>(rXa); 00741 if(iXa<0) iXa=0; 00742 if(iXa>bX) iXa=bX; 00743 G4double dXa=rXa-iXa; 00744 G4double bntal=Il[fE][iXa]; 00745 G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal); 00746 G4double bntar=Il[sE][iXa]; 00747 G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar); 00748 G4double inta=intal+dE*(intar-intal);// interpolated end of the integral 00749 // 00750 // *** Find X using the reversed table *** 00751 G4double intx=inti+(inta-inti)*G4UniformRand(); 00752 G4int intc=static_cast<int>(intx); 00753 if(intc<0) intc=0; 00754 if(intc>bX) intc=bX; 00755 G4double dint=intx-intc; 00756 G4double mXl=Xl[fE][intc]; 00757 G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl); 00758 G4double mXr=Xl[sE][intc]; 00759 G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr); 00760 G4double X=Xlb+dE*(Xrb-Xlb); // interpolated X value 00761 G4double R=shift-std::pow(X,pconv); 00762 G4double Q2=R*Q2ma; 00763 return Q2*GeV*GeV; 00764 }
G4VQCrossSection * G4QNuENuclearCrossSection::GetPointer | ( | ) | [static] |
Definition at line 72 of file G4QNuENuclearCrossSection.cc.
Referenced by G4QInelastic::GetMeanFreePath(), and G4QInelastic::PostStepDoIt().
00073 { 00074 static G4QNuENuclearCrossSection theCrossSection; //**Static body of the Cross Section** 00075 return &theCrossSection; 00076 }
G4double G4QNuENuclearCrossSection::GetQEL_ExchangeQ2 | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 450 of file G4QNuENuclearCrossSection.cc.
References G4UniformRand.
00451 { 00452 static const G4double me=.00051099892; // electron mass in GeV 00453 static const G4double me2=me*me; // Squared mass of an electron in GeV^2 00454 static const G4double hme2=me2/2; // .5*m_e^2 in GeV^2 00455 static const G4double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV) 00456 static const double MN2=MN*MN; // M_N^2 in GeV^2 00457 static const G4double power=-3.5; // direct power for the magic variable 00458 static const G4double pconv=1./power;// conversion power for the magic variable 00459 static const G4int nQ2=101; // #Of point in the Q2l table (in GeV^2) 00460 static const G4int lQ2=nQ2-1; // index of the last in the Q2l table 00461 static const G4int bQ2=lQ2-1; // index of the before last in the Q2 ltable 00462 // Reversed table 00463 static const G4double Xl[nQ2]={1.87905e-10, 00464 .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374, 00465 .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069, 00466 .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718, 00467 .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169, 00468 .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670, 00469 .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720, 00470 .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457, 00471 .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366, 00472 .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128, 00473 .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539}; 00474 // Direct table 00475 static const G4double Xmax=Xl[lQ2]; 00476 static const G4double Xmin=Xl[0]; 00477 static const G4double dX=(Xmax-Xmin)/lQ2; // step in X(Q2, GeV^2) 00478 static const G4double inl[nQ2]={0, 00479 1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811, 00480 15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039, 00481 26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876, 00482 36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950, 00483 46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377, 00484 55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162, 00485 65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852, 00486 74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767, 00487 83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106, 00488 92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000}; 00489 G4double Enu=lastE; // Get energy of the last calculated cross-section 00490 G4double dEnu=Enu+Enu; // doubled energy of nu/anu 00491 G4double Enu2=Enu*Enu; // squared energy of nu/anu 00492 G4double ME=Enu*MN; // M*E 00493 G4double dME=ME+ME; // 2*M*E 00494 G4double dEMN=(dEnu+MN)*ME; 00495 G4double MEm=ME-hme2; 00496 G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2); 00497 G4double E2M=MN*Enu2-(Enu+MN)*hme2; 00498 G4double ymax=(E2M+sqE)/dEMN; 00499 G4double ymin=(E2M-sqE)/dEMN; 00500 G4double rmin=1.-ymin; 00501 G4double rhm2E=hme2/Enu2; 00502 G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu) 00503 G4double Q2ma=dME*ymax; // Q2_max(E_nu) 00504 G4double Xma=std::pow((1.+Q2mi),power); // X_max(E_nu) 00505 G4double Xmi=std::pow((1.+Q2ma),power); // X_min(E_nu) 00506 // Find the integral values integ(Xmi) & integ(Xma) using the direct table 00507 G4double rXi=(Xmi-Xmin)/dX; 00508 G4int iXi=static_cast<int>(rXi); 00509 if(iXi<0) iXi=0; 00510 if(iXi>bQ2) iXi=bQ2; 00511 G4double dXi=rXi-iXi; 00512 G4double bnti=inl[iXi]; 00513 G4double inti=bnti+dXi*(inl[iXi+1]-bnti); 00514 // 00515 G4double rXa=(Xma-Xmin)/dX; 00516 G4int iXa=static_cast<int>(rXa); 00517 if(iXa<0) iXa=0; 00518 if(iXa>bQ2) iXa=bQ2; 00519 G4double dXa=rXa-iXa; 00520 G4double bnta=inl[iXa]; 00521 G4double inta=bnta+dXa*(inl[iXa+1]-bnta); 00522 // *** Find X using the reversed table *** 00523 G4double intx=inti+(inta-inti)*G4UniformRand(); 00524 G4int intc=static_cast<int>(intx); 00525 if(intc<0) intc=0; 00526 if(intc>bQ2) intc=bQ2; // If it is more than max, then the BAD extrapolation 00527 G4double dint=intx-intc; 00528 G4double mX=Xl[intc]; 00529 G4double X=mX+dint*(Xl[intc+1]-mX); 00530 G4double Q2=std::pow(X,pconv)-1.; 00531 return Q2*GeV*GeV; 00532 }
Reimplemented from G4VQCrossSection.
Definition at line 259 of file G4QNuENuclearCrossSection.cc.
Referenced by GetCrossSection().
00260 { 00261 //static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/GeV; 00262 //static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/GeV; 00263 //static const G4double mDeut = G4NucleiProperties::GetNuclearMass(2,1)/GeV/2.; 00264 static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV) 00265 static const G4double dmN=mN+mN; // Doubled nucleon mass (2*AtomicMassUnit, GeV) 00266 static const G4double me=.00051099892; // electron mass in GeV 00267 static const G4double me2=me*me; // Squared mass of an electron in GeV^2 00268 static const G4double thresh=me+me2/dmN; // Universal threshold in GeV 00269 // --------- 00270 //static const G4double infEn = 9.e27; 00271 G4double dN=0.; 00272 if(Z>0||N>0) dN=thresh*GeV; // @@ if upgraded, change it in a total cross section 00273 //@@ "dN=me+me2/G4NucleiProperties::GetNuclearMass(<G4double>(Z+N),<G4double>(Z)/GeV" 00274 return dN; 00275 }