G4QIonIonCrossSection Class Reference

#include <G4QIonIonCrossSection.hh>

Inheritance diagram for G4QIonIonCrossSection:

G4VQCrossSection

Public Member Functions

 ~G4QIonIonCrossSection ()
virtual G4double GetCrossSection (G4bool fCS, G4double pMom, G4int Z, G4int N, G4int PDG)
G4double CalculateCrossSection (G4bool fCS, G4int F, G4int I, G4int PDG, G4int tZ, G4int tN, G4double Momentum)
G4double ThresholdMomentum (G4int pZ, G4int pN, G4int tZ, G4int tN)

Static Public Member Functions

static G4VQCrossSectionGetPointer ()

Protected Member Functions

 G4QIonIonCrossSection ()

Detailed Description

Definition at line 59 of file G4QIonIonCrossSection.hh.


Constructor & Destructor Documentation

G4QIonIonCrossSection::G4QIonIonCrossSection (  )  [inline, protected]

Definition at line 63 of file G4QIonIonCrossSection.hh.

00063 {}

G4QIonIonCrossSection::~G4QIonIonCrossSection (  )  [inline]

Definition at line 67 of file G4QIonIonCrossSection.hh.

00067 {}


Member Function Documentation

G4double G4QIonIonCrossSection::CalculateCrossSection ( G4bool  fCS,
G4int  F,
G4int  I,
G4int  PDG,
G4int  tZ,
G4int  tN,
G4double  Momentum 
) [virtual]

Implements G4VQCrossSection.

Definition at line 252 of file G4QIonIonCrossSection.cc.

References G4VQCrossSection::EquLinearFit(), G4cerr, G4cout, and G4endl.

Referenced by GetCrossSection().

