1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // The lust update: M.V. Kossov, CERN/ITEP(Moscow) 19-Aug-07 |
---|
28 | // GEANT4 tag $Name: geant4-09-03-cand-01 $ |
---|
29 | // |
---|
30 | // |
---|
31 | // G4 Physics class: G4QIonIonCrossSection for gamma+A cross sections |
---|
32 | // Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-03 |
---|
33 | // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Feb-04 |
---|
34 | // -------------------------------------------------------------------------------- |
---|
35 | // **************************************************************************************** |
---|
36 | // ***** This HEADER is a property of the CHIPS hadronic package in Geant4 (M. Kosov) ***** |
---|
37 | // *********** DO NOT MAKE ANY CHANGE without approval of Mikhail.Kossov@cern.ch ********** |
---|
38 | // **************************************************************************************** |
---|
39 | // Short description: CHIPS cross-sectons for Ion-Ion interactions |
---|
40 | // --------------------------------------------------------------- |
---|
41 | // |
---|
42 | //#define debug |
---|
43 | //#define pdebug |
---|
44 | //#define debug3 |
---|
45 | //#define debugn |
---|
46 | //#define debugs |
---|
47 | |
---|
48 | #include "G4QIonIonCrossSection.hh" |
---|
49 | |
---|
50 | // Initialization of the |
---|
51 | G4double* G4QIonIonCrossSection::lastLENI=0;// Pointer to the lastArray of LowEn Inelast CS |
---|
52 | G4double* G4QIonIonCrossSection::lastHENI=0;// Pointer to the lastArray of HighEn InelastCS |
---|
53 | G4double* G4QIonIonCrossSection::lastLENE=0;// Pointer to the lastArray of LowEn Elastic CS |
---|
54 | G4double* G4QIonIonCrossSection::lastHENE=0;// Pointer to the lastArray of HighEn ElasticCS |
---|
55 | G4int G4QIonIonCrossSection::lastPDG=0; // The last PDG code of the projectile |
---|
56 | G4int G4QIonIonCrossSection::lastN=0; // The last N of calculated nucleus |
---|
57 | G4int G4QIonIonCrossSection::lastZ=0; // The last Z of calculated nucleus |
---|
58 | G4double G4QIonIonCrossSection::lastP=0.; // Last used in cross section Momentum |
---|
59 | G4double G4QIonIonCrossSection::lastTH=0.; // Last threshold momentum |
---|
60 | G4double G4QIonIonCrossSection::lastICS=0.;// Last value of the Inelastic Cross Section |
---|
61 | G4double G4QIonIonCrossSection::lastECS=0.;// Last value of the Elastic Cross Section |
---|
62 | G4int G4QIonIonCrossSection::lastI=0; // The last position in the DAMDB |
---|
63 | |
---|
64 | // Returns Pointer to the G4VQCrossSection class |
---|
65 | G4VQCrossSection* G4QIonIonCrossSection::GetPointer() |
---|
66 | { |
---|
67 | static G4QIonIonCrossSection theCrossSection; //**Static body of Cross Section** |
---|
68 | return &theCrossSection; |
---|
69 | } |
---|
70 | |
---|
71 | // The main member function giving the collision cross section (P is in IU, CS is in mb) |
---|
72 | // Make pMom in independent units !(Now it is MeV): fCS=true->Inelastic, fCS=false->Elastic |
---|
73 | G4double G4QIonIonCrossSection::GetCrossSection(G4bool fCS, G4double pMom, G4int tZ, |
---|
74 | G4int tN, G4int pPDG) |
---|
75 | { |
---|
76 | static G4int j; // A#0f records found in DB for this projectile |
---|
77 | static std::vector <G4int> colPDG;// Vector of the projectile PDG code |
---|
78 | static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) |
---|
79 | static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) |
---|
80 | static std::vector <G4double> colP; // Vector of last momenta for the reaction |
---|
81 | static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction |
---|
82 | static std::vector <G4double> colICS;// Vector of last inelastic cross-sections |
---|
83 | static std::vector <G4double> colECS;// Vector of last elastic cross-sections |
---|
84 | // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** |
---|
85 | #ifdef pdebug |
---|
86 | G4cout<<"G4QIICS::GetCS:>>> f="<<fCS<<", Z="<<tZ<<"("<<lastZ<<"), N="<<tN<<"("<<lastN |
---|
87 | <<"),PDG="<<pPDG<<"("<<lastPDG<<"), p="<<pMom<<"("<<lastTH<<")"<<",Sz=" |
---|
88 | <<colN.size()<<G4endl; |
---|
89 | #endif |
---|
90 | if(!pPDG) |
---|
91 | { |
---|
92 | #ifdef pdebug |
---|
93 | G4cout<<"G4QIonIonCS::GetCS: *** Found pPDG="<<pPDG<<" ====> CS=0"<<G4endl; |
---|
94 | #endif |
---|
95 | return 0.; // projectile PDG=0 is a mistake (?!) @@ |
---|
96 | } |
---|
97 | G4bool in=false; // By default the isotope must be found in the AMDB |
---|
98 | if(tN!=lastN || tZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope |
---|
99 | { |
---|
100 | in = false; // By default the isotope haven't be found in AMDB |
---|
101 | lastP = 0.; // New momentum history (nothing to compare with) |
---|
102 | lastPDG = pPDG; // The last PDG of the projectile |
---|
103 | lastN = tN; // The last N of the calculated nucleus |
---|
104 | lastZ = tZ; // The last Z of the calculated nucleus |
---|
105 | lastI = colN.size(); // Size of the Associative Memory DB in the heap |
---|
106 | j = 0; // A#0f records found in DB for this projectile |
---|
107 | #ifdef pdebug |
---|
108 | G4cout<<"G4QIICS::GetCS:FindI="<<lastI<<",pPDG="<<pPDG<<",tN="<<tN<<",tZ="<<tZ<<G4endl; |
---|
109 | #endif |
---|
110 | if(lastI) for(G4int i=0; i<lastI; i++) // Loop over all DB |
---|
111 | { // The nucleus with projPDG is found in AMDB |
---|
112 | #ifdef pdebug |
---|
113 | G4cout<<"G4QII::GCS:P="<<colPDG[i]<<",N="<<colN[i]<<",Z="<<colZ[i]<<",j="<<j<<G4endl; |
---|
114 | #endif |
---|
115 | if(colPDG[i]==pPDG && colN[i]==tN && colZ[i]==tZ) |
---|
116 | { |
---|
117 | lastI=i; |
---|
118 | lastTH =colTH[i]; // Last THreshold (A-dependent) |
---|
119 | #ifdef pdebug |
---|
120 | G4cout<<"G4QIICS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl; |
---|
121 | #endif |
---|
122 | if(pMom<=lastTH) |
---|
123 | { |
---|
124 | #ifdef pdebug |
---|
125 | G4cout<<"G4QIICS::GetCS:Found P="<<pMom<<"<Threshold="<<lastTH<<"->XS=0"<<G4endl; |
---|
126 | #endif |
---|
127 | return 0.; // Energy is below the Threshold value |
---|
128 | } |
---|
129 | lastP =colP [i]; // Last Momentum (A-dependent) |
---|
130 | lastICS=colICS[i]; // Last Inelastic Cross-Section (A-dependent) |
---|
131 | lastECS=colECS[i]; // Last Elastic Cross-Section (A-dependent) |
---|
132 | if(std::fabs(lastP/pMom-1.)<tolerance) |
---|
133 | { |
---|
134 | #ifdef pdebug |
---|
135 | G4cout<<"G4QIonIonCS::GetCS:P="<<pMom<<",InXS="<<lastICS*millibarn<<",ElXS=" |
---|
136 | <<lastECS*millibarn<<G4endl; |
---|
137 | #endif |
---|
138 | CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only |
---|
139 | if(fCS) return lastICS*millibarn; // Use theLastInelasticCS |
---|
140 | return lastECS*millibarn; // Use theLastElasticCS |
---|
141 | } |
---|
142 | in = true; // This is the case when the isotop is found in DB |
---|
143 | // Momentum pMom is in IU ! @@ Units |
---|
144 | #ifdef pdebug |
---|
145 | G4cout<<"G4QIICS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl; |
---|
146 | #endif |
---|
147 | lastICS=CalculateCrossSection( true,-1,j,lastPDG,lastZ,lastN,pMom);// read & update |
---|
148 | lastECS=CalculateCrossSection(false,-1,j,lastPDG,lastZ,lastN,pMom);// read & update |
---|
149 | #ifdef pdebug |
---|
150 | G4cout<<"G4QIonIonCS::GetCS:=>New(inDB) InCS="<<lastICS<<",ElCS="<<lastECS<<G4endl; |
---|
151 | #endif |
---|
152 | if((lastICS<=0. || lastECS<=0.) && pMom>lastTH) // Correct the threshold |
---|
153 | { |
---|
154 | #ifdef pdebug |
---|
155 | G4cout<<"G4QIonIonCS::GetCS:New,T="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl; |
---|
156 | #endif |
---|
157 | lastTH=pMom; |
---|
158 | } |
---|
159 | break; // Go out of the LOOP |
---|
160 | } |
---|
161 | #ifdef pdebug |
---|
162 | G4cout<<"--->G4QIonIonCrossSec::GetCrosSec: pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i] |
---|
163 | <<",Z["<<i<<"]="<<colZ[i]<<",PDG="<<colPDG[i]<<G4endl; |
---|
164 | #endif |
---|
165 | j++; // Increment a#0f records found in DB for this pPDG |
---|
166 | } |
---|
167 | if(!in) // This nucleus has not been calculated previously |
---|
168 | { |
---|
169 | #ifdef pdebug |
---|
170 | G4cout<<"G4QIICS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl; |
---|
171 | #endif |
---|
172 | //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU) |
---|
173 | lastICS=CalculateCrossSection(true ,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create |
---|
174 | lastECS=CalculateCrossSection(false,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create |
---|
175 | if(lastICS<=0. || lastECS<=0.) |
---|
176 | { |
---|
177 | lastTH = ThresholdEnergy(tZ, tN); // Threshold Energy=Mom=0 which is now the last |
---|
178 | #ifdef pdebug |
---|
179 | G4cout<<"G4QIonIonCrossSect::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl; |
---|
180 | #endif |
---|
181 | if(pMom>lastTH) |
---|
182 | { |
---|
183 | #ifdef pdebug |
---|
184 | G4cout<<"G4QIonIonCS::GetCS:1-st,P="<<pMom<<">Thresh="<<lastTH<<"->XS=0"<<G4endl; |
---|
185 | #endif |
---|
186 | lastTH=pMom; |
---|
187 | } |
---|
188 | } |
---|
189 | #ifdef pdebug |
---|
190 | G4cout<<"G4QIICS::GetCS: *New* ICS="<<lastICS<<", ECS="<<lastECS<<",N="<<lastN<<",Z=" |
---|
191 | <<lastZ<<G4endl; |
---|
192 | #endif |
---|
193 | colN.push_back(tN); |
---|
194 | colZ.push_back(tZ); |
---|
195 | colPDG.push_back(pPDG); |
---|
196 | colP.push_back(pMom); |
---|
197 | colTH.push_back(lastTH); |
---|
198 | colICS.push_back(lastICS); |
---|
199 | colECS.push_back(lastECS); |
---|
200 | #ifdef pdebug |
---|
201 | G4cout<<"G4QIICS::GetCS:*1st*, P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn |
---|
202 | <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl; |
---|
203 | #endif |
---|
204 | if(fCS) return lastICS*millibarn; // Use theLastInelasticCS |
---|
205 | return lastECS*millibarn; // Use theLastElasticCS |
---|
206 | } // End of creation of the new set of parameters |
---|
207 | else |
---|
208 | { |
---|
209 | #ifdef pdebug |
---|
210 | G4cout<<"G4QIICS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl; |
---|
211 | #endif |
---|
212 | colP[lastI]=pMom; |
---|
213 | colPDG[lastI]=pPDG; |
---|
214 | colICS[lastI]=lastICS; |
---|
215 | colECS[lastI]=lastECS; |
---|
216 | } |
---|
217 | } // End of parameters udate |
---|
218 | else if(pMom<=lastTH) |
---|
219 | { |
---|
220 | #ifdef pdebug |
---|
221 | G4cout<<"G4QIICS::GetCS: Current T="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl; |
---|
222 | #endif |
---|
223 | return 0.; // Momentum is below the Threshold Value -> CS=0 |
---|
224 | } |
---|
225 | else if(std::fabs(lastP/pMom-1.)<tolerance) |
---|
226 | { |
---|
227 | #ifdef pdebug |
---|
228 | G4cout<<"G4QIICS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", InCS="<<lastICS*millibarn |
---|
229 | <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl; |
---|
230 | #endif |
---|
231 | if(fCS) return lastICS*millibarn; // Use theLastInelasticCS |
---|
232 | return lastECS*millibarn; // Use theLastElasticCS |
---|
233 | } |
---|
234 | else |
---|
235 | { |
---|
236 | #ifdef pdebug |
---|
237 | G4cout<<"G4QIICS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl; |
---|
238 | #endif |
---|
239 | lastICS=CalculateCrossSection( true,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB |
---|
240 | lastECS=CalculateCrossSection(false,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB |
---|
241 | lastP=pMom; |
---|
242 | } |
---|
243 | #ifdef pdebug |
---|
244 | G4cout<<"G4QIICS::GetCroSec:*End*,P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn<<", ElCS=" |
---|
245 | <<lastECS*millibarn<<"(mb)"<<G4endl; |
---|
246 | #endif |
---|
247 | if(fCS) return lastICS*millibarn; // Use theLastInelasticCS |
---|
248 | return lastECS*millibarn; // Use theLastElasticCS |
---|
249 | } |
---|
250 | |
---|
251 | // The main member function giving the A-A cross section (Momentum in MeV, CS in mb) |
---|
252 | G4double G4QIonIonCrossSection::CalculateCrossSection(G4bool XS,G4int F,G4int I,G4int pPDG, |
---|
253 | G4int tZ,G4int tN, G4double TotMom) |
---|
254 | { |
---|
255 | //static const G4double third=1./3.; // power for A^P->R conversion [R=1.1*A^(1/3)] |
---|
256 | //static const G4double conv=38.; // coeff. R2->sig=c*(pR+tR)^2, c=pi*10(mb/fm^2)*1.21 |
---|
257 | // If change the following, please change in ::GetFunctions: |
---|
258 | static const G4double THmin=0.; // @@ start from threshold (?) minimum Energy Threshold |
---|
259 | static const G4double dP=10.; // step for the LEN table |
---|
260 | static const G4int nL=100; // A#of LENesonance points in E (each MeV from 2 to 106) |
---|
261 | static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part |
---|
262 | static const G4double Pmax=300000.; // maxE for the HighE part |
---|
263 | static const G4int nH=100; // A#of HResonance points in lnE |
---|
264 | static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part |
---|
265 | static const G4double malP=std::log(Pmax); // High logarithm energy (each 2.75 percent) |
---|
266 | static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part |
---|
267 | // |
---|
268 | // Associative memory for acceleration |
---|
269 | static std::vector <G4double*> LENI; // Vector of pointers: LowEnIneIonIonCrossSection |
---|
270 | static std::vector <G4double*> HENI; // Vector of pointers: HighEnIneIonIonCrossSection |
---|
271 | static std::vector <G4double*> LENE; // Vector of pointers: LowEnElaIonIonCrossSection |
---|
272 | static std::vector <G4double*> HENE; // Vector of pointers: HighEnElaIonIonCrossSection |
---|
273 | #ifdef debug |
---|
274 | G4cout<<"G4QIonIonCrossSection::CalcCS: Z="<<tZ<<", N="<<tN<<", P="<<TotMom<<G4endl; |
---|
275 | #endif |
---|
276 | G4int dPDG=pPDG/10; // 10SZZZAAA |
---|
277 | G4int zPDG=dPDG/1000; // 10SZZZ (?) |
---|
278 | G4int zA=dPDG%1000; // proj A |
---|
279 | G4int pZ=zPDG%1000; // proj Z (?) |
---|
280 | G4int pN=zA-pZ; // proj N (?) |
---|
281 | G4double Momentum=TotMom/zA; // Momentum per nucleon |
---|
282 | if (Momentum<THmin) return 0.; // @@ This can be dangerouse for the heaviest nuc.! |
---|
283 | G4double sigma=0.; |
---|
284 | if(F&&I) sigma=0.; // @@ *!* Fake line *!* to use F & I !!!Temporary!!! |
---|
285 | G4double tA=tN+tZ; // Target weight |
---|
286 | G4double pA=zA; // Projectile weight |
---|
287 | if(F<=0) // This isotope was not the last used isotop |
---|
288 | { |
---|
289 | if(F<0 || !XS) // This isotope was found in DAMDB or Elast =>RETRIEVE |
---|
290 | { |
---|
291 | lastLENI=LENI[I]; // Pointer to Low Energy inelastic cross sections |
---|
292 | lastHENI=HENI[I]; // Pointer to High Energy inelastic cross sections |
---|
293 | lastLENE=LENE[I]; // Pointer to Low Energy inelastic cross sections |
---|
294 | lastHENE=HENE[I]; // Pointer to High Energy inelastic cross sections |
---|
295 | } |
---|
296 | else // This isotope wasn't calculated previously => CREATE |
---|
297 | { |
---|
298 | lastLENI = new G4double[nL]; // Allocate memory for the new LEN cross sections |
---|
299 | lastHENI = new G4double[nH]; // Allocate memory for the new HEN cross sections |
---|
300 | lastLENE = new G4double[nL]; // Allocate memory for the new LEN cross sections |
---|
301 | lastHENE = new G4double[nH]; // Allocate memory for the new HEN cross sections |
---|
302 | G4int er=GetFunctions(pZ,pN,tZ,tN,lastLENI,lastHENI,lastLENE,lastHENE); |
---|
303 | if(er<1) G4cerr<<"*W*G4QIonIonCroSec::CalcCrossSection: pA="<<tA<<",tA="<<tA<<G4endl; |
---|
304 | #ifdef debug |
---|
305 | G4cout<<"G4QIonIonCrossSection::CalcCS: GetFunctions er="<<er<<",pA="<<pA<<",tA="<<tA |
---|
306 | <<G4endl; |
---|
307 | #endif |
---|
308 | // *** The synchronization check *** |
---|
309 | G4int sync=LENI.size(); |
---|
310 | if(sync!=I) G4cout<<"*W*G4IonIonCrossSec::CalcCrossSect:Sync="<<sync<<"#"<<I<<G4endl; |
---|
311 | LENI.push_back(lastLENI); // added LEN Inelastic |
---|
312 | HENI.push_back(lastHENI); // added HEN Inelastic |
---|
313 | LENE.push_back(lastLENE); // added LEN Elastic |
---|
314 | HENE.push_back(lastHENE); // added HEN Elastic |
---|
315 | } // End of creation of the new set of parameters |
---|
316 | } // End of parameters udate |
---|
317 | // ============================== NOW the Magic Formula ================================= |
---|
318 | if (Momentum<lastTH) return 0.; // It must be already checked in the interface class |
---|
319 | else if (Momentum<Pmin) // LEN region (approximated in E, not in lnE) |
---|
320 | { |
---|
321 | #ifdef debug |
---|
322 | G4cout<<"G4QIICS::CalCS:p="<<pA<<",t="<<tA<<",n="<<nL<<",T="<<THmin<<",d="<<dP<<G4endl; |
---|
323 | #endif |
---|
324 | if(tA<1. || pA<1.) |
---|
325 | { |
---|
326 | G4cout<<"-Warning-G4QIICS::CalcCS: pA="<<pA<<" or tA="<<tA<<" aren't nuclei"<<G4endl; |
---|
327 | sigma=0.; |
---|
328 | } |
---|
329 | else |
---|
330 | { |
---|
331 | G4double dPp=dP*pA; |
---|
332 | if(XS) sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENI); |
---|
333 | else sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENE); |
---|
334 | } |
---|
335 | #ifdef debugn |
---|
336 | if(sigma<0.) G4cout<<"-Warning-G4QIICS::CalcCS:pA="<<pA<<",tA="<<tA<<",XS="<<XS<<",P=" |
---|
337 | <<Momentum<<", Th="<<THmin<<", dP="<<dP<<G4endl; |
---|
338 | #endif |
---|
339 | } |
---|
340 | else if (Momentum<Pmax*pA) // High Energy region |
---|
341 | { |
---|
342 | G4double lP=std::log(Momentum); |
---|
343 | #ifdef debug |
---|
344 | G4cout<<"G4QIonIonCS::CalcCS:before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl; |
---|
345 | #endif |
---|
346 | if(tA<1. || pA<1.) |
---|
347 | { |
---|
348 | G4cout<<"-Warning-G4QIICS::CalCS:pA="<<pA<<" or tA="<<tA<<" aren't composit"<<G4endl; |
---|
349 | sigma=0.; |
---|
350 | } |
---|
351 | else |
---|
352 | { |
---|
353 | G4double milPp=milP+std::log(pA); |
---|
354 | if(XS) sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENI); |
---|
355 | else sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENE); |
---|
356 | } |
---|
357 | } |
---|
358 | else // UltraHighE region (not frequent) |
---|
359 | { |
---|
360 | std::pair<G4double, G4double> inelel = CalculateXS(pZ, pN, tZ, tN, Momentum); |
---|
361 | if(XS) sigma=inelel.first; |
---|
362 | else sigma=inelel.second; |
---|
363 | } |
---|
364 | #ifdef debug |
---|
365 | G4cout<<"G4IonIonCrossSection::CalculateCrossSection: sigma="<<sigma<<G4endl; |
---|
366 | #endif |
---|
367 | if(sigma<0.) return 0.; |
---|
368 | return sigma; |
---|
369 | } |
---|
370 | |
---|
371 | // Linear fit for YN[N] tabulated (from X0 with fixed step DX) function to X point |
---|
372 | |
---|
373 | // Calculate the functions for the log(A) |
---|
374 | G4int G4QIonIonCrossSection::GetFunctions(G4int pZ,G4int pN,G4int tZ,G4int tN,G4double* li, |
---|
375 | G4double* hi, G4double* le, G4double* he) |
---|
376 | { |
---|
377 | // If change the following, please change in ::CalculateCrossSection: |
---|
378 | static const G4double THmin=0.; // @@ start from threshold (?) minimum Energy Threshold |
---|
379 | static const G4double dP=10.; // step for the LEN table |
---|
380 | static const G4int nL=100; // A#of LENesonance points in E (each MeV from 2 to 106) |
---|
381 | static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part |
---|
382 | static const G4double Pmax=300000.; // maxE for the HighE part |
---|
383 | static const G4int nH=100; // A#of HResonance points in lnE |
---|
384 | static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part |
---|
385 | static const G4double malP=std::log(Pmax); // High logarithm energy |
---|
386 | static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part |
---|
387 | static const G4double lP=std::exp(dlP); // Multiplication factor in the HighE part |
---|
388 | // If the cross section approximation formula is changed - replace from file. |
---|
389 | if(pZ<1 || pN<0 || tZ<1 || tN<0) |
---|
390 | { |
---|
391 | G4cout<<"-W-G4QIonIonCS::GetFunct:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl; |
---|
392 | return -1; |
---|
393 | } |
---|
394 | G4int pA=pN+pZ; |
---|
395 | G4double dPp=dP*pA; |
---|
396 | G4double Mom=THmin; |
---|
397 | for(G4int k=0; k<nL; k++) |
---|
398 | { |
---|
399 | std::pair<G4double,G4double> len = CalculateXS(pZ, pN, tZ, tN, Mom); |
---|
400 | li[k]=len.first; |
---|
401 | le[k]=len.second; |
---|
402 | Mom+=dPp; |
---|
403 | } |
---|
404 | G4double lMom=Pmin*pA; |
---|
405 | for(G4int j=0; j<nH; j++) |
---|
406 | { |
---|
407 | std::pair<G4double,G4double> len = CalculateXS(pZ, pN, pZ, pN, lMom); |
---|
408 | hi[j]=len.first; |
---|
409 | he[j]=len.second; |
---|
410 | lMom*=lP; |
---|
411 | } |
---|
412 | #ifdef debug |
---|
413 | G4cout<<"G4QIonIonCS::GetFunctions: pZ="<<pZ<<", tZ="<<tZ<<" pair is calculated"<<G4endl; |
---|
414 | #endif |
---|
415 | return 1; |
---|
416 | } |
---|
417 | |
---|
418 | // Momentum (Mom=p/A) is in MeV/c, first=InelasticXS, second=ElasticXS (mb) |
---|
419 | std::pair<G4double,G4double> G4QIonIonCrossSection::CalculateXS(G4int pZ,G4int pN,G4int tZ, |
---|
420 | G4int tN, G4double Mom) |
---|
421 | { |
---|
422 | static G4VQCrossSection* ElCSman = G4QElasticCrossSection::GetPointer(); |
---|
423 | static G4VQCrossSection* InelPCSman = G4QProtonNuclearCrossSection::GetPointer(); |
---|
424 | static G4VQCrossSection* InelNCSman = G4QNeutronNuclearCrossSection::GetPointer(); |
---|
425 | G4double pA=pZ+pN; |
---|
426 | G4double tA=tZ+tN; |
---|
427 | if(pA<.9 || tA<.9 ||pA>239. || tA>239 || Mom < 0.) return std::make_pair(0.,0.); |
---|
428 | G4double inCS=0.; |
---|
429 | G4double elCS=0.; |
---|
430 | if(pA<1.1 ) // nucleon-ion interaction use NA(in,el) |
---|
431 | { |
---|
432 | if (pZ == 1 && !pN) // proton-nuclear |
---|
433 | { |
---|
434 | inCS=InelPCSman->GetCrossSection(true, Mom, tZ, tN, 2212); |
---|
435 | elCS=ElCSman->GetCrossSection(true, Mom, tZ, tN, 2212); |
---|
436 | } |
---|
437 | else if(pN == 1 && !pZ) // neutron-nuclear |
---|
438 | { |
---|
439 | inCS=InelNCSman->GetCrossSection(true, Mom, tZ, tN, 2112); |
---|
440 | elCS=ElCSman->GetCrossSection(true, Mom, tZ, tN, 2112); |
---|
441 | } |
---|
442 | else G4cerr<<"-Warn-G4QIICS::CaCS:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl; |
---|
443 | } |
---|
444 | else |
---|
445 | { |
---|
446 | G4double T=ThresholdMomentum(pZ, pN, tZ, tN); // @@ Can be cashed as lastTH (?) |
---|
447 | if(Mom<=T) |
---|
448 | { |
---|
449 | elCS=0.; |
---|
450 | inCS=0.; |
---|
451 | } |
---|
452 | else |
---|
453 | { |
---|
454 | G4double P2=Mom*Mom; |
---|
455 | G4double T2=T*T; |
---|
456 | G4double R=1.-T2/P2; // @@ Very rough threshold effect |
---|
457 | //G4double P4=P2*P2; |
---|
458 | //G4double P8=P4*P4; |
---|
459 | //G4double T4=T2*T2; |
---|
460 | //G4double tot=CalculateTotal(pA, tA, Mom)*P8/(P8+T4*T4); // @@ convert to IndepUnits |
---|
461 | G4double tot=R*CalculateTotal(pA, tA, Mom); // @@ convert to IndepUnits |
---|
462 | G4double rat=CalculateElTot(pA, tA, Mom); |
---|
463 | elCS=tot*rat; |
---|
464 | inCS=tot-elCS; |
---|
465 | } |
---|
466 | } |
---|
467 | return std::make_pair(inCS,elCS); |
---|
468 | } |
---|
469 | |
---|
470 | // Total Ion-ion cross-section (mb), Momentum (Mom) here is p/A (MeV/c=IU) (No Threshold) |
---|
471 | G4double G4QIonIonCrossSection::CalculateTotal(G4double pA, G4double tA, G4double Mom) |
---|
472 | { |
---|
473 | G4double y=std::log(Mom/1000.); // Log of momentum in GeV/c |
---|
474 | G4double ab=pA+tA; |
---|
475 | G4double al=std::log(ab); |
---|
476 | G4double ap=std::log(pA*tA); |
---|
477 | G4double e=std::pow(pA,0.1)+std::pow(tA,0.1); |
---|
478 | G4double d=e-1.55/std::pow(al,0.2); |
---|
479 | G4double f=4.; |
---|
480 | if(pA>4. && tA>4.) f=3.3+130./ab/ab+2.25/e; |
---|
481 | G4double c=(385.+11.*ab/(1.+.02*ab*al)+(.5*ab+500./al/al)/al)*std::pow(d,f); |
---|
482 | G4double r=y-3.-4./ap/ap; |
---|
483 | #ifdef pdebug |
---|
484 | G4cout<<"G4QIonIonCS::CalcTot:P="<<Mom<<", stot="<<c+d*r*r<<", d="<<d<<", r="<<r<<G4endl; |
---|
485 | #endif |
---|
486 | return c+d*r*r; |
---|
487 | } |
---|
488 | |
---|
489 | // Ratio elastic/Total, Momentum (Mom) here is p/A (MeV/c=IU) |
---|
490 | G4double G4QIonIonCrossSection::CalculateElTot(G4double pA, G4double tA, G4double Mom) |
---|
491 | { |
---|
492 | G4double y=std::log(Mom/1000.); // Log of momentum in GeV/c |
---|
493 | G4double ab=pA*tA; |
---|
494 | G4double ap=std::log(ab); |
---|
495 | G4double r=y-3.92-1.73/ap/ap; |
---|
496 | G4double d=.00166/(1.+.002*std::pow(ab,1.33333)); |
---|
497 | G4double al=std::log(pA+tA); |
---|
498 | G4double e=1.+.235*(std::fabs(ap-5.78)); |
---|
499 | if (std::fabs(pA-2.)<.1 && std::fabs(tA-2.)<.1) e=2.18; |
---|
500 | else if(std::fabs(pA-2.)<.1 && std::fabs(tA-3.)<.1) e=1.90; |
---|
501 | else if(std::fabs(pA-3.)<.1 && std::fabs(tA-2.)<.1) e=1.90; |
---|
502 | else if(std::fabs(pA-2.)<.1 && std::fabs(tA-4.)<.1) e=1.65; |
---|
503 | else if(std::fabs(pA-4.)<.1 && std::fabs(tA-2.)<.1) e=1.65; |
---|
504 | else if(std::fabs(pA-3.)<.1 && std::fabs(tA-4.)<.1) e=1.32; |
---|
505 | else if(std::fabs(pA-4.)<.1 && std::fabs(tA-3.)<.1) e=1.32; |
---|
506 | else if(std::fabs(pA-4.)<.1 && std::fabs(tA-4.)<.1) e=1.; |
---|
507 | G4double f=.37+.0188*al; |
---|
508 | G4double g=std::log(std::pow(pA,0.35)+std::pow(tA,0.35)); |
---|
509 | G4double h=g*g; |
---|
510 | G4double c=f/(1.+e/h/h); |
---|
511 | #ifdef pdebug |
---|
512 | G4cout<<"G4QIonIonCS::CalcElT:P="<<Mom<<",el/tot="<<c+d*r*r<<",d="<<d<<", r="<<r<<G4endl; |
---|
513 | #endif |
---|
514 | return c+d*r*r; |
---|
515 | } |
---|
516 | |
---|
517 | // Electromagnetic momentum/A-threshold (in MeV/c) |
---|
518 | G4double G4QIonIonCrossSection::ThresholdMomentum(G4int pZ, G4int pN, G4int tZ, G4int tN) |
---|
519 | { |
---|
520 | static const G4double third=1./3.; |
---|
521 | static const G4double pM = G4QPDGCode(2212).GetMass(); // Proton mass in MeV |
---|
522 | static const G4double tpM= pM+pM; // Doubled proton mass (MeV) |
---|
523 | if(pZ<.99 || pN<0. || tZ<.99 || tN<0.) return 0.; |
---|
524 | G4double tA=tZ+tN; |
---|
525 | G4double pA=pZ+pN; |
---|
526 | //G4double dE=1.263*tZ/(1.+std::pow(tA,third)); |
---|
527 | G4double dE=pZ*tZ/(std::pow(pA,third)+std::pow(tA,third))/pA; // dE/pA (per projNucleon) |
---|
528 | return std::sqrt(dE*(tpM+dE)); |
---|
529 | } |
---|