#include <G4QElectronNuclearCrossSection.hh>
Inheritance diagram for G4QElectronNuclearCrossSection:
Public Member Functions | |
~G4QElectronNuclearCrossSection () | |
G4double | ThresholdEnergy (G4int Z, G4int N, G4int PDG=11) |
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 | GetExchangeEnergy () |
G4double | GetVirtualFactor (G4double nu, G4double Q2) |
G4double | GetExchangeQ2 (G4double nu) |
Static Public Member Functions | |
static G4VQCrossSection * | GetPointer () |
Protected Member Functions | |
G4QElectronNuclearCrossSection () |
Definition at line 46 of file G4QElectronNuclearCrossSection.hh.
G4QElectronNuclearCrossSection::G4QElectronNuclearCrossSection | ( | ) | [inline, protected] |
G4QElectronNuclearCrossSection::~G4QElectronNuclearCrossSection | ( | ) |
Definition at line 78 of file G4QElectronNuclearCrossSection.cc.
00079 { 00080 G4int lens=J1->size(); 00081 for(G4int i=0; i<lens; ++i) delete[] (*J1)[i]; 00082 delete J1; 00083 G4int hens=J2->size(); 00084 for(G4int i=0; i<hens; ++i) delete[] (*J2)[i]; 00085 delete J2; 00086 G4int pens=J3->size(); 00087 for(G4int i=0; i<pens; ++i) delete[] (*J3)[i]; 00088 delete J3; 00089 }
G4double G4QElectronNuclearCrossSection::CalculateCrossSection | ( | G4bool | CS, | |
G4int | F, | |||
G4int | I, | |||
G4int | PDG, | |||
G4int | Z, | |||
G4int | N, | |||
G4double | Momentum | |||
) | [virtual] |
Implements G4VQCrossSection.
Definition at line 318 of file G4QElectronNuclearCrossSection.cc.
References G4cerr, G4cout, and G4endl.
Referenced by GetCrossSection().
00320 { 00321 static const G4int nE=336; // !! If change this, change it in GetFunctions() (*.hh) !! 00322 static const G4int mL=nE-1; 00323 static const G4double EMi=2.0612; // Minimum tabulated Energy of the Electron 00324 static const G4double EMa=50000.; // Maximum tabulated Energy of the Electron 00325 static const G4double lEMi=std::log(EMi); // Min tabulated logarithmic Energy of Electron 00326 static const G4double lEMa=std::log(EMa); // Max tabulated logarithmic Energy of Electron 00327 static const G4double dlnE=(lEMa-lEMi)/mL; // Log step in the table for electron Energy 00328 static const G4double alop=1./137.036/3.14159265; //coefficient for Ee>50000 calculations 00329 static const G4double mel=0.5109989; // Mass of the electron in MeV 00330 static const G4double mel2=mel*mel;// Squared Mass of the electron in MeV 00331 static const G4double lmel=std::log(mel); // Log of the electron mass 00332 // *** Begin of the Associative memory for acceleration of the cross section calculations 00333 static std::vector <G4int> colF; // Vector of LastStartPosition in Ji-function tables 00334 static std::vector <G4double> colH;// Vector of HighEnergyCoefficients (functional) 00335 #ifdef pdebug 00336 G4cout<<"G4QElectronNucCrossSection::CalculateCrossSection: ***Called*** "<<J3->size(); 00337 if(J3->size()) G4cout<<", p="<<(*J3)[0]; 00338 G4cout<<G4endl; 00339 #endif 00340 // *** End of Static Definitions (Associative Memory) *** 00341 //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Electron 00342 onlyCS=CS; // Flag to calculate only CS (not Si/Bi) 00343 G4double TotEnergy2=Momentum*Momentum+mel2; 00344 G4double TotEnergy=std::sqrt(TotEnergy2); // Total energy of the electron 00345 lastE=TotEnergy-mel; // Kinetic energy of the electron 00346 #ifdef pdebug 00347 G4cout<<"G4QElectronNucCS::CalcCS: P="<<Momentum<<", F="<<F<<", I="<<I<<", Z="<<targZ; 00348 if(J3->size()) G4cout<<", p="<<(*J3)[0]; 00349 G4cout<<", N="<<targN<<", onlyCS="<<CS<<",E="<<lastE<<",th="<<EMi<<G4endl; 00350 #endif 00351 G4double A=targN+targZ; // New A (can differ from G4double targetAtomicNumber) 00352 if(F<=0) // This isotope was not the last used isotop 00353 { 00354 if(F<0) // This isotope was found in DAMDB =-------=> RETRIEVE 00355 { // ...........................................=------= 00356 if (lastE<=EMi) // Energy is below the minimum energy in the table 00357 { 00358 lastE=0.; 00359 lastG=0.; 00360 lastSig=0.; 00361 #ifdef pdebug 00362 G4cout<<"G4QElectronNucCS::CalcCS: Old CS=0 as lastE="<<lastE<<" < "<<EMi<<G4endl; 00363 #endif 00364 return 0.; 00365 } 00366 lastJ1 =(*J1)[I]; // Pointer to the prepared J1 function 00367 lastJ2 =(*J2)[I]; // Pointer to the prepared J2 function 00368 lastJ3 =(*J3)[I]; // Pointer to the prepared J3 function 00369 lastF =colF[I]; // Last ZeroPosition in the J-functions 00370 lastH =colH[I]; // Last High Energy Coefficient (A-dependent) 00371 } 00372 else // This isotope wasn't calculated previously => CREATE 00373 { 00374 lastJ1 = new G4double[nE]; // Allocate memory for the new J1 function 00375 lastJ2 = new G4double[nE]; // Allocate memory for the new J2 function 00376 lastJ3 = new G4double[nE]; // Allocate memory for the new J3 function 00377 lastF = GetFunctions(A,lastJ1,lastJ2,lastJ3);//newZeroPos and J-functions filling 00378 lastH = alop*A*(1.-.072*std::log(A)); // like lastSP of G4PhotonuclearCrossSection 00379 #ifdef pdebug 00380 G4cout<<"==>G4QElNCS::CalcCS: pJ1="<<lastJ1<<",pJ2="<<lastJ2<<",pJ3="<<lastJ3; 00381 if(lastJ1) G4cout<<", J1="<<lastJ1[0]<<",J2="<<lastJ2[0]<<",J3="<<lastJ3[0]; 00382 G4cout<<G4endl; 00383 #endif 00384 // *** The synchronization check *** 00385 G4int sync=J1->size(); 00386 if(sync!=I) G4cerr<<"***G4QElectNuclearCS::CalcCS: PDG=11 ,S="<<sync<<"#"<<I<<G4endl; 00387 J1->push_back(lastJ1); 00388 J2->push_back(lastJ2); 00389 J3->push_back(lastJ3); 00390 colF.push_back(lastF); 00391 colH.push_back(lastH); 00392 } // End of creation of the new set of parameters 00393 } // End of parameters udate 00394 // =----------------= NOW Calculate the Cross Section =--------------------= 00395 if (lastE<=lastTH || lastE<=EMi) // Check that muKiE is higher than ThreshE 00396 { 00397 lastE=0.; 00398 lastG=0.; 00399 lastSig=0.; 00400 #ifdef pdebug 00401 G4cout<<"G4QElectronNucCS::CalcCS:CS=0 as T="<<lastE<<"<"<<EMi<<" || "<<lastTH<<G4endl; 00402 #endif 00403 return 0.; 00404 } 00405 G4double lE=std::log(lastE); // log(muE) (it is necessary for the fit) 00406 lastG=lE-lmel; // Gamma of the electron (used to recover log(eE)) 00407 G4double dlg1=lastG+lastG-1.; 00408 G4double lgoe=lastG/lastE; 00409 if(lE<lEMa) // Log fit is made explicitly to fix the last bin for the randomization 00410 { 00411 G4double shift=(lE-lEMi)/dlnE; 00412 G4int blast=static_cast<int>(shift); 00413 #ifdef pdebug 00414 G4cout<<"-->G4QElectronNuclearCS::CalcCrossSect:LOGfit b="<<blast<<",max="<<mL<<",lJ1=" 00415 <<lastJ1<<",lJ2="<<lastJ2<<",lJ3="<<lastJ3<<",lEmin="<<lEMi<<",d="<<dlnE<<G4endl; 00416 #endif 00417 if(blast<0) blast=0; 00418 if(blast>=mL) blast=mL-1; 00419 shift-=blast; 00420 lastL=blast+1; 00421 G4double YNi=dlg1*lastJ1[blast] 00422 -lgoe*(lastJ2[blast]+lastJ2[blast]-lastJ3[blast]/lastE); 00423 G4double YNj=dlg1*lastJ1[lastL] 00424 -lgoe*(lastJ2[lastL]+lastJ2[lastL]-lastJ3[lastL]/lastE); 00425 lastSig= YNi+shift*(YNj-YNi); 00426 if(lastSig>YNj)lastSig=YNj; 00427 #ifdef pdebug 00428 G4cout<<"G4QElectNucCS::CalcCS:S="<<lastSig<<",lE="<<lE<<",Yi="<<YNi<<",Yj="<<YNj 00429 <<",J1="<<lastJ1[blast]<<",J2="<<lastJ2[blast]<<",J3="<<lastJ3[blast]<<G4endl; 00430 G4cout<<"G4QElectNucCS::CalcCS:s="<<shift<<",Jb="<<lastJ1[blast]<<",J="<<lastJ1[lastL]; 00431 if(J3->size()) G4cout<<", p="<<(*J3)[0]; 00432 G4cout<<",b="<<blast<<",lEmax="<<lEMa<<",lgoe="<<lgoe<<G4endl; 00433 #endif 00434 } 00435 else 00436 { 00437 #ifdef pdebug 00438 G4cout<<"->G4QElecNucCS::CCS:LOGex="<<lastJ1<<",lJ2="<<lastJ2<<",lJ3="<<lastJ3<<G4endl; 00439 #endif 00440 lastL=mL; 00441 G4double term1=lastJ1[mL]+lastH*HighEnergyJ1(lE); 00442 G4double term2=lastJ2[mL]+lastH*HighEnergyJ2(lE); 00443 G4double term3=lastJ3[mL]+lastH*HighEnergyJ3(lE); 00444 lastSig=dlg1*term1-lgoe*(term2+term2-term3/lastE); 00445 #ifdef pdebug 00446 G4cout<<"G4QElNucCS::CalculateCrossSection:S="<<lastSig<<",lE="<<lE<<",J1=" 00447 <<lastH*HighEnergyJ1(lE)<<",Pm="<<lastJ1[mL]<<",Fm="<<lastJ2[mL]<<",Fh="; 00448 if(J3->size()) G4cout<<", p="<<(*J3)[0]; 00449 G4cout<<lastH*HighEnergyJ2(lE)<<",EM="<<lEMa<<G4endl; 00450 #endif 00451 } 00452 if(lastSig<0.) lastSig = 0.; 00453 #ifdef pdebug 00454 if(J3->size()) G4cout<<"-->>->> G4QElNucCS::CalculateCrossSection: p="<<(*J3)[0]<<G4endl; 00455 #endif 00456 return lastSig; 00457 }
G4double G4QElectronNuclearCrossSection::GetCrossSection | ( | G4bool | fCS, | |
G4double | pMom, | |||
G4int | tgZ, | |||
G4int | tgN, | |||
G4int | pPDG = 0 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 93 of file G4QElectronNuclearCrossSection.cc.
References CalculateCrossSection(), G4cout, G4endl, and ThresholdEnergy().
00095 { 00096 static const G4double mel=0.5109989; // Mass of the electron in MeV 00097 static const G4double mel2=mel*mel; // Squared Mass of the electron in MeV 00098 static G4int j; // A#0f records found in DB for this projectile 00099 static std::vector <G4int> colPDG;// Vector of the projectile PDG code 00100 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) 00101 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) 00102 static std::vector <G4double> colP; // Vector of last momenta for the reaction 00103 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction 00104 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction 00105 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** 00106 G4double pEn=std::sqrt(pMom*pMom+mel2)-mel; // ==> electron/positron kinEnergy 00107 #ifdef pdebug 00108 G4cout<<"G4QENCS::GetCS:->> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN 00109 <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz=" 00110 <<colN.size()<<G4endl; 00111 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00112 #endif 00113 if(std::abs(pPDG)!=11) 00114 { 00115 #ifdef pdebug 00116 G4cout<<"G4QENCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl; 00117 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00118 #endif 00119 return 0.; // projectile PDG=0 is a mistake (?!) @@ 00120 } 00121 G4bool in=false; // By default the isotope must be found in the AMDB 00122 if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope 00123 { 00124 in = false; // By default the isotope haven't be found in AMDB 00125 lastP = 0.; // New momentum history (nothing to compare with) 00126 lastPDG = pPDG; // The last PDG of the projectile 00127 lastN = tgN; // The last N of the calculated nucleus 00128 lastZ = tgZ; // The last Z of the calculated nucleus 00129 lastI = colN.size(); // Size of the Associative Memory DB in the heap 00130 j = 0; // A#0f records found in DB for this projectile 00131 if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found 00132 { // The nucleus with projPDG is found in AMDB 00133 if(colN[i]==tgN && colZ[i]==tgZ) 00134 { 00135 lastI=i; 00136 lastTH =colTH[i]; // Last THreshold (A-dependent) 00137 #ifdef pdebug 00138 G4cout<<"G4QENCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl; 00139 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00140 #endif 00141 if(pEn<=lastTH) 00142 { 00143 #ifdef pdebug 00144 G4cout<<"G4QENCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=0"<<G4endl; 00145 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00146 #endif 00147 return 0.; // Energy is below the Threshold value 00148 } 00149 lastP =colP [i]; // Last Momentum (A-dependent) 00150 lastCS =colCS[i]; // Last CrossSect (A-dependent) 00151 // if(std::fabs(lastP/pMom-1.)<tolerance) // VI (do not use tolerance) 00152 if(lastP == pMom) 00153 { 00154 #ifdef pdebug 00155 G4cout<<"G4QENCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; 00156 #endif 00157 CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only 00158 return lastCS*millibarn; // Use theLastCS 00159 } 00160 in = true; // This is the case when the isotop is found in DB 00161 // Momentum pMom is in IU ! @@ Units 00162 #ifdef pdebug 00163 G4cout<<"G4QENCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl; 00164 #endif 00165 lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update 00166 #ifdef pdebug 00167 G4cout<<"G4QENCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl; 00168 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00169 #endif 00170 if(lastCS<=0. && pEn>lastTH) // Correct the threshold 00171 { 00172 #ifdef pdebug 00173 G4cout<<"G4QENCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00174 #endif 00175 lastTH=pEn; 00176 } 00177 break; // Go out of the LOOP 00178 } 00179 #ifdef pdebug 00180 G4cout<<"---G4QENCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i] 00181 <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl; 00182 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00183 #endif 00184 j++; // Increment a#0f records found in DB for this pPDG 00185 } 00186 if(!in) // This nucleus has not been calculated previously 00187 { 00188 #ifdef pdebug 00189 G4cout<<"G4QENCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl; 00190 #endif 00192 lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create 00193 if(lastCS<=0.) 00194 { 00195 lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last 00196 #ifdef pdebug 00197 G4cout<<"G4QENCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl; 00198 #endif 00199 if(pEn>lastTH) 00200 { 00201 #ifdef pdebug 00202 G4cout<<"G4QENCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00203 #endif 00204 lastTH=pEn; 00205 } 00206 } 00207 #ifdef pdebug 00208 G4cout<<"G4QENCS::GetCrosSec: New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl; 00209 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00210 #endif 00211 colN.push_back(tgN); 00212 colZ.push_back(tgZ); 00213 colPDG.push_back(pPDG); 00214 colP.push_back(pMom); 00215 colTH.push_back(lastTH); 00216 colCS.push_back(lastCS); 00217 #ifdef pdebug 00218 G4cout<<"G4QENCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00219 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00220 #endif 00221 return lastCS*millibarn; 00222 } // End of creation of the new set of parameters 00223 else 00224 { 00225 #ifdef pdebug 00226 G4cout<<"G4QENCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl; 00227 #endif 00228 colP[lastI]=pMom; 00229 colPDG[lastI]=pPDG; 00230 colCS[lastI]=lastCS; 00231 } 00232 } // End of parameters udate 00233 else if(pEn<=lastTH) 00234 { 00235 #ifdef pdebug 00236 G4cout<<"G4QENCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl; 00237 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00238 #endif 00239 return 0.; // Momentum is below the Threshold Value -> CS=0 00240 } 00241 // else if(std::fabs(lastP/pMom-1.)<tolerance) // VI (do not use tolerance) 00242 else if(lastP == pMom) 00243 { 00244 #ifdef pdebug 00245 G4cout<<"G4QENCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", CS="<<lastCS*millibarn<<G4endl; 00246 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00247 #endif 00248 return lastCS*millibarn; // Use theLastCS 00249 } 00250 else 00251 { 00252 #ifdef pdebug 00253 G4cout<<"G4QENCS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl; 00254 #endif 00255 lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB 00256 lastP=pMom; 00257 } 00258 #ifdef pdebug 00259 G4cout<<"G4QENCS::GetCroSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00260 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00261 #endif 00262 return lastCS*millibarn; 00263 }
G4double G4QElectronNuclearCrossSection::GetExchangeEnergy | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 2636 of file G4QElectronNuclearCrossSection.cc.
References G4cerr, G4cout, and G4UniformRand.
02637 { 02638 // @@ All constants are copy of that from PhotonCrossSection funct. => Make them general. 02639 static const G4int nE=336; // !! If change this, change it in GetFunctions() (*.hh) !! 02640 static const G4int mL=nE-1; 02641 static const G4double EMi=2.0612; // Minimum Energy 02642 static const G4double EMa=50000.; // Maximum Energy 02643 static const G4double lEMi=std::log(EMi); // Minimum logarithmic Energy 02644 static const G4double lEMa=std::log(EMa); // Maximum logarithmic Energy 02645 static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic step in Energy 02646 static const G4double mel=0.5109989; // Mass of electron in MeV 02647 static const G4double lmel=std::log(mel); // Log of electron mass 02648 G4double phLE=0.; // Prototype of the log(nu=E_gamma) 02649 G4double Y[nE]; // Prepare the array for randomization 02650 #ifdef debug 02651 G4cout<<"G4QElectronNuclearCroSect::GetExEn:"<<lastF<<",L="<<lastL<<",1="<<lastJ1[lastL] 02652 <<",2="<<lastJ2[lastL]<<",3="<<lastJ3[lastL]<<",S="<<lastSig<<",E="<<lastE<<G4endl; 02653 #endif 02654 G4double lastLE=lastG+lmel; // recover log(eE) from the gamma (lastG) 02655 G4double dlg1=lastG+lastG-1.; 02656 G4double lgoe=lastG/lastE; 02657 for(G4int i=lastF;i<=lastL;i++) 02658 Y[i]=dlg1*lastJ1[i]-lgoe*(lastJ2[i]+lastJ2[i]-lastJ3[i]/lastE); 02659 G4double ris=lastSig*G4UniformRand(); // If Sig > Y[lastL=mL] -> it is in functRegion 02660 #ifdef debug 02661 G4cout<<"G4QElectronNuclearCrossSection::GetExchEnergy: "<<ris<<",Y="<<Y[lastL]<<G4endl; 02662 #endif 02663 if(ris<Y[lastL]) // Search in the table 02664 { 02665 G4int j=lastF; 02666 G4double Yj=Y[j]; // It mast be 0 (some times just very small) 02667 while (ris>Yj && j<lastL) // Associative search 02668 { 02669 j++; 02670 Yj=Y[j]; // High value 02671 } 02672 G4int j1=j-1; 02673 G4double Yi=Y[j1]; // Low value 02674 phLE=lEMi+(j1+(ris-Yi)/(Yj-Yi))*dlnE; 02675 #ifdef debug 02676 G4cout<<"G4QElectronNucCS::GetExchEn="<<phLE<<",l="<<lEMi<<",j="<<j<<",ris="<<ris<<",Yi=" 02677 <<Yi<<",Y="<<Yj<<G4endl; 02678 #endif 02679 } 02680 else // Search with the function 02681 { 02682 if(lastL<mL)G4cerr<<"G4QElNCS::GEE:L="<<lastL<<",S="<<lastSig<<",Y="<<Y[lastL]<<G4endl; 02683 G4double f=(ris-Y[lastL])/lastH; // ScaledResidualValue of the CrossSection integral 02684 #ifdef pdebug 02685 G4cout<<"G4QElNucCS::GetExEn:HighEnergy f="<<f<<",ris="<<ris<<",lastH="<<lastH<<G4endl; 02686 #endif 02687 phLE=SolveTheEquation(f); // Solve equation to find Log(phE) (compare with lastLE) 02688 #ifdef pdebug 02689 G4cout<<"G4QElectronNuclearCS::GetExchangeEnergy: HighEnergy lphE="<<phLE<<G4endl; 02690 #endif 02691 } 02692 if(phLE>lastLE) 02693 { 02694 G4cerr<<"***G4QElNucCS::GetExEn:N="<<lastN<<",Z="<<lastZ<<", lpE"<<phLE<<">leE"<<lastLE 02695 <<",S="<<lastSig<<",rS="<<ris<<",B="<<lastF<<",E="<<lastL<<",Y="<<Y[lastL]<<G4endl; 02696 if(lastLE<7.2) phLE=std::log(std::exp(lastLE)-.511); 02697 else phLE=7.; 02698 } 02699 #ifdef pdebug 02700 G4cout<<"G4QElectronNuclearCS::GetExchangeEnergy: *** DONE ***, lphE="<<phLE<<G4endl; 02701 #endif 02702 return std::exp(phLE); 02703 } // End of GetExchangeEnergy
G4int G4QElectronNuclearCrossSection::GetExchangePDGCode | ( | ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 2806 of file G4QElectronNuclearCrossSection.cc.
Reimplemented from G4VQCrossSection.
Definition at line 2752 of file G4QElectronNuclearCrossSection.cc.
References G4cerr, and G4UniformRand.
02753 { 02754 static const G4double mel=0.5109989; // Mass of electron in MeV 02755 static const G4double mel2=mel*mel; // Squared Mass of electron in MeV 02756 G4double y=nu/lastE; // Part of energy carried by the equivalent pfoton 02757 if(y>=1.-1./(lastG+lastG)) return 0.; // The region where the method does not work 02758 G4double y2=y*y; // Squared photonic part of energy 02759 G4double ye=1.-y; // Part of energy carried by secondary electron 02760 G4double Qi2=mel2*y2/ye; // Minimum Q2 02761 G4double Qa2=4*lastE*lastE*ye; // Maximum Q2 02762 G4double iar=Qi2/Qa2; // Q2min/Q2max ratio 02763 G4double Dy=ye+.5*y2; // D(y) function 02764 G4double Py=ye/Dy; // P(y) function 02765 G4double ePy=1.-std::exp(Py); // 1-exp(P(y)) part 02766 G4double Uy=Py*(1.-iar); // U(y) function 02767 G4double Fy=(ye+ye)*(1.+ye)*iar/y2; // F(y) function 02768 G4double fr=iar/(1.-ePy*iar); // Q-fraction 02769 if(Fy<=-fr) 02770 { 02771 #ifdef edebug 02772 G4cerr<<"**G4QElNucCrossSec::GetExQ2:Fy="<<Fy<<"+fr="<<fr<<" <0"<<",iar="<<iar<<G4endl; 02773 #endif 02774 return 0.; 02775 } 02776 G4double LyQa2=std::log(Fy+fr); // L(y,Q2max) function 02777 G4bool cond=true; 02778 G4int maxTry=3; 02779 G4int cntTry=0; 02780 G4double Q2=Qi2; 02781 while(cond&&cntTry<maxTry) // The loop to avoid x>1. 02782 { 02783 G4double R=G4UniformRand(); // Random number (0,1) 02784 Q2=Qi2*(ePy+1./(std::exp(R*LyQa2-(1.-R)*Uy)-Fy)); 02785 cntTry++; 02786 cond = Q2>1878.*nu; 02787 } 02788 if(Q2<Qi2) 02789 { 02790 #ifdef edebug 02791 G4cerr<<"*G4QElectronNuclearCrossSec::GetExchangeQ2:Q2="<<Q2<<" < Q2min="<<Qi2<<G4endl; 02792 #endif 02793 return Qi2; 02794 } 02795 if(Q2>Qa2) 02796 { 02797 #ifdef edebug 02798 G4cerr<<"*G4QElectronNucCrossSection::GetExchangeQ2:Q2="<<Q2<<" > Q2max="<<Qi2<<G4endl; 02799 #endif 02800 return Qa2; 02801 } 02802 return Q2; 02803 }
G4VQCrossSection * G4QElectronNuclearCrossSection::GetPointer | ( | ) | [static] |
Definition at line 72 of file G4QElectronNuclearCrossSection.cc.
Referenced by G4QInelastic::GetMeanFreePath(), G4QAtomicElectronScattering::GetMeanFreePath(), G4QInelastic::PostStepDoIt(), and G4QAtomicElectronScattering::PostStepDoIt().
00073 { 00074 static G4QElectronNuclearCrossSection theCrossSection; //**Static body of Cross Section** 00075 return &theCrossSection; 00076 }
Reimplemented from G4VQCrossSection.
Definition at line 2808 of file G4QElectronNuclearCrossSection.cc.
References G4cerr.
02809 { 02810 static const G4double dM=938.27+939.57;// MeanDoubleNucleonMass = m_n+m_p(@@ no binding) 02811 static const G4double Q0=843.; // Coefficient of the dipole nucleonic form-factor 02812 static const G4double Q02=Q0*Q0; // SquaredCoefficient of dipoleNucleonicFormFactor 02813 static const G4double blK0=std::log(185.); // Coefficient of the b-function 02814 static const G4double bp=0.85; // Power of the b-function 02815 static const G4double clK0=std::log(1390.);// Coefficient of the c-function 02816 static const G4double cp=3.; // Power of the c-function 02817 //G4double x=Q2/dM/nu; // Direct x definition 02818 G4double K=nu-Q2/dM; // K=nu*(1-x) 02819 if(K<0.) 02820 { 02821 #ifdef edebug 02822 G4cerr<<"**G4QEleNucCS::GetVirtFact:K="<<K<<",nu="<<nu<<",Q2="<<Q2<<",dM="<<dM<<G4endl; 02823 #endif 02824 return 0.; 02825 } 02826 G4double lK=std::log(K); // ln(K) 02827 G4double x=1.-K/nu; // This definitin saves one div. 02828 G4double GD=1.+Q2/Q02; // Reversed nucleonic form-factor 02829 G4double b=std::exp(bp*(lK-blK0)); // b-factor 02830 G4double c=std::exp(cp*(lK-clK0)); // c-factor 02831 G4double r=.5*std::log(Q2+nu*nu)-lK; // r=.5*log((Q^2+nu^2)/K^2) 02832 G4double ef=std::exp(r*(b-c*r*r)); // exponential factor 02833 return (1.-x)*ef/GD/GD; 02834 }
G4double G4QElectronNuclearCrossSection::ThresholdEnergy | ( | G4int | Z, | |
G4int | N, | |||
G4int | PDG = 11 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 271 of file G4QElectronNuclearCrossSection.cc.
References G4cout, G4endl, G4NucleiProperties::GetNuclearMass(), and G4NucleiProperties::IsInStableTable().
Referenced by GetCrossSection().
00272 { 00273 // CHIPS - Direct GEANT 00274 //static const G4double mNeut = G4QPDGCode(2112).GetMass(); 00275 //static const G4double mProt = G4QPDGCode(2212).GetMass(); 00276 static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/MeV; 00277 static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/MeV; 00278 static const G4double mAlph = G4NucleiProperties::GetNuclearMass(4,2)/MeV; 00279 // --------- 00280 static const G4double infEn = 9.e27; 00281 00282 G4int A=Z+N; 00283 if(A<1) return infEn; 00284 else if(A==1) return 144.76; // Pi0 threshold in MeV for the proton: T>m+(m^2+2lm)/2M 00285 // CHIPS - Direct GEANT 00286 //G4double mT= G4QPDGCode(111).GetNuclMass(Z,N,0); 00287 G4double mT= 0.; 00288 if(Z&&G4NucleiProperties::IsInStableTable(A,Z)) 00289 mT = G4NucleiProperties::GetNuclearMass(A,Z)/MeV; 00290 else return 0.; // If it is not in the Table of Stable Nuclei, then the Threshold=0 00291 // --------- 00292 G4double mP= infEn; 00293 //if(Z) mP= G4QPDGCode(111).GetNuclMass(Z-1,N,0); 00294 if(A>1&&Z&&G4NucleiProperties::IsInStableTable(A-1,Z-1)) 00295 mP = G4NucleiProperties::GetNuclearMass(A-1,Z-1)/MeV; // ResNucMass for a proton 00296 G4double mN= infEn; 00297 //if(N) mN= G4QPDGCode(111).GetNuclMass(Z,N-1,0); 00298 if(A>1&&G4NucleiProperties::IsInStableTable(A-1,Z)) 00299 mN = G4NucleiProperties::GetNuclearMass(A-1,Z)/MeV; // ResNucMass for a neutron 00300 00301 G4double mA= infEn; 00302 if(A>4&&Z>1&&G4NucleiProperties::IsInStableTable(A-4,Z-2)) 00303 mA = G4NucleiProperties::GetNuclearMass(A-4.,Z-2.)/MeV;// ResNucMass for an alpha 00304 00305 G4double dP= mP +mProt - mT; 00306 G4double dN= mN +mNeut - mT; 00307 G4double dA= mA +mAlph - mT; 00308 #ifdef pdebug 00309 G4cout<<"G4QElectronNucCS::ThreshEn: mP="<<mP<<",dP="<<dP<<",mN="<<mN<<",dN="<<dN<<",mA=" 00310 <<mA<<",dA="<<dA<<",mT="<<mT<<",A="<<A<<",Z="<<Z<<G4endl; 00311 #endif 00312 if(dP<dN)dN=dP; 00313 if(dA<dN)dN=dA; 00314 return dN; 00315 }