00254 {
00255   //static const G4double third=1./3.; // power for A^P->R conversion [R=1.1*A^(1/3)]
00256   //static const G4double conv=38.; // coeff. R2->sig=c*(pR+tR)^2, c=pi*10(mb/fm^2)*1.21
00257   // If change the following, please change in ::GetFunctions:
00258   static const G4double THmin=0.;  // @@ start from threshold (?) minimum Energy Threshold
00259   static const G4double dP=10.;    // step for the LEN table
00260   static const G4int    nL=100;    // A#of LENesonance points in E (each MeV from 2 to 106)
00261   static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part
00262   static const G4double Pmax=300000.;   // maxE for the HighE part
00263   static const G4int    nH=100;         // A#of HResonance points in lnE
00264   static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part
00265   static const G4double malP=std::log(Pmax); // High logarithm energy (each 2.75 percent)
00266   static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part
00267   //
00268   // Associative memory for acceleration
00269   static std::vector <G4double*> LENI;   // Vector of pointers: LowEnIneIonIonCrossSection
00270   static std::vector <G4double*> HENI;   // Vector of pointers: HighEnIneIonIonCrossSection
00271   static std::vector <G4double*> LENE;   // Vector of pointers: LowEnElaIonIonCrossSection
00272   static std::vector <G4double*> HENE;   // Vector of pointers: HighEnElaIonIonCrossSection
00273 #ifdef debug
00274   G4cout<<"G4QIonIonCrossSection::CalcCS: Z="<<tZ<<", N="<<tN<<", P="<<TotMom<<G4endl;
00275 #endif
00276   G4int dPDG=pPDG/10;                // 10SZZZAAA
00277   G4int zPDG=dPDG/1000;              // 10SZZZ (?)
00278   G4int zA=dPDG%1000;                // proj A
00279   G4int pZ=zPDG%1000;                // proj Z (?)
00280   G4int pN=zA-pZ;                    // proj N (?)
00281   G4double Momentum=TotMom/zA;       // Momentum per nucleon
00282   if (Momentum<THmin) return 0.;     // @@ This can be dangerouse for the heaviest nuc.!
00283   G4double sigma=0.;
00284   if(F&&I) sigma=0.;                 // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
00285   G4double tA=tN+tZ;                 // Target weight
00286   G4double pA=zA;                    // Projectile weight
00287   if(F<=0)                           // This isotope was not the last used isotop
00288   {
00289     if(F<0 || !XS)                   // This isotope was found in DAMDB or Elast =>RETRIEVE
00290     {
00291       lastLENI=LENI[I];              // Pointer to Low Energy inelastic cross sections
00292       lastHENI=HENI[I];              // Pointer to High Energy inelastic cross sections
00293       lastLENE=LENE[I];              // Pointer to Low Energy inelastic cross sections
00294       lastHENE=HENE[I];              // Pointer to High Energy inelastic cross sections
00295     }
00296     else                             // This isotope wasn't calculated previously => CREATE
00297     {
00298       lastLENI = new G4double[nL];   // Allocate memory for the new LEN cross sections
00299       lastHENI = new G4double[nH];   // Allocate memory for the new HEN cross sections
00300       lastLENE = new G4double[nL];   // Allocate memory for the new LEN cross sections
00301       lastHENE = new G4double[nH];   // Allocate memory for the new HEN cross sections
00302       G4int er=GetFunctions(pZ,pN,tZ,tN,lastLENI,lastHENI,lastLENE,lastHENE);
00303       if(er<1) G4cerr<<"*W*G4QIonIonCroSec::CalcCrossSection: pA="<<tA<<",tA="<<tA<<G4endl;
00304 #ifdef debug
00305       G4cout<<"G4QIonIonCrossSection::CalcCS: GetFunctions er="<<er<<",pA="<<pA<<",tA="<<tA
00306             <<G4endl;
00307 #endif
00308       // *** The synchronization check ***
00309       G4int sync=LENI.size();
00310       if(sync!=I) G4cout<<"*W*G4IonIonCrossSec::CalcCrossSect:Sync="<<sync<<"#"<<I<<G4endl;
00311       LENI.push_back(lastLENI);      // added LEN Inelastic
00312       HENI.push_back(lastHENI);      // added HEN Inelastic
00313       LENE.push_back(lastLENE);      // added LEN Elastic
00314       HENE.push_back(lastHENE);      // added HEN Elastic
00315     } // End of creation of the new set of parameters
00316   } // End of parameters udate
00317   // =------------= NOW the Magic Formula =--------------------------=
00318   if (Momentum<lastTH) return 0.;    // It must be already checked in the interface class
00319   else if (Momentum<Pmin)            // LEN region (approximated in E, not in lnE)
00320   {
00321 #ifdef debug
00322     G4cout<<"G4QIICS::CalCS:p="<<pA<<",t="<<tA<<",n="<<nL<<",T="<<THmin<<",d="<<dP<<G4endl;
00323 #endif
00324     if(tA<1. || pA<1.)
00325     {
00326       G4cout<<"-Warning-G4QIICS::CalcCS: pA="<<pA<<" or tA="<<tA<<" aren't nuclei"<<G4endl;
00327       sigma=0.;
00328     }
00329     else
00330     {
00331       G4double dPp=dP*pA;
00332       if(XS) sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENI);
00333       else   sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENE);
00334     }
00335 #ifdef debugn
00336     if(sigma<0.) G4cout<<"-Warning-G4QIICS::CalcCS:pA="<<pA<<",tA="<<tA<<",XS="<<XS<<",P="
00337                        <<Momentum<<", Th="<<THmin<<", dP="<<dP<<G4endl;
00338 #endif
00339   }
00340   else if (Momentum<Pmax*pA)                     // High Energy region
00341   {
00342     G4double lP=std::log(Momentum);
00343 #ifdef debug
00344     G4cout<<"G4QIonIonCS::CalcCS:before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl;
00345 #endif
00346     if(tA<1. || pA<1.)
00347     {
00348       G4cout<<"-Warning-G4QIICS::CalCS:pA="<<pA<<" or tA="<<tA<<" aren't composit"<<G4endl;
00349       sigma=0.;
00350     }
00351     else
00352     {
00353       G4double milPp=milP+std::log(pA);
00354       if(XS) sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENI);
00355       else   sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENE);
00356     }
00357   }
00358   else                                      // UltraHighE region (not frequent)
00359   {
00360     std::pair<G4double, G4double> inelel = CalculateXS(pZ, pN, tZ, tN, Momentum);
00361     if(XS) sigma=inelel.first;
00362     else   sigma=inelel.second;
00363   }
00364 #ifdef debug
00365   G4cout<<"G4IonIonCrossSection::CalculateCrossSection: sigma="<<sigma<<G4endl;
00366 #endif
00367   if(sigma<0.) return 0.;
00368   return sigma;
00369 }

