#include <G4QANuENuclearCrossSection.hh>
Inheritance diagram for G4QANuENuclearCrossSection:
Public Member Functions | |
~G4QANuENuclearCrossSection () | |
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 | |
G4QANuENuclearCrossSection () |
Definition at line 49 of file G4QANuENuclearCrossSection.hh.
G4QANuENuclearCrossSection::G4QANuENuclearCrossSection | ( | ) | [inline, protected] |
G4QANuENuclearCrossSection::~G4QANuENuclearCrossSection | ( | ) |
G4double G4QANuENuclearCrossSection::CalculateCrossSection | ( | G4bool | CS, | |
G4int | F, | |||
G4int | I, | |||
G4int | PDG, | |||
G4int | Z, | |||
G4int | N, | |||
G4double | Momentum | |||
) | [virtual] |
Implements G4VQCrossSection.
Definition at line 278 of file G4QANuENuclearCrossSection.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() !! 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 antiNuEnergy 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 G4double oldL=lastEN[sep]; 00347 if(lastE<=oldL) sep-=newran; 00348 else sep+=newran; 00349 #ifdef pdebug 00350 G4cout<<"G4ANuE::CCS:n="<<newran<<",s="<<sep<<",E="<<lastE<<",oE="<<oldL<<G4endl; 00351 #endif 00352 ran=newran; 00353 chk=chk+chk; 00354 } 00355 if(chk+chk!=mL) G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Table! mL="<<mL<<G4endl; 00356 G4double lowE=lastEN[sep]; 00357 G4double highE=lastEN[sep+1]; 00358 G4double lowTX=lastTX[sep]; 00359 if(lastE<lowE||sep>=mL||lastE>highE) 00360 G4cerr<<"*Warn*G4ANuENuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE 00361 <<", sep="<<sep<<", mL="<<mL<<G4endl; 00362 lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E 00363 if(!onlyCS) // Skip the differential cross-section parameters 00364 { 00365 G4double lowQE=lastQE[sep]; 00366 lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE; 00367 #ifdef pdebug 00368 G4cout<<"G4ANuENuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl; 00369 #endif 00370 } 00371 } 00372 else 00373 { 00374 lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking... 00375 lastQEL=lastQE[mL]; 00376 } 00377 if(lastQEL<0.) lastQEL = 0.; 00378 if(lastSig<0.) lastSig = 0.; 00379 // The cross-sections are expected to be in mb 00380 lastSig*=mb38; 00381 if(!onlyCS) lastQEL*=mb38; 00382 return lastSig; 00383 }
G4double G4QANuENuclearCrossSection::GetCrossSection | ( | G4bool | fCS, | |
G4double | pMom, | |||
G4int | tgZ, | |||
G4int | tgN, | |||
G4int | pPDG = 0 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 90 of file G4QANuENuclearCrossSection.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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"---G4QAENCrossSec::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<<"G4QAENCS::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<<"G4QAENCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl; 00192 #endif 00193 if(pEn>lastTH) 00194 { 00195 #ifdef pdebug 00196 G4cout<<"G4QAENCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00197 #endif 00198 lastTH=pEn; 00199 } 00200 } 00201 #ifdef pdebug 00202 G4cout<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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<<"G4QAENCS::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 770 of file G4QANuENuclearCrossSection.cc.
00771 { 00772 G4double f=Q2/4.62; 00773 G4double ff=f*f; 00774 G4double r=ff*ff; 00775 G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3))); 00776 //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par) 00777 return 1.-s_value*(1.-s_value/2); 00778 }
G4int G4QANuENuclearCrossSection::GetExchangePDGCode | ( | ) | [virtual] |
G4double G4QANuENuclearCrossSection::GetLastQELCS | ( | ) | [inline, virtual] |
G4double G4QANuENuclearCrossSection::GetLastTOTCS | ( | ) | [inline, virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 781 of file G4QANuENuclearCrossSection.cc.
00782 { 00783 return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space 00784 }
G4double G4QANuENuclearCrossSection::GetNQE_ExchangeQ2 | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 538 of file G4QANuENuclearCrossSection.cc.
References G4UniformRand.
00539 { 00540 static const double mpi=.13957018; // charged pi meson mass in GeV 00541 static const G4double me=.00051099892;// electron mass in GeV 00542 static const G4double me2=me*me; // Squared mass of an electron in GeV^2 00543 static const G4double hme2=me2/2; // .5*m_e^2 in GeV^2 00544 static const double MN=.931494043; // Nucleon mass (inside nucleus,atomicMassUnit,GeV) 00545 static const double MN2=MN*MN; // M_N^2 in GeV^2 00546 static const double dMN=MN+MN; // 2*M_N in GeV 00547 static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic 00548 static const G4int power=7; // direct power for the magic variable 00549 static const G4double pconv=1./power; // conversion power for the magic variable 00550 static const G4int nX=21; // #Of point in the Xl table (in GeV^2) 00551 static const G4int lX=nX-1; // index of the last in the Xl table 00552 static const G4int bX=lX-1; // @@ index of the before last in the Xl table 00553 static const G4int nE=20; // #Of point in the El table (in GeV^2) 00554 static const G4int bE=nE-1; // index of the last in the El table 00555 static const G4int pE=bE-1; // index of the before last in the El table 00556 // Reversed table 00557 static const G4double X0[nX]={5.21412e-05, 00558 .437860, .681908, .891529, 1.08434, 1.26751, 1.44494, 1.61915, 1.79198, 1.96493, 2.13937, 00559 2.31664, 2.49816, 2.68559, 2.88097, 3.08705, 3.30774, 3.54917, 3.82233, 4.15131, 4.62182}; 00560 static const G4double X1[nX]={.00102591, 00561 1.00443, 1.55828, 2.03126, 2.46406, 2.87311, 3.26723, 3.65199, 4.03134, 4.40835, 4.78561, 00562 5.16549, 5.55031, 5.94252, 6.34484, 6.76049, 7.19349, 7.64917, 8.13502, 8.66246, 9.25086}; 00563 static const G4double X2[nX]={.0120304, 00564 2.59903, 3.98637, 5.15131, 6.20159, 7.18024, 8.10986, 9.00426, 9.87265, 10.7217, 11.5564, 00565 12.3808, 13.1983, 14.0116, 14.8234, 15.6359, 16.4515, 17.2723, 18.1006, 18.9386, 19.7892}; 00566 static const G4double X3[nX]={.060124, 00567 5.73857, 8.62595, 10.9849, 13.0644, 14.9636, 16.7340, 18.4066, 20.0019, 21.5342, 23.0142, 00568 24.4497, 25.8471, 27.2114, 28.5467, 29.8564, 31.1434, 32.4102, 33.6589, 34.8912, 36.1095}; 00569 static const G4double X4[nX]={.0992363, 00570 8.23746, 12.1036, 15.1740, 17.8231, 20.1992, 22.3792, 24.4092, 26.3198, 28.1320, 29.8615, 00571 31.5200, 33.1169, 34.6594, 36.1536, 37.6044, 39.0160, 40.3920, 41.7353, 43.0485, 44.3354}; 00572 static const G4double X5[nX]={.0561127, 00573 7.33661, 10.5694, 13.0778, 15.2061, 17.0893, 18.7973, 20.3717, 21.8400, 23.2211, 24.5291, 00574 25.7745, 26.9655, 28.1087, 29.2094, 30.2721, 31.3003, 32.2972, 33.2656, 34.2076, 35.1265}; 00575 static const G4double X6[nX]={.0145859, 00576 4.81774, 6.83565, 8.37399, 9.66291, 10.7920, 11.8075, 12.7366, 13.5975, 14.4025, 15.1608, 00577 15.8791, 16.5628, 17.2162, 17.8427, 18.4451, 19.0259, 19.5869, 20.1300, 20.6566, 21.1706}; 00578 static const G4double X7[nX]={.00241155, 00579 2.87095, 4.02492, 4.89243, 5.61207, 6.23747, 6.79613, 7.30433, 7.77270, 8.20858, 8.61732, 00580 9.00296, 9.36863, 9.71682, 10.0495, 10.3684, 10.6749, 10.9701, 11.2550, 11.5306, 11.7982}; 00581 static const G4double X8[nX]={.000316863, 00582 1.76189, 2.44632, 2.95477, 3.37292, 3.73378, 4.05420, 4.34415, 4.61009, 4.85651, 5.08666, 00583 5.30299, 5.50738, 5.70134, 5.88609, 6.06262, 6.23178, 6.39425, 6.55065, 6.70149, 6.84742}; 00584 static const G4double X9[nX]={3.73544e-05, 00585 1.17106, 1.61289, 1.93763, 2.20259, 2.42976, 2.63034, 2.81094, 2.97582, 3.12796, 3.26949, 00586 3.40202, 3.52680, 3.64482, 3.75687, 3.86360, 3.96557, 4.06323, 4.15697, 4.24713, 4.33413}; 00587 static const G4double XA[nX]={4.19131e-06, 00588 .849573, 1.16208, 1.38955, 1.57379, 1.73079, 1.86867, 1.99221, 2.10451, 2.20770, 2.30332, 00589 2.39252, 2.47622, 2.55511, 2.62977, 2.70066, 2.76818, 2.83265, 2.89437, 2.95355, 3.01051}; 00590 static const G4double XB[nX]={4.59981e-07, 00591 .666131, .905836, 1.07880, 1.21796, 1.33587, 1.43890, 1.53080, 1.61399, 1.69011, 1.76040, 00592 1.82573, 1.88682, 1.94421, 1.99834, 2.04959, 2.09824, 2.14457, 2.18878, 2.23107, 2.27162}; 00593 static const G4double XC[nX]={4.99861e-08, 00594 .556280, .752730, .893387, 1.00587, 1.10070, 1.18317, 1.25643, 1.32247, 1.38269, 1.43809, 00595 1.48941, 1.53724, 1.58203, 1.62416, 1.66391, 1.70155, 1.73728, 1.77128, 1.80371, 1.83473}; 00596 static const G4double XD[nX]={5.40832e-09, 00597 .488069, .657650, .778236, .874148, .954621, 1.02432, 1.08599, 1.14138, 1.19172, 1.23787, 00598 1.28049, 1.32008, 1.35705, 1.39172, 1.42434, 1.45514, 1.48429, 1.51197, 1.53829, 1.56339}; 00599 static const G4double XE[nX]={5.84029e-10, 00600 .445057, .597434, .705099, .790298, .861468, .922865, .976982, 1.02542, 1.06930, 1.10939, 00601 1.14630, 1.18050, 1.21233, 1.24208, 1.27001, 1.29630, 1.32113, 1.34462, 1.36691, 1.38812}; 00602 static const G4double XF[nX]={6.30137e-11, 00603 .418735, .560003, .659168, .737230, .802138, .857898, .906854, .950515, .989915, 1.02580, 00604 1.05873, 1.08913, 1.11734, 1.14364, 1.16824, 1.19133, 1.21306, 1.23358, 1.25298, 1.27139}; 00605 static const G4double XG[nX]={6.79627e-12, 00606 .405286, .539651, .633227, .706417, .766929, .818642, .863824, .903931, .939963, .972639, 00607 1.00250, 1.02995, 1.05532, 1.07887, 1.10082, 1.12134, 1.14058, 1.15867, 1.17572, 1.19183}; 00608 static const G4double XH[nX]={7.32882e-13, 00609 .404391, .535199, .625259, .695036, .752243, .800752, .842823, .879906, .912994, .942802, 00610 .969862, .994583, 1.01729, 1.03823, 1.05763, 1.07566, 1.09246, 1.10816, 1.12286, 1.13667}; 00611 static const G4double XI[nX]={7.90251e-14, 00612 .418084, .548382, .636489, .703728, .758106, .803630, .842633, .876608, .906576, .933269, 00613 .957233, .978886, .998556, 1.01651, 1.03295, 1.04807, 1.06201, 1.07489, 1.08683, 1.09792}; 00614 static const G4double XJ[nX]={8.52083e-15, 00615 .447299, .579635, .666780, .731788, .783268, .825512, .861013, .891356, .917626, .940597, 00616 .960842, .978802, .994820, 1.00917, 1.02208, 1.03373, 1.04427, 1.05383, 1.06253, 1.07046}; 00617 // Direct table 00618 static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0], 00619 X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]}; 00620 static const G4double dX[nE]={ 00621 (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX, 00622 (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX, 00623 (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX, 00624 (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX, 00625 (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX}; 00626 static const G4double* Xl[nE]= 00627 {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ}; 00628 static const G4double I0[nX]={0, 00629 .354631, 1.08972, 2.05138, 3.16564, 4.38343, 5.66828, 6.99127, 8.32858, 9.65998, 10.9680, 00630 12.2371, 13.4536, 14.6050, 15.6802, 16.6686, 17.5609, 18.3482, 19.0221, 19.5752, 20.0000}; 00631 static const G4double I1[nX]={0, 00632 .281625, .877354, 1.67084, 2.60566, 3.64420, 4.75838, 5.92589, 7.12829, 8.34989, 9.57708, 00633 10.7978, 12.0014, 13.1781, 14.3190, 15.4162, 16.4620, 17.4496, 18.3724, 19.2245, 20.0000}; 00634 static const G4double I2[nX]={0, 00635 .201909, .642991, 1.24946, 1.98463, 2.82370, 3.74802, 4.74263, 5.79509, 6.89474, 8.03228, 00636 9.19947, 10.3889, 11.5938, 12.8082, 14.0262, 15.2427, 16.4527, 17.6518, 18.8356, 20.0000}; 00637 static const G4double I3[nX]={0, 00638 .140937, .461189, .920216, 1.49706, 2.17728, 2.94985, 3.80580, 4.73758, 5.73867, 6.80331, 00639 7.92637, 9.10316, 10.3294, 11.6013, 12.9150, 14.2672, 15.6548, 17.0746, 18.5239, 20.0000}; 00640 static const G4double I4[nX]={0, 00641 .099161, .337358, .694560, 1.16037, 1.72761, 2.39078, 3.14540, 3.98768, 4.91433, 5.92245, 00642 7.00942, 8.17287, 9.41060, 10.7206, 12.1010, 13.5500, 15.0659, 16.6472, 18.2924, 20.0000}; 00643 static const G4double I5[nX]={0, 00644 .071131, .255084, .543312, .932025, 1.41892, 2.00243, 2.68144, 3.45512, 4.32283, 5.28411, 00645 6.33859, 7.48602, 8.72621, 10.0590, 11.4844, 13.0023, 14.6128, 16.3158, 18.1115, 20.0000}; 00646 static const G4double I6[nX]={0, 00647 .053692, .202354, .443946, .778765, 1.20774, 1.73208, 2.35319, 3.07256, 3.89177, 4.81249, 00648 5.83641, 6.96528, 8.20092, 9.54516, 10.9999, 12.5670, 14.2486, 16.0466, 17.9630, 20.0000}; 00649 static const G4double I7[nX]={0, 00650 .043065, .168099, .376879, .672273, 1.05738, 1.53543, 2.10973, 2.78364, 3.56065, 4.44429, 00651 5.43819, 6.54610, 7.77186, 9.11940, 10.5928, 12.1963, 13.9342, 15.8110, 17.8313, 20.0000}; 00652 static const G4double I8[nX]={0, 00653 .036051, .143997, .327877, .592202, .941572, 1.38068, 1.91433, 2.54746, 3.28517, 4.13277, 00654 5.09574, 6.17984, 7.39106, 8.73568, 10.2203, 11.8519, 13.6377, 15.5854, 17.7033, 20.0000}; 00655 static const G4double I9[nX]={0, 00656 .030977, .125727, .289605, .528146, .846967, 1.25183, 1.74871, 2.34384, 3.04376, 3.85535, 00657 4.78594, 5.84329, 7.03567, 8.37194, 9.86163, 11.5150, 13.3430, 15.3576, 17.5719, 20.0000}; 00658 static const G4double IA[nX]={0, 00659 .027129, .111420, .258935, .475812, .768320, 1.14297, 1.60661, 2.16648, 2.83034, 3.60650, 00660 4.50394, 5.53238, 6.70244, 8.02569, 9.51488, 11.1841, 13.0488, 15.1264, 17.4362, 20.0000}; 00661 static const G4double IB[nX]={0, 00662 .024170, .100153, .234345, .433198, .703363, 1.05184, 1.48607, 2.01409, 2.64459, 3.38708, 00663 4.25198, 5.25084, 6.39647, 7.70319, 9.18708, 10.8663, 12.7617, 14.8968, 17.2990, 20.0000}; 00664 static const G4double IC[nX]={0, 00665 .021877, .091263, .214670, .398677, .650133, .976322, 1.38510, 1.88504, 2.48555, 3.19709, 00666 4.03129, 5.00127, 6.12184, 7.40989, 8.88482, 10.5690, 12.4888, 14.6748, 17.1638, 20.0000}; 00667 static const G4double ID[nX]={0, 00668 .020062, .084127, .198702, .370384, .606100, .913288, 1.30006, 1.77535, 2.34912, 3.03253, 00669 3.83822, 4.78063, 5.87634, 7.14459, 8.60791, 10.2929, 12.2315, 14.4621, 17.0320, 20.0000}; 00670 static const G4double IE[nX]={0, 00671 .018547, .078104, .185102, .346090, .567998, .858331, 1.22535, 1.67824, 2.22735, 2.88443, 00672 3.66294, 4.57845, 5.64911, 6.89637, 8.34578, 10.0282, 11.9812, 14.2519, 16.8993, 20.0000}; 00673 static const G4double IF[nX]={0, 00674 .017143, .072466, .172271, .323007, .531545, .805393, 1.15288, 1.58338, 2.10754, 2.73758, 00675 3.48769, 4.37450, 5.41770, 6.64092, 8.07288, 9.74894, 11.7135, 14.0232, 16.7522, 20.0000}; 00676 static const G4double IG[nX]={0, 00677 .015618, .066285, .158094, .297316, .490692, .745653, 1.07053, 1.47479, 1.96931, 2.56677, 00678 3.28205, 4.13289, 5.14068, 6.33158, 7.73808, 9.40133, 11.3745, 13.7279, 16.5577, 20.0000}; 00679 static const G4double IH[nX]={0, 00680 .013702, .058434, .139923, .264115, .437466, .667179, .961433, 1.32965, 1.78283, 2.33399, 00681 2.99871, 3.79596, 4.74916, 5.88771, 7.24937, 8.88367, 10.8576, 13.2646, 16.2417, 20.0000}; 00682 static const G4double II[nX]={0, 00683 .011264, .048311, .116235, .220381, .366634, .561656, .813132, 1.13008, 1.52322, 2.00554, 00684 2.59296, 3.30542, 4.16834, 5.21490, 6.48964, 8.05434, 9.99835, 12.4580, 15.6567, 20.0000}; 00685 static const G4double IJ[nX]={0, 00686 .008628, .037206, .089928, .171242, .286114, .440251, .640343, .894382, 1.21208, 1.60544, 00687 2.08962, 2.68414, 3.41486, 4.31700, 5.44048, 6.85936, 8.69067, 11.1358, 14.5885, 20.0000}; 00688 static const G4double* Il[nE]= 00689 {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ}; 00690 static const G4double lE[nE]={ 00691 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292, 00692 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218}; 00693 static const G4double lEmi=lE[0]; 00694 static const G4double lEma=lE[nE-1]; 00695 static const G4double dlE=(lEma-lEmi)/bE; 00696 //*************************************************************************************** 00697 G4double Enu=lastE; // Get energy of the last calculated cross-section 00698 G4double lEn=std::log(Enu); // log(E) for interpolation 00699 G4double rE=(lEn-lEmi)/dlE; // Position of the energy 00700 G4int fE=static_cast<int>(rE); // Left bin for interpolation 00701 if(fE<0) fE=0; 00702 if(fE>pE)fE=pE; 00703 G4int sE=fE+1; // Right bin for interpolation 00704 G4double dE=rE-fE; // relative log shift from the left bin 00705 G4double dEnu=Enu+Enu; // doubled energy of nu/anu 00706 G4double Enu2=Enu*Enu; // squared energy of nu/anu 00707 G4double Ee=Enu-me; // Free Energy of neutrino/anti-neutrino 00708 G4double ME=Enu*MN; // M*E 00709 G4double dME=ME+ME; // 2*M*E 00710 G4double dEMN=(dEnu+MN)*ME; 00711 G4double MEm=ME-hme2; 00712 G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2); 00713 G4double E2M=MN*Enu2-(Enu+MN)*hme2; 00714 G4double ymax=(E2M+sqE)/dEMN; 00715 G4double ymin=(E2M-sqE)/dEMN; 00716 G4double rmin=1.-ymin; 00717 G4double rhm2E=hme2/Enu2; 00718 G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu) 00719 G4double Q2ma=dME*ymax; // Q2_max(E_nu) 00720 G4double Q2nq=Ee*dMN-mcV; 00721 if(Q2ma>Q2nq) Q2ma=Q2nq; // Correction for Non Quasi Elastic 00722 // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma --- 00723 G4double Rmi=Q2mi/Q2ma; 00724 G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78); 00725 // --- E-interpolation must be done in a log scale --- 00726 G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu) 00727 G4double Xma=std::pow((shift-1.),power); // X_max(E_nu) 00728 // Find the integral values integ(Xmi) & integ(Xma) using the direct table 00729 G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step 00730 G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum 00731 G4double rXi=(Xmi-iXmi)/idX; 00732 G4int iXi=static_cast<int>(rXi); 00733 if(iXi<0) iXi=0; 00734 if(iXi>bX) iXi=bX; 00735 G4double dXi=rXi-iXi; 00736 G4double bntil=Il[fE][iXi]; 00737 G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil); 00738 G4double bntir=Il[sE][iXi]; 00739 G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir); 00740 G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral 00741 // 00742 G4double rXa=(Xma-iXmi)/idX; 00743 G4int iXa=static_cast<int>(rXa); 00744 if(iXa<0) iXa=0; 00745 if(iXa>bX) iXa=bX; 00746 G4double dXa=rXa-iXa; 00747 G4double bntal=Il[fE][iXa]; 00748 G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal); 00749 G4double bntar=Il[sE][iXa]; 00750 G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar); 00751 G4double inta=intal+dE*(intar-intal);// interpolated end of the integral 00752 // 00753 // *** Find X using the reversed table *** 00754 G4double intx=inti+(inta-inti)*G4UniformRand(); 00755 G4int intc=static_cast<int>(intx); 00756 if(intc<0) intc=0; 00757 if(intc>bX) intc=bX; 00758 G4double dint=intx-intc; 00759 G4double mXl=Xl[fE][intc]; 00760 G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl); 00761 G4double mXr=Xl[sE][intc]; 00762 G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr); 00763 G4double X=Xlb+dE*(Xrb-Xlb); // interpolated X value 00764 G4double R=shift-std::pow(X,pconv); 00765 G4double Q2=R*Q2ma; 00766 return Q2*GeV*GeV; 00767 }
G4VQCrossSection * G4QANuENuclearCrossSection::GetPointer | ( | ) | [static] |
Definition at line 72 of file G4QANuENuclearCrossSection.cc.
Referenced by G4QInelastic::GetMeanFreePath(), and G4QInelastic::PostStepDoIt().
00073 { 00074 static G4QANuENuclearCrossSection theCrossSection;//**Static body of the Cross Section** 00075 return &theCrossSection; 00076 }
G4double G4QANuENuclearCrossSection::GetQEL_ExchangeQ2 | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 453 of file G4QANuENuclearCrossSection.cc.
References G4UniformRand.
00454 { 00455 static const G4double me=.00051099892; // electron mass in GeV 00456 static const G4double me2=me*me; // Squared mass of an electron in GeV^2 00457 static const G4double hme2=me2/2; // .5*m_e^2 in GeV^2 00458 static const double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV) 00459 static const double MN2=MN*MN; // M_N^2 in GeV^2 00460 static const G4double power=-3.5; // direct power for the magic variable 00461 static const G4double pconv=1./power;// conversion power for the magic variable 00462 static const G4int nQ2=101; // #Of point in the Q2l table (in GeV^2) 00463 static const G4int lQ2=nQ2-1; // index of the last in the Q2l table 00464 static const G4int bQ2=lQ2-1; // index of the before last in the Q2 ltable 00465 // Reversed table 00466 static const G4double Xl[nQ2]={5.20224e-16, 00467 .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898, 00468 .093878, .103332, .112844, .122410, .132023, .141680, .151376, .161109, .170875, .180672, 00469 .190499, .200352, .210230, .220131, .230055, .239999, .249963, .259945, .269944, .279960, 00470 .289992, .300039, .310099, .320173, .330260, .340359, .350470, .360592, .370724, .380867, 00471 .391019, .401181, .411352, .421531, .431719, .441915, .452118, .462329, .472547, .482771, 00472 .493003, .503240, .513484, .523734, .533989, .544250, .554517, .564788, .575065, .585346, 00473 .595632, .605923, .616218, .626517, .636820, .647127, .657438, .667753, .678072, .688394, 00474 .698719, .709048, .719380, .729715, .740053, .750394, .760738, .771085, .781434, .791786, 00475 .802140, .812497, .822857, .833219, .843582, .853949, .864317, .874687, .885060, .895434, 00476 .905810, .916188, .926568, .936950, .947333, .957719, .968105, .978493, .988883, .999275}; 00477 // Direct table 00478 static const G4double Xmax=Xl[lQ2]; 00479 static const G4double Xmin=Xl[0]; 00480 static const G4double dX=(Xmax-Xmin)/lQ2; // step in X(Q2, GeV^2) 00481 static const G4double inl[nQ2]={0, 00482 1.52225, 2.77846, 3.96651, 5.11612, 6.23990, 7.34467, 8.43466, 9.51272, 10.5809, 11.6406, 00483 12.6932, 13.7394, 14.7801, 15.8158, 16.8471, 17.8743, 18.8979, 19.9181, 20.9353, 21.9496, 00484 22.9614, 23.9707, 24.9777, 25.9826, 26.9855, 27.9866, 28.9860, 29.9837, 30.9798, 31.9745, 00485 32.9678, 33.9598, 34.9505, 35.9400, 36.9284, 37.9158, 38.9021, 39.8874, 40.8718, 41.8553, 00486 42.8379, 43.8197, 44.8007, 45.7810, 46.7605, 47.7393, 48.7174, 49.6950, 50.6718, 51.6481, 00487 52.6238, 53.5990, 54.5736, 55.5476, 56.5212, 57.4943, 58.4670, 59.4391, 60.4109, 61.3822, 00488 62.3531, 63.3236, 64.2937, 65.2635, 66.2329, 67.2019, 68.1707, 69.1390, 70.1071, 71.0748, 00489 72.0423, 73.0095, 73.9763, 74.9429, 75.9093, 76.8754, 77.8412, 78.8068, 79.7721, 80.7373, 00490 81.7022, 82.6668, 83.6313, 84.5956, 85.5596, 86.5235, 87.4872, 88.4507, 89.4140, 90.3771, 00491 91.3401, 92.3029, 93.2656, 94.2281, 95.1904, 96.1526, 97.1147, 98.0766, 99.0384, 100.000}; 00492 G4double Enu=lastE; // Get energy of the last calculated cross-section 00493 G4double dEnu=Enu+Enu; // doubled energy of nu/anu 00494 G4double Enu2=Enu*Enu; // squared energy of nu/anu 00495 G4double ME=Enu*MN; // M*E 00496 G4double dME=ME+ME; // 2*M*E 00497 G4double dEMN=(dEnu+MN)*ME; 00498 G4double MEm=ME-hme2; 00499 G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2); 00500 G4double E2M=MN*Enu2-(Enu+MN)*hme2; 00501 G4double ymax=(E2M+sqE)/dEMN; 00502 G4double ymin=(E2M-sqE)/dEMN; 00503 G4double rmin=1.-ymin; 00504 G4double rhm2E=hme2/Enu2; 00505 G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu) 00506 G4double Q2ma=dME*ymax; // Q2_max(E_nu) 00507 G4double Xma=std::pow((1.+Q2mi),power); // X_max(E_nu) 00508 G4double Xmi=std::pow((1.+Q2ma),power); // X_min(E_nu) 00509 // Find the integral values integ(Xmi) & integ(Xma) using the direct table 00510 G4double rXi=(Xmi-Xmin)/dX; 00511 G4int iXi=static_cast<int>(rXi); 00512 if(iXi<0) iXi=0; 00513 if(iXi>bQ2) iXi=bQ2; 00514 G4double dXi=rXi-iXi; 00515 G4double bnti=inl[iXi]; 00516 G4double inti=bnti+dXi*(inl[iXi+1]-bnti); 00517 // 00518 G4double rXa=(Xma-Xmin)/dX; 00519 G4int iXa=static_cast<int>(rXa); 00520 if(iXa<0) iXa=0; 00521 if(iXa>bQ2) iXa=bQ2; 00522 G4double dXa=rXa-iXa; 00523 G4double bnta=inl[iXa]; 00524 G4double inta=bnta+dXa*(inl[iXa+1]-bnta); 00525 // *** Find X using the reversed table *** 00526 G4double intx=inti+(inta-inti)*G4UniformRand(); 00527 G4int intc=static_cast<int>(intx); 00528 if(intc<0) intc=0; 00529 if(intc>bQ2) intc=bQ2; // If it is more than max, then the BAD extrapolation 00530 G4double dint=intx-intc; 00531 G4double mX=Xl[intc]; 00532 G4double X=mX+dint*(Xl[intc+1]-mX); 00533 G4double Q2=std::pow(X,pconv)-1.; 00534 return Q2*GeV*GeV; 00535 }
G4double G4QANuENuclearCrossSection::ThresholdEnergy | ( | G4int | Z, | |
G4int | N, | |||
G4int | PDG = -12 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 259 of file G4QANuENuclearCrossSection.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 }