#include <G4QPionPlusNuclearCrossSection.hh>
Inheritance diagram for G4QPionPlusNuclearCrossSection:
Public Member Functions | |
~G4QPionPlusNuclearCrossSection () | |
virtual G4double | GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=211) |
G4double | CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum) |
Static Public Member Functions | |
static G4VQCrossSection * | GetPointer () |
Protected Member Functions | |
G4QPionPlusNuclearCrossSection () |
Definition at line 51 of file G4QPionPlusNuclearCrossSection.hh.
G4QPionPlusNuclearCrossSection::G4QPionPlusNuclearCrossSection | ( | ) | [inline, protected] |
G4QPionPlusNuclearCrossSection::~G4QPionPlusNuclearCrossSection | ( | ) |
G4double G4QPionPlusNuclearCrossSection::CalculateCrossSection | ( | G4bool | CS, | |
G4int | F, | |||
G4int | I, | |||
G4int | PDG, | |||
G4int | Z, | |||
G4int | N, | |||
G4double | Momentum | |||
) | [virtual] |
Implements G4VQCrossSection.
Definition at line 231 of file G4QPionPlusNuclearCrossSection.cc.
References G4VQCrossSection::EquLinearFit(), G4cerr, G4cout, G4endl, and CLHEP::detail::n.
Referenced by GetCrossSection().
00233 { 00234 static const G4double THmin=27.; // default minimum Momentum (MeV/c) Threshold 00235 static const G4double THmiG=THmin*.001; // minimum Momentum (GeV/c) Threshold 00236 static const G4double dP=10.; // step for the LEN (Low ENergy) table MeV/c 00237 static const G4double dPG=dP*.001; // step for the LEN (Low ENergy) table GeV/c 00238 static const G4int nL=105; // A#of LEN points in E (step 10 MeV/c) 00239 static const G4double Pmin=THmin+(nL-1)*dP; // minP for the HighE part with safety 00240 static const G4double Pmax=227000.; // maxP for the HEN (High ENergy) part 227 GeV 00241 static const G4int nH=224; // A#of HEN points in lnE 00242 static const G4double milP=std::log(Pmin);// Low logarithm energy for the HEN part 00243 static const G4double malP=std::log(Pmax);// High logarithm energy (each 2.75 percent) 00244 static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HEN part 00245 static const G4double milPG=std::log(.001*Pmin);// Low logarithmEnergy for HEN part GeV/c 00246 #ifdef debug 00247 G4cout<<"G4QPipNuCS::CalCS:N="<<targN<<",Z="<<targZ<<",P="<<Momentum<<">"<<THmin<<G4endl; 00248 #endif 00249 G4double sigma=0.; 00250 if(F&&I) sigma=0.; // @@ *!* Fake line *!* to use F & I !!!Temporary!!! 00251 //G4double A=targN+targZ; // A of the target 00252 #ifdef debug 00253 G4cout<<"G4QPipNucCS::CalCS: A="<<A<<", F="<<F<<",I="<<I<<",nL="<<nL<<",nH="<<nH<<G4endl; 00254 #endif 00255 if(F<=0) // This isotope was not the last used isotop 00256 { 00257 if(F<0) // This isotope was found in DAMDB =-----=> RETRIEVE 00258 { 00259 G4int sync=LEN->size(); 00260 if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl; 00261 lastLEN=(*LEN)[I]; // Pointer to prepared LowEnergy cross sections 00262 lastHEN=(*HEN)[I]; // Pointer to prepared High Energy cross sections 00263 } 00264 else // This isotope wasn't calculated before => CREATE 00265 { 00266 lastLEN = new G4double[nL]; // Allocate memory for the new LEN cross sections 00267 lastHEN = new G4double[nH]; // Allocate memory for the new HEN cross sections 00268 // --- Instead of making a separate function --- 00269 G4double P=THmiG; // Table threshold in GeV/c 00270 for(G4int n=0; n<nL; n++) 00271 { 00272 lastLEN[n] = CrossSectionLin(targZ, targN, P); 00273 P+=dPG; 00274 } 00275 G4double lP=milPG; 00276 for(G4int n=0; n<nH; n++) 00277 { 00278 lastHEN[n] = CrossSectionLog(targZ, targN, lP); 00279 lP+=dlP; 00280 } 00281 #ifdef debug 00282 G4cout<<"-*->G4QPipNucCS::CalcCS:Tab for Z="<<targZ<<", N="<<targN<<",I="<<I<<G4endl; 00283 #endif 00284 // --- End of possible separate function 00285 // *** The synchronization check *** 00286 G4int sync=LEN->size(); 00287 if(sync!=I) 00288 { 00289 G4cerr<<"***G4QPiMinusNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ 00290 <<", N="<<targN<<", F="<<F<<G4endl; 00291 //G4Exception("G4PiMinusNuclearCS::CalculateCS:","39",FatalException,"DBoverflow"); 00292 } 00293 LEN->push_back(lastLEN); // remember the Low Energy Table 00294 HEN->push_back(lastHEN); // remember the High Energy Table 00295 } // End of creation of the new set of parameters 00296 } // End of parameters udate 00297 // =-----------------= NOW the Magic Formula =-------------------------= 00298 #ifdef debug 00299 G4cout<<"G4QPipNCS::CalcCS:lTH="<<lastTH<<",Pi="<<Pmin<<",dP="<<dP<<",dlP="<<dlP<<G4endl; 00300 #endif 00301 if (Momentum<lastTH) return 0.; // It must be already checked in the interface class 00302 else if (Momentum<Pmin) // High Energy region 00303 { 00304 #ifdef debug 00305 G4cout<<"G4QPipNCS::CalcCS:bLEN nL="<<nL<<",TH="<<THmin<<",dP="<<dP<<G4endl; 00306 #endif 00307 sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN); 00308 #ifdef debugn 00309 if(sigma<0.) 00310 G4cout<<"G4QPipNuCS::CalCS: E="<<Momentum<<",T="<<THmin<<",dP="<<dP<<G4endl; 00311 #endif 00312 } 00313 else if (Momentum<Pmax) // High Energy region 00314 { 00315 G4double lP=std::log(Momentum); 00316 #ifdef debug 00317 G4cout<<"G4QPipNucCS::CalcCS: before HEN nH="<<nH<<", iE="<<milP<<",dlP="<<dlP<<G4endl; 00318 #endif 00319 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN); 00320 } 00321 else // UHE region (calculation, not frequent) 00322 { 00323 G4double P=0.001*Momentum; // Approximation formula is for P in GeV/c 00324 sigma=CrossSectionFormula(targZ, targN, P, std::log(P)); 00325 } 00326 #ifdef debug 00327 G4cout<<"G4QPionPlusNuclearCrossSection::CalcCS: CS="<<sigma<<G4endl; 00328 #endif 00329 if(sigma<0.) return 0.; 00330 return sigma; 00331 }
G4double G4QPionPlusNuclearCrossSection::GetCrossSection | ( | G4bool | fCS, | |
G4double | pMom, | |||
G4int | tgZ, | |||
G4int | tgN, | |||
G4int | pPDG = 211 | |||
) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 81 of file G4QPionPlusNuclearCrossSection.cc.
References CalculateCrossSection(), G4cout, G4endl, G4VQCrossSection::ThresholdEnergy(), and G4VQCrossSection::tolerance.
00083 { 00084 //A.R.23-Oct-2012 Shadowed variable static G4double tolerance=0.001; // Tolerance (0.1%) to consider as "the same mom" 00085 static G4int j; // A#0f Z/N-records already tested in AMDB 00086 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) 00087 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) 00088 static std::vector <G4double> colP; // Vector of last momenta for the reaction 00089 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction 00090 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction 00091 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** 00092 #ifdef debug 00093 G4cout<<"G4QPipCS::GetCS:>>>f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN 00094 <<"("<<lastN<<"),PDG=211, thresh="<<lastTH<<",Sz="<<colN.size()<<G4endl; 00095 #endif 00096 if(PDG!=211) G4cout<<"-Warning-G4QPionPlusCS::GetCS:**Not a PiPlus**,PDG="<<PDG<<G4endl; 00097 G4bool in=false; // By default the isotope must be found in the AMDB 00098 if(tgN!=lastN || tgZ!=lastZ) // 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 lastN = tgN; // The last N of the calculated nucleus 00103 lastZ = tgZ; // The last Z of the calculated nucleus 00104 lastI = colN.size(); // Size of the Associative Memory DB in the heap 00105 j = 0; // A#0f records found in DB for this projectile 00106 #ifdef debug 00107 G4cout<<"G4QPipCS::GetCS: the amount of records in the AMDB lastI="<<lastI<<G4endl; 00108 #endif 00109 if(lastI) for(G4int i=0; i<lastI; i++) // AMDB exists, try to find the (Z,N) isotope 00110 { 00111 if(colN[i]==tgN && colZ[i]==tgZ) // Try the record "i" in the AMDB 00112 { 00113 lastI=i; // Remember the index for future fast/last use 00114 lastTH =colTH[i]; // The last THreshold (A-dependent) 00115 #ifdef debug 00116 G4cout<<"G4QPipCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl; 00117 #endif 00118 if(pMom<=lastTH) 00119 { 00120 #ifdef debug 00121 G4cout<<"G4QPipCS::GCS:Found,P="<<pMom<<" < Threshold="<<lastTH<<",CS=0"<<G4endl; 00122 #endif 00123 return 0.; // Energy is below the Threshold value 00124 } 00125 lastP =colP [i]; // Last Momentum (A-dependent) 00126 lastCS =colCS[i]; // Last CrossSect (A-dependent) 00127 if(std::fabs(lastP-pMom)<tolerance*pMom) 00128 //if(lastP==pMom) // VI do not use tolerance 00129 { 00130 #ifdef debug 00131 G4cout<<"..G4QPipCS::GetCS:DoNothing,P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; 00132 #endif 00133 //CalculateCrossSection(fCS,-1,j,211,lastZ,lastN,pMom); // Update param's only 00134 return lastCS*millibarn; // Use theLastCS 00135 } 00136 in = true; // This is the case when the isotop is found in DB 00137 // Momentum pMom is in IU ! @@ Units 00138 #ifdef debug 00139 G4cout<<"G4QPipCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl; 00140 #endif 00141 lastCS=CalculateCrossSection(fCS,-1,j,211,lastZ,lastN,pMom); // read & update 00142 #ifdef debug 00143 G4cout<<"G4QPipNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl; 00144 #endif 00145 if(lastCS<=0. && pMom>lastTH) // Correct the threshold (@@ No intermediate Zeros) 00146 { 00147 #ifdef debug 00148 G4cout<<"G4QPipNCS::GetCS: New P="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00149 #endif 00150 lastCS=0.; 00151 lastTH=pMom; 00152 } 00153 break; // Go out of the LOOP 00154 } 00155 #ifdef debug 00156 G4cout<<"-->G4QPipNucCrossSec::GetCrosSec: pPDG=211, j="<<j<<", N="<<colN[i] 00157 <<",Z["<<i<<"]="<<colZ[i]<<G4endl; 00158 #endif 00159 j++; // Increment a#0f records found in DB 00160 } 00161 #ifdef debug 00162 G4cout<<"-?-G4QPipCS::GetCS:R,Z="<<tgZ<<",N="<<tgN<<",in="<<in<<",j="<<j<<" ?"<<G4endl; 00163 #endif 00164 if(!in) // This isotope has not been calculated previously 00165 { 00166 #ifdef debug 00167 G4cout<<"^^^G4QPipCS::GeCS:CalcNew P="<<pMom<<", f="<<fCS<<", lastI="<<lastI<<G4endl; 00168 #endif 00170 lastCS=CalculateCrossSection(fCS,0,j,211,lastZ,lastN,pMom); //calculate & create 00171 //if(lastCS>0.) // It means that the AMBD was initialized 00172 //{ 00173 00174 lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last 00175 #ifdef debug 00176 G4cout<<"G4QPipCrossSection::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl; 00177 #endif 00178 colN.push_back(tgN); 00179 colZ.push_back(tgZ); 00180 colP.push_back(pMom); 00181 colTH.push_back(lastTH); 00182 colCS.push_back(lastCS); 00183 #ifdef debug 00184 G4cout<<"G4QPipNCS::GetCrosSec:lCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl; 00185 #endif 00186 //} // M.K. Presence of H1 with high threshold breaks the syncronization 00187 #ifdef pdebug 00188 G4cout<<"G4QPipCS::GeCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00189 #endif 00190 return lastCS*millibarn; 00191 } // End of creation of the new set of parameters 00192 else 00193 { 00194 #ifdef debug 00195 G4cout<<"G4QPipNucCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl; 00196 #endif 00197 colP[lastI]=pMom; 00198 colCS[lastI]=lastCS; 00199 } 00200 } // End of parameters udate 00201 else if(pMom<=lastTH) 00202 { 00203 #ifdef debug 00204 G4cout<<"G4QPipCS::GetCS: Current P="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl; 00205 #endif 00206 return 0.; // Momentum is below the Threshold Value -> CS=0 00207 } 00208 else if(std::fabs(lastP-pMom)<tolerance*pMom) 00209 //else if(lastP==pMom) // VI do not use tolerance 00210 { 00211 #ifdef debug 00212 G4cout<<"..G4QPCS::GetCS:OldNZ&P="<<lastP<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; 00213 #endif 00214 return lastCS*millibarn; // Use theLastCS 00215 } 00216 else // It is the last used -> use the current tables 00217 { 00218 #ifdef debug 00219 G4cout<<"-!-G4QPCS::GetCS:UseCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl; 00220 #endif 00221 lastCS=CalculateCrossSection(fCS,1,j,211,lastZ,lastN,pMom); // Only read and UpdateDB 00222 lastP=pMom; 00223 } 00224 #ifdef debug 00225 G4cout<<"==>G4QPipCS::GetCroSec:P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00226 #endif 00227 return lastCS*millibarn; 00228 }
G4VQCrossSection * G4QPionPlusNuclearCrossSection::GetPointer | ( | ) | [static] |
Definition at line 63 of file G4QPionPlusNuclearCrossSection.cc.
Referenced by G4QHadronInelasticDataSet::GetIsoCrossSection(), and G4QInelastic::GetMeanFreePath().
00064 { 00065 static G4QPionPlusNuclearCrossSection theCrossSection; //**Static body of Cross Section** 00066 return &theCrossSection; 00067 }