G4double G4QIonIonCrossSection::GetCrossSection ( G4bool  fCS,
G4double  pMom,
G4int  Z,
G4int  N,
G4int  PDG 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 73 of file G4QIonIonCrossSection.cc.

References CalculateCrossSection(), G4cout, G4endl, G4VQCrossSection::ThresholdEnergy(), and G4VQCrossSection::tolerance.

00075 {
00076   static G4int j;                      // A#0f records found in DB for this projectile
00077   static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
00078   static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
00079   static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
00080   static std::vector <G4double> colP;  // Vector of last momenta for the reaction
00081   static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
00082   static std::vector <G4double> colICS;// Vector of last inelastic cross-sections
00083   static std::vector <G4double> colECS;// Vector of last elastic cross-sections
00084   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
00085 #ifdef pdebug
00086   G4cout<<"G4QIICS::GetCS:>>> f="<<fCS<<", Z="<<tZ<<"("<<lastZ<<"), N="<<tN<<"("<<lastN
00087         <<"),PDG="<<pPDG<<"("<<lastPDG<<"), p="<<pMom<<"("<<lastTH<<")"<<",Sz="
00088         <<colN.size()<<G4endl;
00089 #endif
00090   if(!pPDG)
00091   {
00092 #ifdef pdebug
00093     G4cout<<"G4QIonIonCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00094 #endif
00095     return 0.;                         // projectile PDG=0 is a mistake (?!) @@
00096   }
00097   G4bool in=false;                     // By default the isotope must be found in the AMDB
00098   if(tN!=lastN || tZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
00099   {
00100     in = false;                        // By default the isotope haven't be found in AMDB  
00101     lastP   = 0.;                      // New momentum history (nothing to compare with)
00102     lastPDG = pPDG;                    // The last PDG of the projectile
00103     lastN   = tN;                      // The last N of the calculated nucleus
00104     lastZ   = tZ;                      // The last Z of the calculated nucleus
00105     lastI   = colN.size();             // Size of the Associative Memory DB in the heap
00106     j  = 0;                            // A#0f records found in DB for this projectile
00107 #ifdef pdebug
00108     G4cout<<"G4QIICS::GetCS:FindI="<<lastI<<",pPDG="<<pPDG<<",tN="<<tN<<",tZ="<<tZ<<G4endl;
00109 #endif
00110     if(lastI) for(G4int i=0; i<lastI; i++) // Loop over all DB
00111     {                                  // The nucleus with projPDG is found in AMDB
00112 #ifdef pdebug
00113       G4cout<<"G4QII::GCS:P="<<colPDG[i]<<",N="<<colN[i]<<",Z="<<colZ[i]<<",j="<<j<<G4endl;
00114 #endif
00115       if(colPDG[i]==pPDG && colN[i]==tN && colZ[i]==tZ)
00116       {
00117         lastI=i;
00118         lastTH =colTH[i];                // Last THreshold (A-dependent)
00119 #ifdef pdebug
00120         G4cout<<"G4QIICS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00121 #endif
00122         if(pMom<=lastTH)
00123         {
00124 #ifdef pdebug
00125           G4cout<<"G4QIICS::GetCS:Found P="<<pMom<<"<Threshold="<<lastTH<<"->XS=0"<<G4endl;
00126 #endif
00127           return 0.;                     // Energy is below the Threshold value
00128         }
00129         lastP  =colP [i];                // Last Momentum  (A-dependent)
00130         lastICS=colICS[i];               // Last Inelastic Cross-Section (A-dependent)
00131         lastECS=colECS[i];               // Last Elastic Cross-Section (A-dependent)
00132         if(std::fabs(lastP/pMom-1.)<tolerance)
00133         {
00134 #ifdef pdebug
00135           G4cout<<"G4QIonIonCS::GetCS:P="<<pMom<<",InXS="<<lastICS*millibarn<<",ElXS="
00136                 <<lastECS*millibarn<<G4endl;
00137 #endif
00138           CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only
00139           if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
00140           return         lastECS*millibarn;     // Use theLastElasticCS
00141         }
00142         in = true;                       // This is the case when the isotop is found in DB
00143         // Momentum pMom is in IU ! @@ Units
00144 #ifdef pdebug
00145         G4cout<<"G4QIICS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00146 #endif
00147         lastICS=CalculateCrossSection( true,-1,j,lastPDG,lastZ,lastN,pMom);// read & update
00148         lastECS=CalculateCrossSection(false,-1,j,lastPDG,lastZ,lastN,pMom);// read & update
00149 #ifdef pdebug
00150         G4cout<<"G4QIonIonCS::GetCS:=>New(inDB) InCS="<<lastICS<<",ElCS="<<lastECS<<G4endl;
00151 #endif
00152         if((lastICS<=0. || lastECS<=0.) && pMom>lastTH) // Correct the threshold
00153         {
00154 #ifdef pdebug
00155           G4cout<<"G4QIonIonCS::GetCS:New,T="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00156 #endif
00157           lastTH=pMom;
00158         }
00159         break;                           // Go out of the LOOP
00160       }
00161 #ifdef pdebug
00162       G4cout<<"--->G4QIonIonCrossSec::GetCrosSec: pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
00163             <<",Z["<<i<<"]="<<colZ[i]<<",PDG="<<colPDG[i]<<G4endl;
00164 #endif
00165       j++;                             // Increment a#0f records found in DB for this pPDG
00166     }
00167     if(!in)                            // This nucleus has not been calculated previously
00168     {
00169 #ifdef pdebug
00170       G4cout<<"G4QIICS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
00171 #endif
00173       lastICS=CalculateCrossSection(true ,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
00174       lastECS=CalculateCrossSection(false,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
00175       if(lastICS<=0. || lastECS<=0.)
00176       {
00177         lastTH = ThresholdEnergy(tZ, tN); // Threshold Energy=Mom=0 which is now the last
00178 #ifdef pdebug
00179         G4cout<<"G4QIonIonCrossSect::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
00180 #endif
00181         if(pMom>lastTH)
00182         {
00183 #ifdef pdebug
00184           G4cout<<"G4QIonIonCS::GetCS:1-st,P="<<pMom<<">Thresh="<<lastTH<<"->XS=0"<<G4endl;
00185 #endif
00186           lastTH=pMom;
00187         }
00188       }
00189 #ifdef pdebug
00190       G4cout<<"G4QIICS::GetCS: *New* ICS="<<lastICS<<", ECS="<<lastECS<<",N="<<lastN<<",Z="
00191             <<lastZ<<G4endl;
00192 #endif
00193       colN.push_back(tN);
00194       colZ.push_back(tZ);
00195       colPDG.push_back(pPDG);
00196       colP.push_back(pMom);
00197       colTH.push_back(lastTH);
00198       colICS.push_back(lastICS);
00199       colECS.push_back(lastECS);
00200 #ifdef pdebug
00201       G4cout<<"G4QIICS::GetCS:*1st*, P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn
00202             <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
00203 #endif
00204       if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
00205       return         lastECS*millibarn;     // Use theLastElasticCS
00206     } // End of creation of the new set of parameters
00207     else
00208     {
00209 #ifdef pdebug
00210       G4cout<<"G4QIICS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00211 #endif
00212       colP[lastI]=pMom;
00213       colPDG[lastI]=pPDG;
00214       colICS[lastI]=lastICS;
00215       colECS[lastI]=lastECS;
00216     }
00217   } // End of parameters udate
00218   else if(pMom<=lastTH)
00219   {
00220 #ifdef pdebug
00221     G4cout<<"G4QIICS::GetCS: Current T="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00222 #endif
00223     return 0.;                         // Momentum is below the Threshold Value -> CS=0
00224   }
00225   else if(std::fabs(lastP/pMom-1.)<tolerance)
00226   {
00227 #ifdef pdebug
00228     G4cout<<"G4QIICS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", InCS="<<lastICS*millibarn
00229           <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
00230 #endif
00231     if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
00232     return         lastECS*millibarn;     // Use theLastElasticCS
00233   }
00234   else
00235   {
00236 #ifdef pdebug
00237     G4cout<<"G4QIICS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00238 #endif
00239     lastICS=CalculateCrossSection( true,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00240     lastECS=CalculateCrossSection(false,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00241     lastP=pMom;
00242   }
00243 #ifdef pdebug
00244   G4cout<<"G4QIICS::GetCroSec:*End*,P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn<<", ElCS="
00245         <<lastECS*millibarn<<"(mb)"<<G4endl;
00246 #endif
00247     if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
00248     return         lastECS*millibarn;     // Use theLastElasticCS
00249 }

G4VQCrossSection * G4QIonIonCrossSection::GetPointer (  )  [static]

Definition at line 65 of file G4QIonIonCrossSection.cc.

Referenced by G4QLowEnergy::GetMeanFreePath(), G4QIonIonElastic::GetMeanFreePath(), and G4QIonIonElastic::PostStepDoIt().

00066 {
00067   static G4QIonIonCrossSection theCrossSection; //**Static body of Cross Section**
00068   return &theCrossSection;
00069 }

G4double G4QIonIonCrossSection::ThresholdMomentum ( G4int  pZ,
G4int  pN,
G4int  tZ,
G4int  tN 
)

Definition at line 519 of file G4QIonIonCrossSection.cc.

00520 {
00521   static const G4double third=1./3.;
00522   static const G4double pM = G4QPDGCode(2212).GetMass(); // Proton mass in MeV
00523   static const G4double tpM= pM+pM;       // Doubled proton mass (MeV)
00524   if(pZ<.99 || pN<0. || tZ<.99 || tN<0.) return 0.;
00525   G4double tA=tZ+tN;
00526   G4double pA=pZ+pN;
00527   //G4double dE=1.263*tZ/(1.+std::pow(tA,third));
00528   G4double dE=pZ*tZ/(std::pow(pA,third)+std::pow(tA,third))/pA; // dE/pA (per projNucleon)
00529   return std::sqrt(dE*(tpM+dE));
00530 }


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