#include <G4ElectroNuclearCrossSection.hh>
Inheritance diagram for G4ElectroNuclearCrossSection:
Public Member Functions | |
G4ElectroNuclearCrossSection (const G4String &name="ElectroNuclearXS") | |
virtual | ~G4ElectroNuclearCrossSection () |
virtual void | CrossSectionDescription (std::ostream &) const |
virtual G4bool | IsIsoApplicable (const G4DynamicParticle *aParticle, G4int, G4int, const G4Element *, const G4Material *) |
virtual G4double | GetIsoCrossSection (const G4DynamicParticle *aParticle, G4int, G4int, const G4Isotope *, const G4Element *, const G4Material *) |
G4double | GetEquivalentPhotonEnergy () |
G4double | GetVirtualFactor (G4double nu, G4double Q2) |
G4double | GetEquivalentPhotonQ2 (G4double nu) |
Definition at line 45 of file G4ElectroNuclearCrossSection.hh.
G4ElectroNuclearCrossSection::G4ElectroNuclearCrossSection | ( | const G4String & | name = "ElectroNuclearXS" |
) |
Definition at line 76 of file G4ElectroNuclearCrossSection.cc.
00077 : G4VCrossSectionDataSet(nam) 00078 {}
G4ElectroNuclearCrossSection::~G4ElectroNuclearCrossSection | ( | ) | [virtual] |
Definition at line 80 of file G4ElectroNuclearCrossSection.cc.
00081 { 00082 std::vector<G4double*>::iterator pos; 00083 for(pos=J1.begin(); pos<J1.end(); pos++) 00084 { delete [] *pos; } 00085 J1.clear(); 00086 for(pos=J2.begin(); pos<J2.end(); pos++) 00087 { delete [] *pos; } 00088 J2.clear(); 00089 for(pos=J3.begin(); pos<J3.end(); pos++) 00090 { delete [] *pos; } 00091 J3.clear(); 00092 }
void G4ElectroNuclearCrossSection::CrossSectionDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4VCrossSectionDataSet.
Definition at line 95 of file G4ElectroNuclearCrossSection.cc.
00096 { 00097 outFile << "G4ElectroNuclearCrossSection provides the total inelastic\n" 00098 << "cross section for e- and e+ interactions with nuclei. The\n" 00099 << "cross sections are retrieved from a table which is\n" 00100 << "generated using the equivalent photon approximation. In\n" 00101 << "this approximation real gammas are produced from the virtual\n" 00102 << "ones generated at the electromagnetic vertex. This cross\n" 00103 << "section set is valid for incident electrons and positrons at\n" 00104 << "all energies.\n"; 00105 }
G4double G4ElectroNuclearCrossSection::GetEquivalentPhotonEnergy | ( | ) |
Definition at line 2395 of file G4ElectroNuclearCrossSection.cc.
References G4cerr, G4cout, G4endl, and G4UniformRand.
Referenced by G4ElectroVDNuclearModel::ApplyYourself(), and G4ElectroNuclearReaction::ApplyYourself().
02396 { 02397 if(lastSig <= 0.0) { return 0.0; } // VI 02398 02399 // All constants are the copy of that from GetCrossSection funct. 02400 // => Make them general. 02401 static const G4int nE=336; // !! If you change this, change it in 02402 // GetFunctions() (*.hh) !! 02403 static const G4int mL=nE-1; 02404 static const G4double EMi=2.0612; // Minimum Energy 02405 static const G4double EMa=50000.; // Maximum Energy 02406 static const G4double lEMi=std::log(EMi); // Minimum logarithmic Energy 02407 static const G4double lEMa=std::log(EMa); // Maximum logarithmic Energy 02408 static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic step in Energy 02409 static const G4double mel=0.5109989; // Mass of electron in MeV 02410 static const G4double lmel=std::log(mel); // Log of electron mass 02411 G4double phLE = 0.; // Prototype of the std::log(nu=E_gamma) 02412 G4double Y[nE] = {0.0}; // Prepare the array for randomization 02413 02414 #ifdef debug 02415 G4cout << "G4ElectroNuclearCrossSection::GetEguPhotE:B=" 02416 << lastF<<",l=" << lastL << ",J1=" << lastJ1[lastL] 02417 << ",J2=" << lastJ2[lastL] << ",J3=" << lastJ3[lastL] << ",S=" 02418 << lastSig << ",E=" << lastE << G4endl; 02419 #endif 02420 02421 G4double lastLE=lastG+lmel; // recover std::log(eE) from the gamma (lastG) 02422 G4double dlg1=lastG+lastG-1.; 02423 G4double lgoe=lastG/lastE; 02424 for (G4int i=lastF;i<=lastL;i++) { 02425 Y[i] = dlg1*lastJ1[i]-lgoe*(lastJ2[i]+lastJ2[i]-lastJ3[i]/lastE); 02426 if(Y[i] < 0.0) { Y[i] = 0.0; } 02427 } 02428 // Tempory IF of H.P.: delete it if the *HP* err message does not 02429 // show up M.K. 02430 if(lastSig>0.99*Y[lastL] && lastL<mL && Y[lastL]<1.E-30) 02431 { 02432 G4cerr << "*HP*G4ElNucCS::GetEqPhotE:S=" << lastSig <<">" << Y[lastL] 02433 << ",l=" << lastL << ">" << mL << G4endl; 02434 if(lastSig <= 0.0) { return 0.0; } // VI 02435 } 02436 G4double ris = lastSig*G4UniformRand(); // Sig can be > Y[lastL = mL], then it 02437 // is in the funct. region 02438 #ifdef debug 02439 G4cout << "G4ElectroNuclearCrossSection::GetEquivalentPhotonEnergy: " << ris 02440 << ",Y=" << Y[lastL] << G4endl; 02441 #endif 02442 if (ris < Y[lastL]) { // Search the table 02443 G4int j = lastF; 02444 G4double Yj = Y[j]; // It must be 0 (sometimes just very small) 02445 while (ris > Yj && j < lastL) { // Associative search 02446 j++; 02447 Yj = Y[j]; // Yj is first value above ris 02448 } 02449 G4int j1 = j-1; 02450 G4double Yi = Y[j1]; // Previous value is below ris 02451 phLE = lEMi + (j1 + (ris-Yi)/(Yj-Yi) )*dlnE; 02452 #ifdef debug 02453 G4cout << "G4EleNucCS::E=" << phLE << ",l=" << lEMi << ",j=" << j << ",ris=" 02454 << ris << ",Yi=" << Yi << ",Y=" << Yj << G4endl; 02455 #endif 02456 } else { // Search with the function 02457 if (lastL < mL) G4cerr << "**G4EleNucCS::GetEfPhE:L=" << lastL << ",S=" 02458 << lastSig << ",Y=" << Y[lastL] << G4endl; 02459 G4double f = (ris-Y[lastL])/lastH; // The scaled residual value of the cross-section integral 02460 #ifdef pdebug 02461 G4cout << "G4EleNucCS::GetEfPhE:HighEnergy f=" << f << ",ris=" << ris 02462 << ",lastH=" << lastH << G4endl; 02463 #endif 02464 phLE=SolveTheEquation(f); // Solve the equation to find theLog(phE) (compare with lastLE) 02465 #ifdef pdebug 02466 G4cout << "G4EleNucCS::GetEfPhE:HighEnergy lphE=" << phLE << G4endl; 02467 #endif 02468 } 02469 02470 if (phLE>lastLE) { 02471 G4cerr << "***G4ElectroNuclearCS::GetEquPhotE:N=" << lastN << ",Z=" 02472 << lastZ << ", lpE" << phLE << ">leE" << lastLE << ",Sig=" 02473 << lastSig << ",rndSig=" << ris << ",Beg=" << lastF << ",End=" 02474 << lastL << ",Y=" << Y[lastL] << G4endl; 02475 if(lastLE<7.2) phLE=std::log(std::exp(lastLE)-.511); 02476 else phLE=7.; 02477 } 02478 return std::exp(phLE); 02479 }
Definition at line 2527 of file G4ElectroNuclearCrossSection.cc.
References G4cerr, and G4UniformRand.
Referenced by G4ElectroVDNuclearModel::ApplyYourself(), and G4ElectroNuclearReaction::ApplyYourself().
02528 { 02529 if(lastG <= 0.0 || lastE <= 0.0) { return 0.; } // VI 02530 if(lastSig <= 0.0) { return 0.0; } // VI 02531 static const G4double mel=0.5109989; // Mass of electron in MeV 02532 static const G4double mel2=mel*mel; // Squared Mass of electron in MeV 02533 G4double y=nu/lastE; // Part of energy carried by the equivalent pfoton 02534 if(y>=1.-1./(lastG+lastG)) return 0.; // The region where the method does not work 02535 G4double y2=y*y; // Squared photonic part of energy 02536 G4double ye=1.-y; // Part of energy carried by the secondary electron 02537 G4double Qi2=mel2*y2/ye; // Minimum Q2 02538 G4double Qa2=4*lastE*lastE*ye; // Maximum Q2 02539 G4double iar=Qi2/Qa2; // Q2min/Q2max ratio 02540 G4double Dy=ye+.5*y2; // D(y) function 02541 G4double Py=ye/Dy; // P(y) function 02542 G4double ePy=1.-std::exp(Py); // 1-std::exp(P(y)) part 02543 G4double Uy=Py*(1.-iar); // U(y) function 02544 G4double Fy=(ye+ye)*(1.+ye)*iar/y2; // F(y) function 02545 G4double fr=iar/(1.-ePy*iar); // Q-fraction 02546 if(Fy<=-fr) 02547 { 02548 #ifdef edebug 02549 G4cerr<<"***G4ElectroNucCrossSec::GetEquPhQ2:Fy="<<Fy<<"+fr="<<fr<<" <0"<<",iar="<<iar<<G4endl; 02550 #endif 02551 return 0.; 02552 } 02553 G4double LyQa2=std::log(Fy+fr); // L(y,Q2max) function 02554 G4bool cond=true; 02555 G4int maxTry=3; 02556 G4int cntTry=0; 02557 G4double Q2=Qi2; 02558 while(cond&&cntTry<maxTry) // The loop to avoid x>1. 02559 { 02560 G4double R=G4UniformRand(); // Random number (0,1) 02561 Q2=Qi2*(ePy+1./(std::exp(R*LyQa2-(1.-R)*Uy)-Fy)); 02562 cntTry++; 02563 cond = Q2>1878.*nu; 02564 } 02565 if(Q2<Qi2) 02566 { 02567 #ifdef edebug 02568 G4cerr<<"***G4ElectroNucCrossSec::GetEquPhQ2:Q2="<<Q2<<" < Q2min="<<Qi2<<G4endl; 02569 #endif 02570 return Qi2; 02571 } 02572 if(Q2>Qa2) 02573 { 02574 #ifdef edebug 02575 G4cerr<<"***G4ElectroNucCrossSec::GetEquPhQ2:Q2="<<Q2<<" > Q2max="<<Qi2<<G4endl; 02576 #endif 02577 return Qa2; 02578 } 02579 return Q2; 02580 }
G4double G4ElectroNuclearCrossSection::GetIsoCrossSection | ( | const G4DynamicParticle * | aParticle, | |
G4int | , | |||
G4int | , | |||
const G4Isotope * | , | |||
const G4Element * | , | |||
const G4Material * | ||||
) | [virtual] |
Reimplemented from G4VCrossSectionDataSet.
Definition at line 121 of file G4ElectroNuclearCrossSection.cc.
References G4cout, G4endl, G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4ParticleDefinition::GetPDGEncoding(), and CLHEP::detail::n.
Referenced by G4ElectroVDNuclearModel::ApplyYourself().
00125 { 00126 static const G4int nE=336; // !! If you change this, change it in GetFunctions() (*.hh) !! 00127 static const G4int mL=nE-1; 00128 static const G4double EMi=2.0612; // Minimum tabulated Energy of the Electron 00129 static const G4double EMa=50000.; // Maximum tabulated Energy of the Electron 00130 static const G4double lEMi=std::log(EMi); // Minimum tabulated logarithmic Energy of the Electron 00131 static const G4double lEMa=std::log(EMa); // Maximum tabulated logarithmic Energy of the Electron 00132 static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic step in the table for the electron Energy 00133 static const G4double alop=1./137.036/3.14159265; //coef. for the calculated functions (Ee>50000.) 00134 static const G4double mel=0.5109989; // Mass of the electron in MeV 00135 static const G4double lmel=std::log(mel); // Log of the electron mass 00136 // *** Begin of the Associative memory for acceleration of the cross section calculations 00137 static std::vector <G4int> colN; // Vector of N for calculated nucleus (isotop) 00138 static std::vector <G4int> colZ; // Vector of Z for calculated nucleus (isotop) 00139 static std::vector <G4int> colF; // Vector of Last StartPosition in the Ji-function tables 00140 static std::vector <G4double> colTH; // Vector of the energy thresholds for the eA->eX reactions 00141 static std::vector <G4double> colH; // Vector of HighEnergyCoefficients (functional calculations) 00142 // *** End of Static Definitions (Associative Memory) *** 00143 00144 const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the electron 00145 const G4int targetAtomicNumber = AA; 00146 const G4int targZ = ZZ; 00147 const G4int targN = targetAtomicNumber-targZ; // @@ Get isotops (can change initial A) 00148 if (Energy<=EMi) return 0.; // Energy is below the minimum energy in the table 00149 00150 G4int PDG=aPart->GetDefinition()->GetPDGEncoding(); 00151 if (PDG == 11 || PDG == -11) // @@ Now only for electrons, but can be fo muons 00152 { 00153 G4double A = targN + targZ; // New A (can differ from G4double targetAtomicNumber) 00154 if(targN!=lastN || targZ!=lastZ) // This nucleus was not the last used isotop 00155 { 00156 lastE = 0.; // New history in the electron Energy 00157 lastG = 0.; // New history in the photon Energy 00158 lastN = targN; // The last N of calculated nucleus 00159 lastZ = targZ; // The last Z of calculated nucleus 00160 G4int n=colN.size(); // Size of the Associative Memory DB in the heap 00161 G4bool in=false; // "Found in AMDB" flag 00162 if(n) for(G4int i=0; i<n; i++) if(colN[i]==targN && colZ[i]==targZ) // Calculated nuclei 00163 { // The nucleus is found in AMDB 00164 in=true; // Rais the "Found in AMDB" flag 00165 lastTH =colTH[i]; // Last Energy threshold (A-dependent) 00166 lastF =colF[i]; // Last ZeroPosition in the J-functions 00167 lastH =colH[i]; // Last High Energy Coefficient (A-dependent) 00168 lastJ1 =J1[i]; // Pointer to the prepared J1 function 00169 lastJ2 =J2[i]; // Pointer to the prepared J2 function 00170 lastJ3 =J3[i]; // Pointer to the prepared J3 function 00171 } 00172 if(!in) // This nucleus has not been calculated previously 00173 { 00174 lastJ1 = new G4double[nE]; // Allocate memory for the new J1 function 00175 lastJ2 = new G4double[nE]; // Allocate memory for the new J2 function 00176 lastJ3 = new G4double[nE]; // Allocate memory for the new J3 function 00177 lastF = GetFunctions(A,lastJ1,lastJ2,lastJ3); // new ZeroPos and filling of J-functions 00178 lastH = alop*A*(1.-.072*std::log(A));// corresponds to lastSP from G4PhotonuclearCrossSection 00179 lastTH = ThresholdEnergy(targZ, targN); // The last Threshold Energy 00180 #ifdef pdebug 00181 G4cout<<"G4ElNucCS::GetCrossSection: lastH="<<lastH<<",A="<<A<<G4endl; 00182 #endif 00183 colN.push_back(targN); 00184 colZ.push_back(targZ); 00185 colF.push_back(lastF); 00186 J1.push_back(lastJ1); 00187 J2.push_back(lastJ2); 00188 J3.push_back(lastJ3); 00189 colH.push_back(lastH); 00190 colTH.push_back(lastTH); 00191 } // End of creation of the new set of parameters 00192 } // End of parameters udate 00193 00194 // else if(std::abs((lastE-Energy)/Energy)<.001) return lastSig*millibarn; // Don't calc. same CS twice 00195 else if(lastE == Energy) return lastSig*millibarn; // Don't calc. same CS twice 00196 // ============================== NOW Calculate the Cross Section ========================== 00197 lastE=Energy; // lastE - the electron energy 00198 if (Energy<=lastTH) // Once more check that the eE is higher than the ThreshE 00199 { 00200 lastSig=0.; 00201 return 0.; 00202 } 00203 G4double lE=std::log(Energy); // std::log(eE) (it is necessary at this point for the fit) 00204 lastG=lE-lmel; // Gamma of the electron (used to recover std::log(eE)) 00205 G4double dlg1=lastG+lastG-1.; 00206 G4double lgoe=lastG/lastE; 00207 if(lE<lEMa) // Linear fit is made explicitly to fix the last bin for the randomization 00208 { 00209 G4double shift=(lE-lEMi)/dlnE; 00210 G4int blast=static_cast<int>(shift); 00211 if(blast<0) blast=0; 00212 if(blast>=mL) blast=mL-1; 00213 shift-=blast; 00214 lastL=blast+1; 00215 G4double YNi=dlg1*lastJ1[blast]-lgoe*(lastJ2[blast]+lastJ2[blast]-lastJ3[blast]/lastE); 00216 G4double YNj=dlg1*lastJ1[lastL]-lgoe*(lastJ2[lastL]+lastJ2[lastL]-lastJ3[lastL]/lastE); 00217 lastSig= YNi+shift*(YNj-YNi); 00218 if(lastSig>YNj)lastSig=YNj; 00219 #ifdef pdebug 00220 G4cout<<"G4ElNucCS::GetCS:S="<<lastSig<<",E="<<lE<<",Yi="<<YNi<<",Yj="<<YNj<<",M="<<lEMa<<G4endl; 00221 G4cout<<"G4EN::GCS:s="<<shift<<",Jb="<<lastJ1[blast]<<",J="<<lastJ1[lastL]<<",b="<<blast<<G4endl; 00222 #endif 00223 } 00224 else 00225 { 00226 lastL=mL; 00227 G4double term1=lastJ1[mL]+lastH*HighEnergyJ1(lE); 00228 G4double term2=lastJ2[mL]+lastH*HighEnergyJ2(lE); 00229 G4double term3=lastJ3[mL]+lastH*HighEnergyJ3(lE); 00230 lastSig=dlg1*term1-lgoe*(term2+term2-term3/lastE); 00231 #ifdef pdebug 00232 G4cout<<"G4ElNucCS::GetCrossSec:S="<<lastSig<<",lE="<<lE<<",J1="<<lastH*HighEnergyJ1(lE)<<",Pm=" 00233 <<lastJ1[mL]<<",Fm="<<lastJ2[mL]<<",Fh="<<lastH*HighEnergyJ2(lE)<<",EM="<<lEMa<<G4endl; 00234 #endif 00235 } 00236 } // End of "sigma" calculation 00237 else return 0.; 00238 if(lastSig<0.) lastSig = 0.; 00239 lastE=Energy; 00240 return lastSig*millibarn; 00241 }
Definition at line 2582 of file G4ElectroNuclearCrossSection.cc.
References G4cerr.
Referenced by G4ElectroNuclearReaction::ApplyYourself().
02583 { 02584 if(nu <= 0.0 || Q2 <= 0.0) { return 0.0; } 02585 static const G4double dM=938.27+939.57; // Mean double nucleon mass = m_n+m_p (@@ no binding) 02586 static const G4double Q0=843.; // Coefficient of the dipole nucleonic form-factor 02587 static const G4double Q02=Q0*Q0; // Squared coefficient of the dipole nucleonic form-factor 02588 static const G4double blK0=std::log(185.); // Coefficient of the b-function 02589 static const G4double bp=0.85; // Power of the b-function 02590 static const G4double clK0=std::log(1390.); // Coefficient of the c-function 02591 static const G4double cp=3.; // Power of the c-function 02592 //G4double x=Q2/dM/nu; // Direct x definition 02593 G4double K=nu-Q2/dM; // K=nu*(1-x) 02594 if(K <= 0.) // VI 02595 { 02596 #ifdef edebug 02597 G4cerr<<"**G4ElectroNucCrossSec::GetVirtFact:K="<<K<<",nu="<<nu<<",Q2="<<Q2<<",dM="<<dM<<G4endl; 02598 #endif 02599 return 0.; 02600 } 02601 G4double lK=std::log(K); // ln(K) 02602 G4double x=1.-K/nu; // This definitin saves one div. 02603 G4double GD=1.+Q2/Q02; // Reversed nucleonic form-factor 02604 G4double b=std::exp(bp*(lK-blK0)); // b-factor 02605 G4double c=std::exp(cp*(lK-clK0)); // c-factor 02606 G4double r=.5*std::log(Q2+nu*nu)-lK; // r=.5*std::log((Q^2+nu^2)/K^2) 02607 G4double ef=std::exp(r*(b-c*r*r)); // exponential factor 02608 return (1.-x)*ef/GD/GD; 02609 }
G4bool G4ElectroNuclearCrossSection::IsIsoApplicable | ( | const G4DynamicParticle * | aParticle, | |
G4int | , | |||
G4int | , | |||
const G4Element * | , | |||
const G4Material * | ||||
) | [virtual] |
Reimplemented from G4VCrossSectionDataSet.
Definition at line 108 of file G4ElectroNuclearCrossSection.cc.
References G4Electron::ElectronDefinition(), G4DynamicParticle::GetDefinition(), and G4Positron::PositronDefinition().
00111 { 00112 G4bool result = false; 00113 if (aParticle->GetDefinition() == G4Electron::ElectronDefinition()) 00114 result = true; 00115 if (aParticle->GetDefinition() == G4Positron::PositronDefinition()) 00116 result = true; 00117 return result; 00118 }