#include <G4QNeutronCaptureRatio.hh>
Public Member Functions | |
~G4QNeutronCaptureRatio () | |
G4double | GetRatio (G4double pIU, G4int tgZ, G4int tgN) |
Static Public Member Functions | |
static G4QNeutronCaptureRatio * | GetPointer () |
Protected Member Functions | |
G4QNeutronCaptureRatio () |
Definition at line 53 of file G4QNeutronCaptureRatio.hh.
G4QNeutronCaptureRatio::G4QNeutronCaptureRatio | ( | ) | [inline, protected] |
G4QNeutronCaptureRatio::~G4QNeutronCaptureRatio | ( | ) | [inline] |
G4QNeutronCaptureRatio * G4QNeutronCaptureRatio::GetPointer | ( | ) | [static] |
Definition at line 48 of file G4QNeutronCaptureRatio.cc.
Referenced by G4QNGamma::GetMeanFreePath(), and G4QInelastic::GetMeanFreePath().
00049 { 00050 static G4QNeutronCaptureRatio theRatios; // *** Static body of the NeutCaptureRatio *** 00051 return &theRatios; 00052 }
Definition at line 55 of file G4QNeutronCaptureRatio.cc.
References G4cout, G4endl, and CLHEP::detail::n.
Referenced by G4QNGamma::GetMeanFreePath(), and G4QInelastic::GetMeanFreePath().
00056 { 00057 // Table parameters 00058 static const G4int npp=100; // Number of steps in the R(s) LinTable 00059 static const G4int mpp=npp+1; // Number of elements in the R(s) LinTable 00060 static const G4double pma=6.; // The first LinTabEl(mom=0)=1., mom>pma -> logTab 00061 static const G4double dp=pma/npp; // Step of the linear Table 00062 static const G4int nls=150; // Number of steps in the R(lns) logTable 00063 static const G4int mls=nls+1; // Number of elements in the R(lns) logTable 00064 static const G4double lpi=1.79; // The min ln(p) logTabEl(p=5.99 < pma=6.) 00065 static const G4double lpa=8.; // The max ln(p) logTabEl(p=5.99 - 2981 GeV) 00066 static const G4double mi=std::exp(lpi);// The min mom of logTabEl(~ 5.99 GeV) 00067 static const G4double max_s=std::exp(lpa);// The max mom of logTabEl(~ 2981 GeV) 00068 static const G4double dl=(lpa-lpi)/nls;// Step of the logarithmic Table 00069 static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table 00070 static const G4double toler=.0001; // Tolarence (GeV) defining the same momentum 00071 static G4double lastP=0.; // Last mometum value for which R was calculated 00072 static G4double lastR=0.; // Last ratio R which was calculated 00073 // Local Associative Data Base: 00074 static std::vector<G4int> vZ; // Vector of calculated Z (target) 00075 static std::vector<G4int> vN; // Vector of calculated N (target) 00076 static std::vector<G4double> vH; // Vector of max mom initialized in the LinTable 00077 static std::vector<G4int> vJ; // Vector of topBin number initialized in LinTable 00078 static std::vector<G4double> vM; // Vector of relMax ln(p) initialized in LogTable 00079 static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable 00080 static std::vector<G4double*> vT; // Vector of pointers to LinTable in C++ heap 00081 static std::vector<G4double*> vL; // Vector of pointers to LogTable in C++ heap 00082 // Last values of the Associative Data Base: 00083 static G4int lastZ=0; // theLast of calculated A 00084 static G4int lastN=0; // theLast of calculated A 00085 static G4double lastH=0.; // theLast of max mom initialized in the LinTable 00086 static G4int lastJ=0; // theLast of topBinNumber initialized in LinTable 00087 static G4double lastM=0.; // theLast of relMax ln(p) initialized in LogTab. 00088 static G4int lastK=0; // theLast of topBinNumber initialized in LogTable 00089 static G4double* lastT=0; // theLast of pointer to LinTable in the C++ heap 00090 static G4double* lastL=0; // theLast of pointer to LogTable in the C++ heap 00091 // LogTable is created only if necessary. R(p>2981GeV) calcul by formula for any nuclei 00092 G4int A=tgN+tgZ; 00093 if(pIU > 50) return 0.; 00094 if(pIU > 30 && ((tgN==1 && tgZ==1) || (tgN==8 && tgZ==7))) return 0.; 00095 if(pIU > 20 && tgN==2 && tgZ==1) return 0.; 00096 if(pIU > 15 && ((tgN==1 && tgZ==2) || (tgN==8 && tgZ==8))) return 0.; 00097 if(pIU<toler || A<1) return 1.; // Fake use of toler as non zero number 00098 if(A>247) 00099 { 00100 G4cout<<"-*-Warning-*-G4NeutronCaptureRatio::GetRatio:A="<<A<<">247, return 0"<<G4endl; 00101 return 0.; 00102 } 00103 G4int nDB=vZ.size(); // A number of nuclei already initialized in AMDB 00104 if(nDB && lastZ==tgZ && lastN==tgN && std::fabs(pIU-lastP)<toler) return lastR; 00105 if(pIU>max_s) 00106 { 00107 lastR=CalcCap2In_Ratio(s,tgZ,tgN); //@@ Probably user ought to be notified about bigP 00108 return lastR; 00109 } 00110 G4bool found=false; 00111 G4int i=-1; 00112 if(nDB) for (i=0; i<nDB; i++) if(tgZ==vZ[i] && tgN==vN[i]) // Sirch for this Z,N in AMDB 00113 { 00114 found=true; // The (Z,N) is found 00115 break; 00116 } 00117 if(!nDB || !found) // Create new line in the AMDB 00118 { 00119 lastZ = tgZ; 00120 lastN = tgN; 00121 lastT = new G4double[mpp]; // Create the linear Table 00122 lastJ = static_cast<int>(pIU/dp)+1; // MaxBin to be initialized 00123 if(lastJ>npp) 00124 { 00125 lastJ=npp; 00126 lastH=pma; 00127 } 00128 else lastH = lastJ*dp; // Calculate max initialized s for LinTab 00129 G4double pv=0; 00130 lastT[0]=1.; 00131 for(G4int j=1; j<=lastJ; j++) // Calculate LinTab values 00132 { 00133 pv+=dp; 00134 lastT[j]=CalcCap2In_Ratio(pv,tgZ,tgN); // ?? 00135 } 00136 lastL=new G4double[mls]; // Create the logarithmic Table 00137 G4double ls=std::log(s); 00138 lastK = static_cast<int>((ls-lpi)/dl)+1; // MaxBin to be initialized in LogTaB 00139 if(lastK>nls) 00140 { 00141 lastK=nls; 00142 lastM=lpa-lpi; 00143 } 00144 else lastM = lastK*dl; // Calculate max initialized ln(s)-lpi for LogTab 00145 pv=mi; 00146 for(G4int j=0; j<=lastK; j++) // Calculate LogTab values 00147 { 00148 lastL[j]=CalcCap2In_Ratio(pv,tgZ,tgN); 00149 if(j!=lastK) pv*=edl; 00150 } 00151 i++; // Make a new record to AMDB and position on it 00152 vZ.push_back(lastZ); 00153 vN.push_back(lastN); 00154 vH.push_back(lastH); 00155 vJ.push_back(lastJ); 00156 vM.push_back(lastM); 00157 vK.push_back(lastK); 00158 vT.push_back(lastT); 00159 vL.push_back(lastL); 00160 } 00161 else // The A value was found in AMDB 00162 { 00163 lastZ=vZ[i]; 00164 lastN=vN[i]; 00165 lastH=vH[i]; 00166 lastJ=vJ[i]; 00167 lastM=vM[i]; 00168 lastK=vK[i]; 00169 lastT=vT[i]; 00170 lastL=vL[i]; 00171 if(s>lastH) // At least LinTab must be updated 00172 { 00173 G4int nextN=lastJ+1; // The next bin to be initialized 00174 if(lastJ<npp) 00175 { 00176 lastJ = static_cast<int>(pIU/dp)+1;// MaxBin to be initialized 00177 G4double pv=lastH; 00178 if(lastJ>npp) 00179 { 00180 lastJ=npp; 00181 lastH=pma; 00182 } 00183 else lastH = lastJ*dp; // Calculate max initialized s for LinTab 00184 for(G4int j=nextN; j<=lastJ; j++)// Calculate LogTab values 00185 { 00186 pv+=dp; 00187 lastT[j]=CalcCap2In_Ratio(pv,tgZ,tgN); 00188 } 00189 } // End of LinTab update 00190 if(lastJ>=nextN) 00191 { 00192 vH[i]=lastH; 00193 vJ[i]=lastJ; 00194 } 00195 G4int nextK=lastK+1; 00196 if(pIU>pma && lastK<nls) // LogTab must be updated 00197 { 00198 G4double pv=std::exp(lastM+lpi); // Define starting poit (lastM will be changed) 00199 G4double ls=std::log(s); 00200 lastK = static_cast<int>((ls-lpi)/dl)+1; // MaxBin to be initialized in LogTaB 00201 if(lastK>nls) 00202 { 00203 lastK=nls; 00204 lastM=lpa-lpi; 00205 } 00206 else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab 00207 for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values 00208 { 00209 pv*=edl; 00210 lastL[j]=CalcCap2In_Ratio(pv,tgZ,tgN); 00211 } 00212 } // End of LogTab update 00213 if(lastK>=nextK) 00214 { 00215 vM[i]=lastM; 00216 vK[i]=lastK; 00217 } 00218 } 00219 } 00220 // Now one can use tabeles to calculate the value 00221 if(pIU<pma) // Use linear table 00222 { 00223 G4int n=static_cast<int>(pIU/dp); // Low edge number of the bin 00224 G4double d=s-n*dp; // Linear shift 00225 G4double v=lastT[n]; // Base 00226 lastR=v+d*(lastT[n+1]-v)/dp; // Result 00227 } 00228 else // Use log table 00229 { 00230 G4double ls=std::log(pIU)-lpi; // ln(p)-l_min 00231 G4int n=static_cast<int>(ls/dl); // Low edge number of the bin 00232 G4double d=ls-n*dl; // Log shift 00233 G4double v=lastL[n]; // Base 00234 lastR=v+d*(lastL[n+1]-v)/dl; // Result 00235 } 00236 if(lastR<0.) lastR=0.; 00237 if(lastR>1.) lastR=1.; 00238 return lastR; 00239 } // End of GetRatio