G4ElectroNuclearCrossSection Class Reference

#include <G4ElectroNuclearCrossSection.hh>

Inheritance diagram for G4ElectroNuclearCrossSection:

G4VCrossSectionDataSet

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)

Detailed Description

Definition at line 45 of file G4ElectroNuclearCrossSection.hh.


Constructor & Destructor Documentation

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 }


Member Function Documentation

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 }

G4double G4ElectroNuclearCrossSection::GetEquivalentPhotonQ2 ( G4double  nu  ) 

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 }

G4double G4ElectroNuclearCrossSection::GetVirtualFactor ( G4double  nu,
G4double  Q2 
)

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 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:51 2013 for Geant4 by  doxygen 1.4.7