#include <G4QIonIonCrossSection.hh>
Inheritance diagram for G4QIonIonCrossSection:
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 G4VQCrossSection * | GetPointer () |
Protected Member Functions | |
G4QIonIonCrossSection () |
Definition at line 59 of file G4QIonIonCrossSection.hh.
G4QIonIonCrossSection::G4QIonIonCrossSection | ( | ) | [inline, protected] |
G4QIonIonCrossSection::~G4QIonIonCrossSection | ( | ) | [inline] |
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 }
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 }