#include <G4QPhotonNuclearCrossSection.hh>
Inheritance diagram for G4QPhotonNuclearCrossSection:
Public Member Functions | |
| ~G4QPhotonNuclearCrossSection () | |
| virtual G4double | GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=0) |
| G4double | CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum) |
| G4double | ThresholdEnergy (G4int Z, G4int N, G4int PDG=22) |
Static Public Member Functions | |
| static G4VQCrossSection * | GetPointer () |
Protected Member Functions | |
| G4QPhotonNuclearCrossSection () | |
Definition at line 51 of file G4QPhotonNuclearCrossSection.hh.
| G4QPhotonNuclearCrossSection::G4QPhotonNuclearCrossSection | ( | ) | [inline, protected] |
| G4QPhotonNuclearCrossSection::~G4QPhotonNuclearCrossSection | ( | ) |
| G4double G4QPhotonNuclearCrossSection::CalculateCrossSection | ( | G4bool | CS, | |
| G4int | F, | |||
| G4int | I, | |||
| G4int | PDG, | |||
| G4int | Z, | |||
| G4int | N, | |||
| G4double | Momentum | |||
| ) | [virtual] |
Implements G4VQCrossSection.
Definition at line 308 of file G4QPhotonNuclearCrossSection.cc.
References G4VQCrossSection::EquLinearFit(), G4cerr, G4cout, and G4endl.
Referenced by GetCrossSection().
00310 { 00311 #ifdef debug 00312 G4cout<<"G4QPhotonNucCrossSection::CalculateCrossSection: ***Called***"<<G4endl; 00313 #endif 00314 static const G4double THmin=2.; // minimum Energy Threshold 00315 static const G4double dE=1.; // step for the GDR table 00316 static const G4int nL=105; // A#of GDResonance points in E (1 MeV from 2 to 106) 00317 static const G4double Emin=THmin+(nL-1)*dE; // minE for the HighE part 00318 static const G4double Emax=50000.; // maxE for the HighE part 00319 static const G4int nH=224; // A#of HResonance points in lnE 00320 static const G4double milE=std::log(Emin); // Low logarithm energy for the HighE part 00321 static const G4double malE=std::log(Emax); // High logarithm energy (each 2.75 %) 00322 static const G4double dlE=(malE-milE)/(nH-1); // Step in log energy in the HighE part 00323 // 00324 //static const G4double shd=1.075-.0023*log(2.); // HE PomShadowing(D) 00325 static const G4double shd=1.0734; // HE PomShadowing(D) 00326 static const G4double shc=0.072; // HE Shadowing constant 00327 static const G4double poc=0.0375; // HE Pomeron coefficient 00328 static const G4double pos=16.5; // HE Pomeron shift 00329 static const G4double reg=.11; // HE Reggeon slope 00330 static const G4double shp=1.075; // HE PomShadowing(P) 00331 // 00332 // Associative memory for acceleration 00333 static std::vector <G4double> spA; // shadowing coefficients (A-dependent) 00334 // 00335 onlyCS=CS; // Flag to calculate only CS (not Si/Bi) 00336 #ifdef pdebug 00337 G4cout<<"G4QPhotonNucCS::CalcCS: P="<<Energy<<", F="<<F<<", I="<<I<<", Z="<<targZ 00338 <<", N="<<targN<<", onlyCS="<<CS<<",E="<<Energy<<",th="<<THmin<<G4endl; 00339 if(F==-27) return 0.; 00340 #endif 00341 //if (Energy<THmin) 00342 //{ 00343 // lastE=0.; 00344 // lastSig=0.; 00345 #ifdef debug 00346 // G4cout<<"---> G4QMuonNucCS::CalcCS: CS=0 as E="<<Energy<<" < "<<THmin<<G4endl; 00347 #endif 00348 // return 0.; // @@ This can be dangerouse for the heaviest nuc.! 00349 //} 00350 G4double sigma=0.; 00351 G4double A=targN+targZ; 00352 if(F<=0) // This isotope was not the last used isotop 00353 { 00354 if(F<0) // This isotope was found in DAMDB =-------=> RETRIEVE 00355 { 00356 lastGDR=(*GDR)[I]; // Pointer to prepared GDR cross sections 00357 lastHEN=(*HEN)[I]; // Pointer to prepared High Energy cross sections 00358 lastSP =spA[I]; // Shadowing coefficient for UHE 00359 } 00360 else // This isotope wasn't calculated previously => CREATE 00361 { 00362 G4double lnA=std::log(A); // The nucleus is not found in DB. It is new. 00363 if(A==1.) lastSP=1.; // The Reggeon shadowing (A=1) 00364 else lastSP=A*(1.-shc*lnA); // The Reggeon shadowing 00365 #ifdef debug 00366 G4cout<<">>>G4QPhotonNuclearCrossSect::CalcCS:lnA="<<lnA<<",lastSP="<<lastSP<<G4endl; 00367 #endif 00368 #ifdef debug3 00369 if(A==3) G4cout<<"G4QPhotonNuclearCrossSection::CalcCS: lastSP="<<lastSP<<G4endl; 00370 #endif 00371 lastGDR = new G4double[nL]; // Allocate memory for the new GDR cross sections 00372 lastHEN = new G4double[nH]; // Allocate memory for the new HEN cross sections 00373 G4int er=GetFunctions(A,lastGDR,lastHEN);// set newZeroPosition and fill theFunctions 00374 if(er<1) G4cerr<<"***G4QPhotNucCrosSec::CalcCrossSection: A="<<A<<" failed"<<G4endl; 00375 #ifdef pdebug 00376 G4cout<<">>G4QPhotNucCS::CalcCS:**GDR/HEN're made for A="<<A<<"** GFEr="<<er<<G4endl; 00377 #endif 00378 // *** The synchronization check *** 00379 G4int sync=GDR->size(); 00380 if(sync!=I) G4cerr<<"***G4QPhotoNuclearCS::CalcCS: PDG=22, S="<<sync<<"#"<<I<<G4endl; 00381 GDR->push_back(lastGDR); // added GDR, found by AH 10/7/02 00382 HEN->push_back(lastHEN); // added HEN, found by AH 10/7/02 00383 spA.push_back(lastSP); // Pomeron Shadowing 00384 } // End of creation of the new set of parameters 00385 } // End of parameters udate 00386 // =-------------------= NOW the Magic Formula =--------------------------= 00387 if (Energy<lastTH) return 0.; // It must be already checked in the interface 00388 else if (Energy<=Emin) // GDR region (approximated in E, not in lnE) 00389 { 00390 #ifdef debug 00391 G4cout<<"G4QPhNCS::CalcCS:bGDR A="<<A<<", nL="<<nL<<",TH="<<THmin<<",dE="<<dE<<G4endl; 00392 #endif 00393 if(A<=1. || dE <= THmin) sigma=0.; // No GDR for A=1 00394 else sigma=EquLinearFit(Energy,nL,THmin,dE,lastGDR); 00395 #ifdef debugn 00396 if(sigma<0.) 00397 G4cout<<"G4QPhoNucCS::CalcCS:A="<<A<<",E="<<Energy<<",T="<<THmin<<",dE="<<dE<<G4endl; 00398 #endif 00399 } 00400 else if (Energy<Emax) // High Energy region 00401 { 00402 G4double lE=std::log(Energy); 00403 #ifdef pdebug 00404 G4cout<<"G4QPhotNucCS::CalcCS:lE="<<milE<<",n="<<nH<<",iE="<<milE<<",dE="<<dlE<<",HEN=" 00405 <<lastHEN<<",A="<<A<<G4endl; 00406 #endif 00407 if(lE > milE) sigma=EquLinearFit(lE,nH,milE,dlE,lastHEN); 00408 else sigma=0.; 00409 } 00410 else // UHE region (calculation, not frequent) 00411 { 00412 G4double lE=std::log(Energy); 00413 G4double sh=shd; 00414 if(A==1.)sh=shp; 00415 sigma=lastSP*(poc*(lE-pos)+sh*std::exp(-reg*lE)); 00416 } 00417 #ifdef debug 00418 G4cout<<"G4QPhotonNuclearCrossSection::CalcCS: sigma="<<sigma<<G4endl; 00419 #endif 00420 #ifdef debug 00421 if(Energy>45000.&&Energy<60000.) 00422 G4cout<<"G4QPhotoNucCS::GetCS: A="<<A<<", E="<<Energy<<",CS="<<sigma<<G4endl; 00423 #endif 00424 if(sigma<0.) return 0.; 00425 return sigma; 00426 }
| G4double G4QPhotonNuclearCrossSection::GetCrossSection | ( | G4bool | fCS, | |
| G4double | pMom, | |||
| G4int | tgZ, | |||
| G4int | tgN, | |||
| G4int | pPDG = 0 | |||
| ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 88 of file G4QPhotonNuclearCrossSection.cc.
References CalculateCrossSection(), G4cout, G4endl, and ThresholdEnergy().
00090 { 00091 static G4int j; // A#0f records found in DB for this projectile 00092 static std::vector <G4int> colPDG;// Vector of the projectile PDG code 00093 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) 00094 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) 00095 static std::vector <G4double> colP; // Vector of last momenta for the reaction 00096 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction 00097 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction 00098 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** 00099 G4double pEn=pMom; 00100 #ifdef debug 00101 G4cout<<"G4QPhCS::GetCS:>>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN 00102 <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz=" 00103 <<colN.size()<<G4endl; 00104 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00105 #endif 00106 if(!pPDG) 00107 { 00108 #ifdef debug 00109 G4cout<<"G4QPhCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl; 00110 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00111 #endif 00112 return 0.; // projectile PDG=0 is a mistake (?!) @@ 00113 } 00114 if(pPDG==22 && tgZ==1 && !tgN && pMom<150.*MeV) return 0.; // Examle of pre-threshold (A) 00115 G4bool in=false; // By default the isotope must be found in the AMDB 00116 if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope 00117 { 00118 in = false; // By default the isotope haven't be found in AMDB 00119 lastP = 0.; // New momentum history (nothing to compare with) 00120 lastPDG = pPDG; // The last PDG of the projectile 00121 lastN = tgN; // The last N of the calculated nucleus 00122 lastZ = tgZ; // The last Z of the calculated nucleus 00123 lastI = colN.size(); // Size of the Associative Memory DB in the heap 00124 j = 0; // A#0f records found in DB for this projectile 00125 if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found 00126 { // The nucleus with projPDG is found in AMDB 00127 if(colN[i]==tgN && colZ[i]==tgZ) 00128 { 00129 lastI=i; 00130 lastTH =colTH[i]; // Last THreshold (A-dependent) 00131 #ifdef debug 00132 G4cout<<"G4QPhCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl; 00133 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00134 #endif 00135 if(pEn<=lastTH) 00136 { 00137 #ifdef debug 00138 G4cout<<"G4QPhCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=0"<<G4endl; 00139 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00140 #endif 00141 return 0.; // Energy is below the Threshold value 00142 } 00143 lastP =colP [i]; // Last Momentum (A-dependent) 00144 lastCS =colCS[i]; // Last CrossSect (A-dependent) 00145 if(lastP == pMom) 00146 { 00147 #ifdef debug 00148 G4cout<<"G4QPhCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; 00149 #endif 00150 CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only 00151 return lastCS*millibarn; // Use theLastCS 00152 } 00153 in = true; // This is the case when the isotop is found in DB 00154 // Momentum pMom is in IU ! @@ Units 00155 #ifdef pdebug 00156 G4cout<<"G4QPhCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl; 00157 #endif 00158 lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update 00159 #ifdef debug 00160 G4cout<<"G4QPhCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl; 00161 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00162 #endif 00163 if(lastCS<=0. && pEn>lastTH) // Correct the threshold 00164 { 00165 #ifdef debug 00166 G4cout<<"G4QPhCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00167 #endif 00168 lastTH=pEn; 00169 } 00170 break; // Go out of the LOOP 00171 } 00172 #ifdef debug 00173 G4cout<<"---G4QPhCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i] 00174 <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl; 00175 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00176 #endif 00177 j++; // Increment a#0f records found in DB for this pPDG 00178 } 00179 if(!in) // This nucleus has not been calculated previously 00180 { 00181 #ifdef debug 00182 G4cout<<"G4QPhCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl; 00183 #endif 00185 lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create 00186 if(lastCS<=0.) 00187 { 00188 lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last 00189 #ifdef debug 00190 G4cout<<"G4QPhCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl; 00191 #endif 00192 if(pEn>lastTH) 00193 { 00194 #ifdef debug 00195 G4cout<<"G4QPhCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; 00196 #endif 00197 lastTH=pEn; 00198 } 00199 } 00200 #ifdef debug 00201 G4cout<<"G4QPhCS::GetCrosSec: New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl; 00202 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00203 #endif 00204 colN.push_back(tgN); 00205 colZ.push_back(tgZ); 00206 colPDG.push_back(pPDG); 00207 colP.push_back(pMom); 00208 colTH.push_back(lastTH); 00209 colCS.push_back(lastCS); 00210 #ifdef debug 00211 G4cout<<"G4QPhCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00212 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00213 #endif 00214 return lastCS*millibarn; 00215 } // End of creation of the new set of parameters 00216 else 00217 { 00218 #ifdef debug 00219 G4cout<<"G4QPrCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl; 00220 #endif 00221 colP[lastI]=pMom; 00222 colPDG[lastI]=pPDG; 00223 colCS[lastI]=lastCS; 00224 } 00225 } // End of parameters udate 00226 else if(pEn<=lastTH) 00227 { 00228 #ifdef debug 00229 G4cout<<"G4QPhCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl; 00230 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00231 #endif 00232 return 0.; // Momentum is below the Threshold Value -> CS=0 00233 } 00234 else if(lastP == pMom) 00235 { 00236 #ifdef debug 00237 G4cout<<"G4QPhCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", CS="<<lastCS*millibarn<<G4endl; 00238 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00239 #endif 00240 return lastCS*millibarn; // Use theLastCS 00241 } 00242 else 00243 { 00244 #ifdef debug 00245 G4cout<<"G4QPhCS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl; 00246 #endif 00247 lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB 00248 lastP=pMom; 00249 } 00250 #ifdef debug 00251 G4cout<<"G4QPhCS::GetCroSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; 00252 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST 00253 #endif 00254 return lastCS*millibarn; 00255 }
| G4VQCrossSection * G4QPhotonNuclearCrossSection::GetPointer | ( | ) | [static] |
Definition at line 70 of file G4QPhotonNuclearCrossSection.cc.
Referenced by G4QInelastic::GetMeanFreePath(), G4QInelastic::PostStepDoIt(), and G4QAtomicElectronScattering::PostStepDoIt().
00071 { 00072 static G4QPhotonNuclearCrossSection theCrossSection; //**Static body of Cross Section** 00073 return &theCrossSection; 00074 }
| G4double G4QPhotonNuclearCrossSection::ThresholdEnergy | ( | G4int | Z, | |
| G4int | N, | |||
| G4int | PDG = 22 | |||
| ) | [virtual] |
Reimplemented from G4VQCrossSection.
Definition at line 263 of file G4QPhotonNuclearCrossSection.cc.
References G4cout, G4endl, G4NucleiProperties::GetNuclearMass(), and G4NucleiProperties::IsInStableTable().
Referenced by GetCrossSection().
00264 { 00265 // CHIPS - Direct GEANT 00266 //static const G4double mNeut = G4QPDGCode(2112).GetMass(); 00267 //static const G4double mProt = G4QPDGCode(2212).GetMass(); 00268 static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/MeV; 00269 static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/MeV; 00270 static const G4double mAlph = G4NucleiProperties::GetNuclearMass(4,2)/MeV; 00271 // --------- 00272 static const G4double infEn = 9.e27; 00273 00274 G4int A=Z+N; 00275 if(A<1) return infEn; 00276 else if(A==1) return 134.9766; // Pi0 threshold for the nucleon 00277 // CHIPS - Direct GEANT 00278 //G4double mT= G4QPDGCode(111).GetNuclMass(Z,N,0); 00279 G4double mT= 0.; 00280 if(Z&&G4NucleiProperties::IsInStableTable(A,Z)) 00281 mT = G4NucleiProperties::GetNuclearMass(A,Z)/MeV; 00282 else return 0.; // If it is not in the Table of Stable Nuclei, then the Threshold=0 00283 G4double mP= infEn; 00284 if(A>1&&Z&&G4NucleiProperties::IsInStableTable(A-1,Z-1)) 00285 mP = G4NucleiProperties::GetNuclearMass(A-1,Z-1)/MeV;// ResNucMass for a proton 00286 00287 G4double mN= infEn; 00288 if(A>1&&G4NucleiProperties::IsInStableTable(A-1,Z)) 00289 mN = G4NucleiProperties::GetNuclearMass(A-1,Z)/MeV; // ResNucMass for a neutron 00290 00291 G4double mA= infEn; 00292 if(A>4&&Z>1&&G4NucleiProperties::IsInStableTable(A-4,Z-2)) 00293 mA=G4NucleiProperties::GetNuclearMass(A-4,Z-2)/MeV; // ResNucMass for an alpha 00294 00295 G4double dP= mP +mProt - mT; 00296 G4double dN= mN +mNeut - mT; 00297 G4double dA= mA +mAlph - mT; 00298 #ifdef debug 00299 G4cout<<"G4QPhotoNucCS::ThreshEn: mP="<<mP<<",dP="<<dP<<",mN="<<mN<<",dN="<<dN<<",mA=" 00300 <<mA<<",dA="<<dA<<",mT="<<mT<<",A="<<A<<",Z="<<Z<<G4endl; 00301 #endif 00302 if(dP<dN)dN=dP; 00303 if(dA<dN)dN=dA; 00304 return dN; 00305 }
1.4.7