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 | // $Id: G4QuasiFreeRatios.cc,v 1.3 2010/01/22 17:02:48 mkossov Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ |
---|
29 | // |
---|
30 | // |
---|
31 | // G4 Physics class: G4QuasiFreeRatios for N+A elastic cross sections |
---|
32 | // Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01 |
---|
33 | // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Oct-06 |
---|
34 | // |
---|
35 | //======================================================================= |
---|
36 | // Short description: Provides percentage of quasi-free and quasi-elastic |
---|
37 | // reactions in the inelastic reactions. |
---|
38 | // ---------------------------------------------------------------------- |
---|
39 | |
---|
40 | //#define debug |
---|
41 | //#define pdebug |
---|
42 | //#define ppdebug |
---|
43 | //#define nandebug |
---|
44 | |
---|
45 | #include "G4QuasiFreeRatios.hh" |
---|
46 | |
---|
47 | // initialisation of statics |
---|
48 | std::vector<G4double*> G4QuasiFreeRatios::vT; // Vector of pointers to LinTable in C++ heap |
---|
49 | std::vector<G4double*> G4QuasiFreeRatios::vL; // Vector of pointers to LogTable in C++ heap |
---|
50 | std::vector<std::pair<G4double,G4double>*> G4QuasiFreeRatios::vX; // ETPointers to LogTable |
---|
51 | |
---|
52 | G4QuasiFreeRatios::G4QuasiFreeRatios() |
---|
53 | { |
---|
54 | #ifdef pdebug |
---|
55 | G4cout<<"***^^^*** G4QuasiFreeRatios singletone is created ***^^^***"<<G4endl; |
---|
56 | #endif |
---|
57 | } |
---|
58 | |
---|
59 | G4QuasiFreeRatios::~G4QuasiFreeRatios() |
---|
60 | { |
---|
61 | std::vector<G4double*>::iterator pos; |
---|
62 | for(pos=vT.begin(); pos<vT.end(); pos++) |
---|
63 | { delete [] *pos; } |
---|
64 | vT.clear(); |
---|
65 | for(pos=vL.begin(); pos<vL.end(); pos++) |
---|
66 | { delete [] *pos; } |
---|
67 | vL.clear(); |
---|
68 | |
---|
69 | std::vector<std::pair<G4double,G4double>*>::iterator pos2; |
---|
70 | for(pos2=vX.begin(); pos2<vX.end(); pos2++) |
---|
71 | { delete [] *pos2; } |
---|
72 | vX.clear(); |
---|
73 | } |
---|
74 | |
---|
75 | // Returns Pointer to the G4VQCrossSection class |
---|
76 | G4QuasiFreeRatios* G4QuasiFreeRatios::GetPointer() |
---|
77 | { |
---|
78 | static G4QuasiFreeRatios theRatios; // *** Static body of the QEl Cross Section *** |
---|
79 | return &theRatios; |
---|
80 | } |
---|
81 | |
---|
82 | // Calculation of pair(QuasiFree/Inelastic,QuasiElastic/QuasiFree) |
---|
83 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetRatios(G4double pIU, G4int pPDG, |
---|
84 | G4int tgZ, G4int tgN) |
---|
85 | { |
---|
86 | #ifdef pdebug |
---|
87 | G4cout<<">>>IN>>>G4QFRat::GetQF:P="<<pIU<<",pPDG="<<pPDG<<",Z="<<tgZ<<",N="<<tgN<<G4endl; |
---|
88 | #endif |
---|
89 | G4double R=0.; |
---|
90 | G4double QF2In=1.; // Prototype of QuasiFree/Inel ratio for hN_tot |
---|
91 | G4int tgA=tgZ+tgN; |
---|
92 | if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon |
---|
93 | std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean hN El&Tot(IU) |
---|
94 | //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE |
---|
95 | if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.; // To accelerate @lowE |
---|
96 | else if(ElTot.second>0.) |
---|
97 | { |
---|
98 | R=ElTot.first/ElTot.second; // El/Total ratio (does not depend on units |
---|
99 | QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN); // QuasiFree/Inelastic ratio |
---|
100 | } |
---|
101 | #ifdef pdebug |
---|
102 | G4cout<<">>>OUT>>>G4QuasiFreeRatio::GetQF2IN_Ratio: QF2In="<<QF2In<<", R="<<R<<G4endl; |
---|
103 | #endif |
---|
104 | return std::make_pair(QF2In,R); |
---|
105 | } |
---|
106 | |
---|
107 | // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section (mb) and A |
---|
108 | G4double G4QuasiFreeRatios::GetQF2IN_Ratio(G4double s, G4int A) |
---|
109 | { |
---|
110 | static const G4int nps=150; // Number of steps in the R(s) LinTable |
---|
111 | static const G4int mps=nps+1; // Number of elements in the R(s) LinTable |
---|
112 | static const G4double sma=150.; // The first LinTabEl(s=0)=1., s>sma -> logTab |
---|
113 | static const G4double ds=sma/nps; // Step of the linear Table |
---|
114 | static const G4int nls=100; // Number of steps in the R(lns) logTable |
---|
115 | static const G4int mls=nls+1; // Number of elements in the R(lns) logTable |
---|
116 | static const G4double lsi=5.; // The min ln(s) logTabEl(s=148.4 < sma=150.) |
---|
117 | static const G4double lsa=9.; // The max ln(s) logTabEl(s=148.4 - 8103. mb) |
---|
118 | static const G4double mi=std::exp(lsi);// The min s of logTabEl(~ 148.4 mb) |
---|
119 | static const G4double ms=std::exp(lsa);// The max s of logTabEl(~ 8103. mb) |
---|
120 | static const G4double dl=(lsa-lsi)/nls;// Step of the logarithmic Table |
---|
121 | static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table |
---|
122 | static const G4double toler=.01; // The tolarence mb defining the same cross-section |
---|
123 | static G4double lastS=0.; // The last sigma value for which R was calculated |
---|
124 | static G4double lastR=0.; // The last ratio R which was calculated |
---|
125 | // Local Associative Data Base: |
---|
126 | static std::vector<G4int> vA; // Vector of calculated A |
---|
127 | static std::vector<G4double> vH; // Vector of max s initialized in the LinTable |
---|
128 | static std::vector<G4int> vN; // Vector of topBin number initialized in LinTable |
---|
129 | static std::vector<G4double> vM; // Vector of rel max ln(s) initialized in LogTable |
---|
130 | static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable |
---|
131 | // Last values of the Associative Data Base: |
---|
132 | static G4int lastA=0; // theLast of calculated A |
---|
133 | static G4double lastH=0.; // theLast of max s initialized in the LinTable |
---|
134 | static G4int lastN=0; // theLast of topBin number initialized in LinTable |
---|
135 | static G4double lastM=0.; // theLast of rel max ln(s) initialized in LogTable |
---|
136 | static G4int lastK=0; // theLast of topBin number initialized in LogTable |
---|
137 | static G4double* lastT=0; // theLast of pointer to LinTable in the C++ heap |
---|
138 | static G4double* lastL=0; // theLast of pointer to LogTable in the C++ heap |
---|
139 | // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei |
---|
140 | #ifdef pdebug |
---|
141 | G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<", s="<<s<<G4endl; |
---|
142 | #endif |
---|
143 | if(s<toler || A<2) return 1.; |
---|
144 | if(s>ms) return 0.; |
---|
145 | if(A>238) |
---|
146 | { |
---|
147 | G4cout<<"-Warning-G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl; |
---|
148 | return 0.; |
---|
149 | } |
---|
150 | G4int nDB=vA.size(); // A number of nuclei already initialized in AMDB |
---|
151 | // if(nDB && lastA==A && std::fabs(s-lastS)<toler) return lastR; |
---|
152 | if(nDB && lastA==A && s==lastS) return lastR; // VI do not use tolerance |
---|
153 | G4bool found=false; |
---|
154 | G4int i=-1; |
---|
155 | if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Sirch for this A in AMDB |
---|
156 | { |
---|
157 | found=true; // The A value is found |
---|
158 | break; |
---|
159 | } |
---|
160 | #ifdef pdebug |
---|
161 | G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio: nDB="<<nDB<<", found="<<found<<G4endl; |
---|
162 | #endif |
---|
163 | if(!nDB || !found) // Create new line in the AMDB |
---|
164 | { |
---|
165 | lastA = A; |
---|
166 | #ifdef pdebug |
---|
167 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewT, A="<<A<<", nDB="<<nDB<<G4endl; |
---|
168 | #endif |
---|
169 | lastT = new G4double[mps]; // Create the linear Table |
---|
170 | lastN = static_cast<int>(s/ds)+1; // MaxBin to be initialized |
---|
171 | if(lastN>nps) |
---|
172 | { |
---|
173 | lastN=nps; |
---|
174 | lastH=sma; |
---|
175 | } |
---|
176 | else lastH = lastN*ds; // Calculate max initialized s for LinTab |
---|
177 | G4double sv=0; |
---|
178 | lastT[0]=1.; |
---|
179 | for(G4int j=1; j<=lastN; j++) // Calculate LogTab values |
---|
180 | { |
---|
181 | sv+=ds; |
---|
182 | lastT[j]=CalcQF2IN_Ratio(sv,A); |
---|
183 | } |
---|
184 | lastL=new G4double[mls]; // Create the logarithmic Table |
---|
185 | if(s>sma) // Initialize the logarithmic Table |
---|
186 | { |
---|
187 | #ifdef pdebug |
---|
188 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewL, A="<<A<<", nDB="<<nDB<<G4endl; |
---|
189 | #endif |
---|
190 | G4double ls=std::log(s); |
---|
191 | lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB |
---|
192 | if(lastK>nls) |
---|
193 | { |
---|
194 | lastK=nls; |
---|
195 | lastM=lsa-lsi; |
---|
196 | } |
---|
197 | else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab |
---|
198 | sv=mi; |
---|
199 | for(G4int j=0; j<=lastK; j++) // Calculate LogTab values |
---|
200 | { |
---|
201 | lastL[j]=CalcQF2IN_Ratio(sv,A); |
---|
202 | if(j!=lastK) sv*=edl; |
---|
203 | } |
---|
204 | } |
---|
205 | else // LogTab is not initialized |
---|
206 | { |
---|
207 | lastK = 0; |
---|
208 | lastM = 0.; |
---|
209 | } |
---|
210 | i++; // Make a new record to AMDB and position on it |
---|
211 | vA.push_back(lastA); |
---|
212 | vH.push_back(lastH); |
---|
213 | vN.push_back(lastN); |
---|
214 | vM.push_back(lastM); |
---|
215 | vK.push_back(lastK); |
---|
216 | vT.push_back(lastT); |
---|
217 | vL.push_back(lastL); |
---|
218 | } |
---|
219 | else // The A value was found in AMDB |
---|
220 | { |
---|
221 | lastA=vA[i]; |
---|
222 | lastH=vH[i]; |
---|
223 | lastN=vN[i]; |
---|
224 | lastM=vM[i]; |
---|
225 | lastK=vK[i]; |
---|
226 | lastT=vT[i]; |
---|
227 | lastL=vL[i]; |
---|
228 | #ifdef pdebug |
---|
229 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: Found, s="<<s<<", lastM="<<lastM<<G4endl; |
---|
230 | #endif |
---|
231 | if(s>lastM) // At least LinTab must be updated |
---|
232 | { |
---|
233 | G4int nextN=lastN+1; // The next bin to be initialized |
---|
234 | #ifdef pdebug |
---|
235 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: lastN="<<lastN<<" ?< nps="<<nps<<G4endl; |
---|
236 | #endif |
---|
237 | if(lastN<nps) |
---|
238 | { |
---|
239 | lastN = static_cast<int>(s/ds)+1;// MaxBin to be initialized |
---|
240 | if(lastN>nps) |
---|
241 | { |
---|
242 | lastN=nps; |
---|
243 | lastH=sma; |
---|
244 | } |
---|
245 | else lastH = lastN*ds; // Calculate max initialized s for LinTab |
---|
246 | G4double sv=lastM; |
---|
247 | for(G4int j=nextN; j<=lastN; j++)// Calculate LogTab values |
---|
248 | { |
---|
249 | sv+=ds; |
---|
250 | lastT[j]=CalcQF2IN_Ratio(sv,A); |
---|
251 | } |
---|
252 | #ifdef pdebug |
---|
253 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl; |
---|
254 | #endif |
---|
255 | } // End of LinTab update |
---|
256 | #ifdef pdebug |
---|
257 | G4cout<<"G4QFRatios::GetQF2IN_Ratio: lN="<<lastN<<", nN="<<nextN<<", i="<<i<<G4endl; |
---|
258 | #endif |
---|
259 | if(lastN>=nextN) |
---|
260 | { |
---|
261 | vH[i]=lastH; |
---|
262 | vN[i]=lastN; |
---|
263 | } |
---|
264 | G4int nextK=lastK+1; |
---|
265 | if(!lastK) nextK=0; |
---|
266 | #ifdef pdebug |
---|
267 | G4cout<<"G4QFRat::GetQF2IN_Ratio: sma="<<sma<<", lastK="<<lastK<<" < "<<nls<<G4endl; |
---|
268 | #endif |
---|
269 | if(s>sma && lastK<nls) // LogTab must be updated |
---|
270 | { |
---|
271 | G4double sv=std::exp(lastM+lsi); // Define starting poit (lastM will be changed) |
---|
272 | G4double ls=std::log(s); |
---|
273 | lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB |
---|
274 | if(lastK>nls) |
---|
275 | { |
---|
276 | lastK=nls; |
---|
277 | lastM=lsa-lsi; |
---|
278 | } |
---|
279 | else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab |
---|
280 | #ifdef pdebug |
---|
281 | G4cout<<"G4QFRat::GetQF2IN_Ratio: nK="<<nextK<<", lK="<<lastK<<", sv="<<sv<<G4endl; |
---|
282 | #endif |
---|
283 | for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values |
---|
284 | { |
---|
285 | sv*=edl; |
---|
286 | #ifdef pdebug |
---|
287 | G4cout<<"G4QFRat::GetQF2IN_Ratio: j="<<j<<", sv="<<sv<<", A="<<A<<G4endl; |
---|
288 | #endif |
---|
289 | lastL[j]=CalcQF2IN_Ratio(sv,A); |
---|
290 | } |
---|
291 | #ifdef pdebug |
---|
292 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl; |
---|
293 | #endif |
---|
294 | } // End of LogTab update |
---|
295 | #ifdef pdebug |
---|
296 | G4cout<<"G4QFRatios::GetQF2IN_Ratio: lK="<<lastK<<", nK="<<nextK<<", i="<<i<<G4endl; |
---|
297 | #endif |
---|
298 | if(lastK>=nextK) |
---|
299 | { |
---|
300 | vM[i]=lastM; |
---|
301 | vK[i]=lastK; |
---|
302 | } |
---|
303 | } |
---|
304 | } |
---|
305 | #ifdef pdebug |
---|
306 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeTab s="<<s<<", sma="<<sma<<G4endl; |
---|
307 | #endif |
---|
308 | // Now one can use tabeles to calculate the value |
---|
309 | if(s<sma) // Use linear table |
---|
310 | { |
---|
311 | G4int n=static_cast<int>(s/ds); // Low edge number of the bin |
---|
312 | G4double d=s-n*ds; // Linear shift |
---|
313 | G4double v=lastT[n]; // Base |
---|
314 | lastR=v+d*(lastT[n+1]-v)/ds; // Result |
---|
315 | } |
---|
316 | else // Use log table |
---|
317 | { |
---|
318 | G4double ls=std::log(s)-lsi; // ln(s)-l_min |
---|
319 | G4int n=static_cast<int>(ls/dl); // Low edge number of the bin |
---|
320 | G4double d=ls-n*dl; // Log shift |
---|
321 | G4double v=lastL[n]; // Base |
---|
322 | lastR=v+d*(lastL[n+1]-v)/dl; // Result |
---|
323 | } |
---|
324 | if(lastR<0.) lastR=0.; |
---|
325 | if(lastR>1.) lastR=1.; |
---|
326 | #ifdef pdebug |
---|
327 | G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeRet lastR="<<lastR<<G4endl; |
---|
328 | #endif |
---|
329 | return lastR; |
---|
330 | } // End of CalcQF2IN_Ratio |
---|
331 | |
---|
332 | // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section and A |
---|
333 | G4double G4QuasiFreeRatios::CalcQF2IN_Ratio(G4double s, G4int A) |
---|
334 | { |
---|
335 | static const G4double C=1.246; |
---|
336 | G4double s2=s*s; |
---|
337 | G4double s4=s2*s2; |
---|
338 | G4double ss=std::sqrt(std::sqrt(s)); |
---|
339 | G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2); |
---|
340 | G4double E=.2644+.016/(1.+std::exp((29.54-s)/2.49)); |
---|
341 | G4double F=ss*.1526*std::exp(-s2*ss*.0000859); |
---|
342 | return C*std::exp(-E*std::pow(G4double(A-1.),F))/std::pow(G4double(A),P); |
---|
343 | } // End of CalcQF2IN_Ratio |
---|
344 | |
---|
345 | // Calculatio pair(hN_el,hN_tot) (mb): p in GeV/c, index(PDG,F) (see FetchElTot) |
---|
346 | std::pair<G4double,G4double> G4QuasiFreeRatios::CalcElTot(G4double p, G4int I) |
---|
347 | { |
---|
348 | // ---------> Each parameter set can have not more than nPoints=128 parameters |
---|
349 | static const G4double lmi=3.5; // min of (lnP-lmi)^2 parabola |
---|
350 | static const G4double pbe=.0557; // elastic (lnP-lmi)^2 parabola coefficient |
---|
351 | static const G4double pbt=.3; // total (lnP-lmi)^2 parabola coefficient |
---|
352 | static const G4double pmi=.1; // Below that fast LE calculation is made |
---|
353 | static const G4double pma=1000.; // Above that fast HE calculation is made |
---|
354 | // ====================================================================================== |
---|
355 | G4double El=0.; // prototype of the elastic hN cross-section |
---|
356 | G4double To=0.; // prototype of the total hN cross-section |
---|
357 | if(p<=0.) |
---|
358 | { |
---|
359 | G4cout<<"-Warning-G4QuasiFreeRatios::CalcElTot: p="<<p<<" is zero or negative"<<G4endl; |
---|
360 | return std::make_pair(El,To); |
---|
361 | } |
---|
362 | if (!I) // pp/nn |
---|
363 | { |
---|
364 | #ifdef debug |
---|
365 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0, p="<<p<<", pmi="<<pmi<<", pma="<<pma<<G4endl; |
---|
366 | #endif |
---|
367 | if(p<pmi) |
---|
368 | { |
---|
369 | G4double p2=p*p; |
---|
370 | El=1./(.00012+p2*.2); |
---|
371 | To=El; |
---|
372 | #ifdef debug |
---|
373 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0i, El="<<El<<", To="<<To<<", p2="<<p2<<G4endl; |
---|
374 | #endif |
---|
375 | } |
---|
376 | else if(p>pma) |
---|
377 | { |
---|
378 | G4double lp=std::log(p)-lmi; |
---|
379 | G4double lp2=lp*lp; |
---|
380 | El=pbe*lp2+6.72; |
---|
381 | To=pbt*lp2+38.2; |
---|
382 | #ifdef debug |
---|
383 | G4cout<<"G4QuasiFreeR::CalcElTot:I=0a, El="<<El<<", To="<<To<<", lp2="<<lp2<<G4endl; |
---|
384 | #endif |
---|
385 | } |
---|
386 | else |
---|
387 | { |
---|
388 | G4double p2=p*p; |
---|
389 | G4double LE=1./(.00012+p2*.2); |
---|
390 | G4double lp=std::log(p)-lmi; |
---|
391 | G4double lp2=lp*lp; |
---|
392 | G4double rp2=1./p2; |
---|
393 | El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p); |
---|
394 | To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2); |
---|
395 | #ifdef debug |
---|
396 | G4cout<<"G4QuasiFreeR::CalcElTot:0,E="<<El<<",T="<<To<<",s="<<p2<<",l="<<lp2<<G4endl; |
---|
397 | #endif |
---|
398 | } |
---|
399 | } |
---|
400 | else if(I==1) // np/pn |
---|
401 | { |
---|
402 | if(p<pmi) |
---|
403 | { |
---|
404 | G4double p2=p*p; |
---|
405 | El=1./(.00012+p2*(.051+.1*p2)); |
---|
406 | To=El; |
---|
407 | } |
---|
408 | else if(p>pma) |
---|
409 | { |
---|
410 | G4double lp=std::log(p)-lmi; |
---|
411 | G4double lp2=lp*lp; |
---|
412 | El=pbe*lp2+6.72; |
---|
413 | To=pbt*lp2+38.2; |
---|
414 | } |
---|
415 | else |
---|
416 | { |
---|
417 | G4double p2=p*p; |
---|
418 | G4double LE=1./(.00012+p2*(.051+.1*p2)); |
---|
419 | G4double lp=std::log(p)-lmi; |
---|
420 | G4double lp2=lp*lp; |
---|
421 | G4double rp2=1./p2; |
---|
422 | El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p); |
---|
423 | To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2); |
---|
424 | } |
---|
425 | } |
---|
426 | else if(I==2) // pimp/pipn |
---|
427 | { |
---|
428 | G4double lp=std::log(p); |
---|
429 | if(p<pmi) |
---|
430 | { |
---|
431 | G4double lr=lp+1.27; |
---|
432 | El=1.53/(lr*lr+.0676); |
---|
433 | To=El*3; |
---|
434 | } |
---|
435 | else if(p>pma) |
---|
436 | { |
---|
437 | G4double ld=lp-lmi; |
---|
438 | G4double ld2=ld*ld; |
---|
439 | G4double sp=std::sqrt(p); |
---|
440 | El=pbe*ld2+2.4+7./sp; |
---|
441 | To=pbt*ld2+22.3+12./sp; |
---|
442 | } |
---|
443 | else |
---|
444 | { |
---|
445 | G4double lr=lp+1.27; // p1 |
---|
446 | G4double LE=1.53/(lr*lr+.0676); // p2, p3 |
---|
447 | G4double ld=lp-lmi; // p4 (lmi=3.5) |
---|
448 | G4double ld2=ld*ld; |
---|
449 | G4double p2=p*p; |
---|
450 | G4double p4=p2*p2; |
---|
451 | G4double sp=std::sqrt(p); |
---|
452 | G4double lm=lp+.36; // p5 |
---|
453 | G4double md=lm*lm+.04; // p6 |
---|
454 | G4double lh=lp-.017; // p7 |
---|
455 | G4double hd=lh*lh+.0025; // p8 |
---|
456 | El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14 |
---|
457 | To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd; |
---|
458 | } |
---|
459 | } |
---|
460 | else if(I==3) // pipp/pimn |
---|
461 | { |
---|
462 | G4double lp=std::log(p); |
---|
463 | if(p<pmi) |
---|
464 | { |
---|
465 | G4double lr=lp+1.27; |
---|
466 | G4double lr2=lr*lr; |
---|
467 | El=13./(lr2+lr2*lr2+.0676); |
---|
468 | To=El; |
---|
469 | } |
---|
470 | else if(p>pma) |
---|
471 | { |
---|
472 | G4double ld=lp-lmi; |
---|
473 | G4double ld2=ld*ld; |
---|
474 | G4double sp=std::sqrt(p); |
---|
475 | El=pbe*ld2+2.4+6./sp; |
---|
476 | To=pbt*ld2+22.3+5./sp; |
---|
477 | } |
---|
478 | else |
---|
479 | { |
---|
480 | G4double lr=lp+1.27; // p1 |
---|
481 | G4double lr2=lr*lr; |
---|
482 | G4double LE=13./(lr2+lr2*lr2+.0676); // p2, p3 |
---|
483 | G4double ld=lp-lmi; // p4 (lmi=3.5) |
---|
484 | G4double ld2=ld*ld; |
---|
485 | G4double p2=p*p; |
---|
486 | G4double p4=p2*p2; |
---|
487 | G4double sp=std::sqrt(p); |
---|
488 | G4double lm=lp-.32; // p5 |
---|
489 | G4double md=lm*lm+.0576; // p6 |
---|
490 | El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11 |
---|
491 | To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md; |
---|
492 | } |
---|
493 | } |
---|
494 | else if(I==4) // Kmp/Kmn/K0p/K0n |
---|
495 | { |
---|
496 | |
---|
497 | if(p<pmi) |
---|
498 | { |
---|
499 | G4double psp=p*std::sqrt(p); |
---|
500 | El=5.2/psp; |
---|
501 | To=14./psp; |
---|
502 | } |
---|
503 | else if(p>pma) |
---|
504 | { |
---|
505 | G4double ld=std::log(p)-lmi; |
---|
506 | G4double ld2=ld*ld; |
---|
507 | El=pbe*ld2+2.23; |
---|
508 | To=pbt*ld2+19.5; |
---|
509 | } |
---|
510 | else |
---|
511 | { |
---|
512 | G4double ld=std::log(p)-lmi; |
---|
513 | G4double ld2=ld*ld; |
---|
514 | G4double sp=std::sqrt(p); |
---|
515 | G4double psp=p*sp; |
---|
516 | G4double p2=p*p; |
---|
517 | G4double p4=p2*p2; |
---|
518 | G4double lm=p-.39; |
---|
519 | G4double md=lm*lm+.000156; |
---|
520 | G4double lh=p-1.; |
---|
521 | G4double hd=lh*lh+.0156; |
---|
522 | El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd; |
---|
523 | To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd; |
---|
524 | } |
---|
525 | } |
---|
526 | else if(I==5) // Kpp/Kpn/aKp/aKn |
---|
527 | { |
---|
528 | if(p<pmi) |
---|
529 | { |
---|
530 | G4double lr=p-.38; |
---|
531 | G4double lm=p-1.; |
---|
532 | G4double md=lm*lm+.372; |
---|
533 | El=.7/(lr*lr+.0676)+2./md; |
---|
534 | To=El+.6/md; |
---|
535 | } |
---|
536 | else if(p>pma) |
---|
537 | { |
---|
538 | G4double ld=std::log(p)-lmi; |
---|
539 | G4double ld2=ld*ld; |
---|
540 | El=pbe*ld2+2.23; |
---|
541 | To=pbt*ld2+19.5; |
---|
542 | } |
---|
543 | else |
---|
544 | { |
---|
545 | G4double ld=std::log(p)-lmi; |
---|
546 | G4double ld2=ld*ld; |
---|
547 | G4double lr=p-.38; |
---|
548 | G4double LE=.7/(lr*lr+.0676); |
---|
549 | G4double sp=std::sqrt(p); |
---|
550 | G4double p2=p*p; |
---|
551 | G4double p4=p2*p2; |
---|
552 | G4double lm=p-1.; |
---|
553 | G4double md=lm*lm+.372; |
---|
554 | El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md; |
---|
555 | To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md; |
---|
556 | } |
---|
557 | } |
---|
558 | else if(I==6) // hyperon-N |
---|
559 | { |
---|
560 | if(p<pmi) |
---|
561 | { |
---|
562 | G4double p2=p*p; |
---|
563 | El=1./(.002+p2*(.12+p2)); |
---|
564 | To=El; |
---|
565 | } |
---|
566 | else if(p>pma) |
---|
567 | { |
---|
568 | G4double lp=std::log(p)-lmi; |
---|
569 | G4double lp2=lp*lp; |
---|
570 | G4double sp=std::sqrt(p); |
---|
571 | El=(pbe*lp2+6.72)/(1.+2./sp); |
---|
572 | To=(pbt*lp2+38.2+900./sp)/(1.+27./sp); |
---|
573 | } |
---|
574 | else |
---|
575 | { |
---|
576 | G4double p2=p*p; |
---|
577 | G4double LE=1./(.002+p2*(.12+p2)); |
---|
578 | G4double lp=std::log(p)-lmi; |
---|
579 | G4double lp2=lp*lp; |
---|
580 | G4double p4=p2*p2; |
---|
581 | G4double sp=std::sqrt(p); |
---|
582 | El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4); |
---|
583 | To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4); |
---|
584 | } |
---|
585 | } |
---|
586 | else if(I==7) // antibaryon-N |
---|
587 | { |
---|
588 | if(p>pma) |
---|
589 | { |
---|
590 | G4double lp=std::log(p)-lmi; |
---|
591 | G4double lp2=lp*lp; |
---|
592 | El=pbe*lp2+6.72; |
---|
593 | To=pbt*lp2+38.2; |
---|
594 | } |
---|
595 | else |
---|
596 | { |
---|
597 | G4double ye=std::pow(p,1.25); |
---|
598 | G4double yt=std::pow(p,.35); |
---|
599 | G4double lp=std::log(p)-lmi; |
---|
600 | G4double lp2=lp*lp; |
---|
601 | El=80./(ye+1.)+pbe*lp2+6.72; |
---|
602 | To=(80./yt+.3)/yt+pbt*lp2+38.2; |
---|
603 | } |
---|
604 | } |
---|
605 | else |
---|
606 | { |
---|
607 | G4cout<<"*Error*G4QuasiFreeRatios::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl; |
---|
608 | G4Exception("G4QuasiFreeRatios::CalcElTot:","23",FatalException,"CHIPScrash"); |
---|
609 | } |
---|
610 | if(El>To) El=To; |
---|
611 | return std::make_pair(El,To); |
---|
612 | } // End of CalcElTot |
---|
613 | |
---|
614 | // For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb) |
---|
615 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTotXS(G4double p, G4int PDG, G4bool F) |
---|
616 | { |
---|
617 | G4int ind=0; // Prototype of the reaction index |
---|
618 | G4bool kfl=true; // Flag of K0/aK0 oscillation |
---|
619 | G4bool kf=false; |
---|
620 | if(PDG==130||PDG==310) |
---|
621 | { |
---|
622 | kf=true; |
---|
623 | if(G4UniformRand()>.5) kfl=false; |
---|
624 | } |
---|
625 | if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn |
---|
626 | else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn |
---|
627 | else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn |
---|
628 | else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn |
---|
629 | else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N |
---|
630 | else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N |
---|
631 | else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda |
---|
632 | else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n) |
---|
633 | else { |
---|
634 | G4cout<<"*Error*G4QuasiFreeRatios::CalcElTotXS: PDG="<<PDG |
---|
635 | <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; |
---|
636 | G4Exception("G4QuasiFreeRatio::CalcElTotXS:","22",FatalException,"CHIPScrash"); |
---|
637 | } |
---|
638 | return CalcElTot(p,ind); |
---|
639 | } |
---|
640 | |
---|
641 | // Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron |
---|
642 | std::pair<G4double,G4double> G4QuasiFreeRatios::FetchElTot(G4double p, G4int PDG, G4bool F) |
---|
643 | { |
---|
644 | static const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step) |
---|
645 | static const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable |
---|
646 | static const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) |
---|
647 | static const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) |
---|
648 | static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c) |
---|
649 | static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV) |
---|
650 | static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table |
---|
651 | static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table |
---|
652 | //static const G4double toler=.001; // Relative Tolarence defining "theSameMomentum" |
---|
653 | static G4double lastP=0.; // The last momentum for which XS was calculated |
---|
654 | static G4int lastH=0; // The last projPDG for which XS was calculated |
---|
655 | static G4bool lastF=true; // The last nucleon for which XS was calculated |
---|
656 | static std::pair<G4double,G4double> lastR=std::make_pair(0.,0.); // The last result |
---|
657 | // Local Associative Data Base: |
---|
658 | static std::vector<G4int> vI; // Vector of index for which XS was calculated |
---|
659 | static std::vector<G4double> vM; // Vector of rel max ln(p) initialized in LogTable |
---|
660 | static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable |
---|
661 | // Last values of the Associative Data Base: |
---|
662 | static G4int lastI=0; // The Last index for which XS was calculated |
---|
663 | static G4double lastM=0.; // The Last rel max ln(p) initialized in LogTable |
---|
664 | static G4int lastK=0; // The Last topBin number initialized in LogTable |
---|
665 | static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap |
---|
666 | // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei |
---|
667 | G4int nDB=vI.size(); // A number of hadrons already initialized in AMDB |
---|
668 | #ifdef pdebug |
---|
669 | G4cout<<"G4QuasiFreeR::FetchElTot:p="<<p<<",PDG="<<PDG<<",F="<<F<<",nDB="<<nDB<<G4endl; |
---|
670 | #endif |
---|
671 | if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler. |
---|
672 | // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR; |
---|
673 | lastH=PDG; |
---|
674 | lastF=F; |
---|
675 | G4int ind=-1; // Prototipe of the index of the PDG/F combination |
---|
676 | // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p), |
---|
677 | // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp) |
---|
678 | G4bool kfl=true; // Flag of K0/aK0 oscillation |
---|
679 | G4bool kf=false; |
---|
680 | if(PDG==130||PDG==310) |
---|
681 | { |
---|
682 | kf=true; |
---|
683 | if(G4UniformRand()>.5) kfl=false; |
---|
684 | } |
---|
685 | if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn |
---|
686 | else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn |
---|
687 | else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn |
---|
688 | else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn |
---|
689 | else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N |
---|
690 | else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N |
---|
691 | else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda |
---|
692 | else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n) |
---|
693 | else { |
---|
694 | G4cout<<"*Error*G4QuasiFreeRatios::FetchElTot: PDG="<<PDG |
---|
695 | <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; |
---|
696 | G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash"); |
---|
697 | } |
---|
698 | if(nDB && lastI==ind && p>0. && p==lastP) return lastR; // VI do not use toler |
---|
699 | // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR; |
---|
700 | if(p<=mi || p>=ma) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?) |
---|
701 | G4bool found=false; |
---|
702 | G4int i=-1; |
---|
703 | if(nDB) for (i=0; i<nDB; i++) if(ind==vI[i]) // Sirch for this index in AMDB |
---|
704 | { |
---|
705 | found=true; // The index is found |
---|
706 | break; |
---|
707 | } |
---|
708 | G4double lp=std::log(p); |
---|
709 | #ifdef pdebug |
---|
710 | G4cout<<"G4QuasiFreeR::FetchElTot:I="<<ind<<",i="<<i<<",fd="<<found<<",lp="<<lp<<G4endl; |
---|
711 | #endif |
---|
712 | if(!nDB || !found) // Create new line in the AMDB |
---|
713 | { |
---|
714 | #ifdef pdebug |
---|
715 | G4cout<<"G4QuasiFreeRatios::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl; |
---|
716 | #endif |
---|
717 | lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot |
---|
718 | lastI = ind; // Remember the initialized inex |
---|
719 | lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTaB |
---|
720 | if(lastK>nlp) |
---|
721 | { |
---|
722 | lastK=nlp; |
---|
723 | lastM=lpa-lpi; |
---|
724 | } |
---|
725 | else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab |
---|
726 | G4double pv=mi; |
---|
727 | for(G4int j=0; j<=lastK; j++) // Calculate LogTab values |
---|
728 | { |
---|
729 | lastX[j]=CalcElTot(pv,ind); |
---|
730 | #ifdef pdebug |
---|
731 | G4cout<<"G4QuasiFreeR::FetchElTot:I,j="<<j<<",pv="<<pv<<",E="<<lastX[j].first<<",T=" |
---|
732 | <<lastX[j].second<<G4endl; |
---|
733 | #endif |
---|
734 | if(j!=lastK) pv*=edl; |
---|
735 | } |
---|
736 | i++; // Make a new record to AMDB and position on it |
---|
737 | vI.push_back(lastI); |
---|
738 | vM.push_back(lastM); |
---|
739 | vK.push_back(lastK); |
---|
740 | vX.push_back(lastX); |
---|
741 | } |
---|
742 | else // The A value was found in AMDB |
---|
743 | { |
---|
744 | lastI=vI[i]; |
---|
745 | lastM=vM[i]; |
---|
746 | lastK=vK[i]; |
---|
747 | lastX=vX[i]; |
---|
748 | G4int nextK=lastK+1; |
---|
749 | G4double lpM=lastM+lpi; |
---|
750 | #ifdef pdebug |
---|
751 | G4cout<<"G4QuasiFreeR::FetchElTo:M="<<lpM<<",l="<<lp<<",K="<<lastK<<",n="<<nlp<<G4endl; |
---|
752 | #endif |
---|
753 | if(lp>lpM && lastK<nlp) // LogTab must be updated |
---|
754 | { |
---|
755 | lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab |
---|
756 | #ifdef pdebug |
---|
757 | G4cout<<"G4QuasiFreeR::FetET:K="<<lastK<<",lp="<<lp<<",li="<<lpi<<",dl="<<dl<<G4endl; |
---|
758 | #endif |
---|
759 | if(lastK>nlp) |
---|
760 | { |
---|
761 | lastK=nlp; |
---|
762 | lastM=lpa-lpi; |
---|
763 | } |
---|
764 | else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab |
---|
765 | G4double pv=std::exp(lpM); // momentum of the last calculated beam |
---|
766 | for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values |
---|
767 | { |
---|
768 | pv*=edl; |
---|
769 | lastX[j]=CalcElTot(pv,ind); |
---|
770 | #ifdef pdebug |
---|
771 | G4cout<<"G4QuasiFreeR::FetchElTot:U:j="<<j<<",p="<<pv<<",E="<<lastX[j].first<<",T=" |
---|
772 | <<lastX[j].second<<G4endl; |
---|
773 | #endif |
---|
774 | } |
---|
775 | } // End of LogTab update |
---|
776 | if(lastK>=nextK) // The AMDB was apdated |
---|
777 | { |
---|
778 | vM[i]=lastM; |
---|
779 | vK[i]=lastK; |
---|
780 | } |
---|
781 | } |
---|
782 | // Now one can use tabeles to calculate the value |
---|
783 | G4double dlp=lp-lpi; // Shifted log(p) value |
---|
784 | G4int n=static_cast<int>(dlp/dl); // Low edge number of the bin |
---|
785 | G4double d=dlp-n*dl; // Log shift |
---|
786 | G4double e=lastX[n].first; // E-Base |
---|
787 | lastR.first=e+d*(lastX[n+1].first-e)/dl; // E-Result |
---|
788 | if(lastR.first<0.) lastR.first = 0.; |
---|
789 | G4double t=lastX[n].second; // T-Base |
---|
790 | lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result |
---|
791 | if(lastR.second<0.) lastR.second= 0.; |
---|
792 | #ifdef pdebug |
---|
793 | G4cout<<"=O=>G4QuasiFreeR::FetchElTot:1st="<<lastR.first<<", 2nd="<<lastR.second<<G4endl; |
---|
794 | #endif |
---|
795 | if(lastR.first>lastR.second) lastR.first = lastR.second; |
---|
796 | return lastR; |
---|
797 | } // End of FetchElTot |
---|
798 | |
---|
799 | // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c] |
---|
800 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTot(G4double pIU, G4int hPDG, |
---|
801 | G4int Z, G4int N) |
---|
802 | { |
---|
803 | G4double pGeV=pIU/gigaelectronvolt; |
---|
804 | #ifdef pdebug |
---|
805 | G4cout<<"-->G4QuasiFreeR::GetElTot: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl; |
---|
806 | #endif |
---|
807 | if(Z<1 && N<1) |
---|
808 | { |
---|
809 | G4cout<<"-Warning-G4QuasiFreeRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl; |
---|
810 | return std::make_pair(0.,0.); |
---|
811 | } |
---|
812 | std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); |
---|
813 | std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); |
---|
814 | #ifdef pdebug |
---|
815 | G4cout<<"-OUT->G4QFRat::GetElTot: hp("<<hp.first<<","<<hp.second<<"), hn("<<hn.first<<"," |
---|
816 | <<hn.second<<")"<<G4endl; |
---|
817 | #endif |
---|
818 | G4double A=(Z+N)/millibarn; // To make the result in independent units(IU) |
---|
819 | return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A); |
---|
820 | } // End of GetElTot |
---|
821 | |
---|
822 | // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c] |
---|
823 | std::pair<G4double,G4double> G4QuasiFreeRatios::GetChExFactor(G4double pIU, G4int hPDG, |
---|
824 | G4int Z, G4int N) |
---|
825 | { |
---|
826 | G4double pGeV=pIU/gigaelectronvolt; |
---|
827 | G4double resP=0.; |
---|
828 | G4double resN=0.; |
---|
829 | if(Z<1 && N<1) |
---|
830 | { |
---|
831 | G4cout<<"-Warning-G4QuasiFreeRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl; |
---|
832 | return std::make_pair(resP,resN); |
---|
833 | } |
---|
834 | G4double A=Z+N; |
---|
835 | G4double pf=0.; // Possibility to interact with a proton |
---|
836 | G4double nf=0.; // Possibility to interact with a neutron |
---|
837 | if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N); |
---|
838 | else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z); |
---|
839 | else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310) |
---|
840 | { |
---|
841 | G4double dA=A+A; |
---|
842 | pf=Z/(dA+N+N); |
---|
843 | nf=N/(dA+Z+Z); |
---|
844 | } |
---|
845 | G4double mult=1.; // Factor of increasing multiplicity ( ? @@) |
---|
846 | if(pGeV>.5) |
---|
847 | { |
---|
848 | mult=1./(1.+std::log(pGeV+pGeV))/pGeV; |
---|
849 | if(mult>1.) mult=1.; |
---|
850 | } |
---|
851 | if(pf) |
---|
852 | { |
---|
853 | std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); |
---|
854 | resP=pf*(hp.second/hp.first-1.)*mult; |
---|
855 | } |
---|
856 | if(nf) |
---|
857 | { |
---|
858 | std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); |
---|
859 | resN=nf*(hn.second/hn.first-1.)*mult; |
---|
860 | } |
---|
861 | return std::make_pair(resP,resN); |
---|
862 | } // End of GetChExFactor |
---|
863 | |
---|
864 | // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M) |
---|
865 | // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened |
---|
866 | std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::Scatter(G4int NPDG, |
---|
867 | G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) |
---|
868 | { |
---|
869 | static const G4double mNeut= G4QPDGCode(2112).GetMass(); |
---|
870 | static const G4double mProt= G4QPDGCode(2212).GetMass(); |
---|
871 | static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron |
---|
872 | static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium |
---|
873 | static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3 |
---|
874 | static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha |
---|
875 | G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M) |
---|
876 | N4M/=megaelectronvolt; |
---|
877 | G4LorentzVector tot4M=N4M+p4M; |
---|
878 | #ifdef ppdebug |
---|
879 | G4cerr<<"->G4QFR::Scat:p4M="<<pr4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl; |
---|
880 | #endif |
---|
881 | G4double mT=mNeut; |
---|
882 | G4int Z=0; |
---|
883 | G4int N=1; |
---|
884 | if(NPDG==2212||NPDG==90001000) |
---|
885 | { |
---|
886 | mT=mProt; |
---|
887 | Z=1; |
---|
888 | N=0; |
---|
889 | } |
---|
890 | else if(NPDG==90001001) |
---|
891 | { |
---|
892 | mT=mDeut; |
---|
893 | Z=1; |
---|
894 | N=1; |
---|
895 | } |
---|
896 | else if(NPDG==90002001) |
---|
897 | { |
---|
898 | mT=mHel3; |
---|
899 | Z=2; |
---|
900 | N=1; |
---|
901 | } |
---|
902 | else if(NPDG==90001002) |
---|
903 | { |
---|
904 | mT=mTrit; |
---|
905 | Z=1; |
---|
906 | N=2; |
---|
907 | } |
---|
908 | else if(NPDG==90002002) |
---|
909 | { |
---|
910 | mT=mAlph; |
---|
911 | Z=2; |
---|
912 | N=2; |
---|
913 | } |
---|
914 | else if(NPDG!=2112&&NPDG!=90000001) |
---|
915 | { |
---|
916 | G4cout<<"Error:G4QuasiFreeRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; |
---|
917 | G4Exception("G4QuasiFreeRatios::Scatter:","21",FatalException,"CHIPScomplain"); |
---|
918 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
919 | } |
---|
920 | G4double mT2=mT*mT; |
---|
921 | G4double mP2=pr4M.m2(); |
---|
922 | G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); |
---|
923 | #ifdef pdebug |
---|
924 | G4cerr<<"G4QFR::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl; |
---|
925 | #endif |
---|
926 | G4double E2=E*E; |
---|
927 | if(E<0. || E2<mP2) |
---|
928 | { |
---|
929 | #ifdef ppdebug |
---|
930 | G4cerr<<"-Warning-G4QFR::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl; |
---|
931 | #endif |
---|
932 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
933 | } |
---|
934 | G4double P=std::sqrt(E2-mP2); // Momentum in pseudo laboratory system |
---|
935 | G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer(); |
---|
936 | G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer(); |
---|
937 | #ifdef ppdebug |
---|
938 | G4cout<<"G4QFR::Scatter: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl; |
---|
939 | #endif |
---|
940 | // @@ Temporary NN t-dependence for all hadrons |
---|
941 | if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QElast::Scatter: pPDG="<<pPDG<<G4endl; |
---|
942 | G4int PDG=2212; // *TMP* instead of pPDG |
---|
943 | if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG |
---|
944 | if(!Z && N==1) // Change for Quasi-Elastic on neutron |
---|
945 | { |
---|
946 | Z=1; |
---|
947 | N=0; |
---|
948 | if (PDG==2212) PDG=2112; |
---|
949 | else if(PDG==2112) PDG=2212; |
---|
950 | } |
---|
951 | G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP* |
---|
952 | if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP* |
---|
953 | else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP* |
---|
954 | #ifdef ppdebug |
---|
955 | G4cout<<"G4QElast::Scatter:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl; |
---|
956 | #endif |
---|
957 | #ifdef nandebug |
---|
958 | if(xSec>0. || xSec<0. || xSec==0); |
---|
959 | else G4cout<<"*Warning*G4QElast::Scatter: xSec="<<xSec/millibarn<<G4endl; |
---|
960 | #endif |
---|
961 | // @@ check a possibility to separate p, n, or alpha (!) |
---|
962 | if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing |
---|
963 | { |
---|
964 | #ifdef ppdebug |
---|
965 | G4cerr<<"-Warning-G4QFR::Scat:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl; |
---|
966 | #endif |
---|
967 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action |
---|
968 | } |
---|
969 | G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP* |
---|
970 | if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP* |
---|
971 | else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP* |
---|
972 | G4double maxt=0.; // Prototype of max possible -t |
---|
973 | if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t |
---|
974 | else maxt=NCSmanager->GetHMaxT(); // max possible -t |
---|
975 | #ifdef ppdebug |
---|
976 | G4cout<<"G4QFR::Scat:PDG="<<PDG<<",P="<<P<<",X="<<xSec<<",-t="<<mint<<"<"<<maxt<<", Z=" |
---|
977 | <<Z<<",N="<<N<<G4endl; |
---|
978 | #endif |
---|
979 | #ifdef nandebug |
---|
980 | if(mint>-.0000001); |
---|
981 | else G4cout<<"*Warning*G4QFR::Scat: -t="<<mint<<G4endl; |
---|
982 | #endif |
---|
983 | G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS |
---|
984 | #ifdef ppdebug |
---|
985 | G4cout<<"G4QFR::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl; |
---|
986 | #endif |
---|
987 | if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) |
---|
988 | { |
---|
989 | if (cost>1.) cost=1.; |
---|
990 | else if(cost<-1.) cost=-1.; |
---|
991 | else |
---|
992 | { |
---|
993 | G4double tm=0.; |
---|
994 | if(PDG==2212) tm=PCSmanager->GetHMaxT(); |
---|
995 | else tm=NCSmanager->GetHMaxT(); |
---|
996 | G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl; |
---|
997 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
998 | } |
---|
999 | } |
---|
1000 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon |
---|
1001 | G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); |
---|
1002 | if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost)) |
---|
1003 | { |
---|
1004 | G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl; |
---|
1005 | //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp"); |
---|
1006 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
1007 | } |
---|
1008 | #ifdef ppdebug |
---|
1009 | G4cout<<"G4QFR::Scat:p4M="<<pr4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl; |
---|
1010 | #endif |
---|
1011 | return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result |
---|
1012 | } // End of Scatter |
---|
1013 | |
---|
1014 | // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M) |
---|
1015 | // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened |
---|
1016 | // User should himself change the charge (PDG) (e.g. pn->np, pi+n->pi0p, pi-p->pi0n etc.) |
---|
1017 | std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::ChExer(G4int NPDG, |
---|
1018 | G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) |
---|
1019 | { |
---|
1020 | static const G4double mNeut= G4QPDGCode(2112).GetMass(); |
---|
1021 | static const G4double mProt= G4QPDGCode(2212).GetMass(); |
---|
1022 | G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV(keep p4M) |
---|
1023 | N4M/=megaelectronvolt; |
---|
1024 | G4LorentzVector tot4M=N4M+p4M; |
---|
1025 | G4int Z=0; |
---|
1026 | G4int N=1; |
---|
1027 | G4int RPDG=2212; // PDG of the recoil nucleon |
---|
1028 | G4int sPDG=0; // PDG code of the scattered hadron |
---|
1029 | G4double mS=0.; // proto of mass of scattered hadron |
---|
1030 | G4double mT=mProt; // mass of the recoil nucleon |
---|
1031 | if(NPDG==2212) |
---|
1032 | { |
---|
1033 | mT=mNeut; |
---|
1034 | Z=1; |
---|
1035 | N=0; |
---|
1036 | RPDG=2112; // PDG of the recoil nucleon |
---|
1037 | if(pPDG==-211) sPDG=111; // pi+ -> pi0 |
---|
1038 | else if(pPDG==-321) |
---|
1039 | { |
---|
1040 | sPDG=310; // K+ -> K0S |
---|
1041 | if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L |
---|
1042 | } |
---|
1043 | else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321; // K0 -> K+ (?) |
---|
1044 | else if(pPDG==3112) sPDG=3212; // Sigma- -> Sigma0 |
---|
1045 | else if(pPDG==3212) sPDG=3222; // Sigma0 -> Sigma+ |
---|
1046 | else if(pPDG==3312) sPDG=3322; // Xi- -> Xi0 |
---|
1047 | } |
---|
1048 | else if(NPDG==2112) // Default |
---|
1049 | { |
---|
1050 | if(pPDG==211) sPDG=111; // pi+ -> pi0 |
---|
1051 | else if(pPDG==321) |
---|
1052 | { |
---|
1053 | sPDG=310; // K+ -> K0S |
---|
1054 | if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L |
---|
1055 | } |
---|
1056 | else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0 -> K- (?) |
---|
1057 | else if(pPDG==3222) sPDG=3212; // Sigma+ -> Sigma0 |
---|
1058 | else if(pPDG==3212) sPDG=3112; // Sigma0 -> Sigma- |
---|
1059 | else if(pPDG==3322) sPDG=3312; // Xi0 -> Xi- |
---|
1060 | } |
---|
1061 | else |
---|
1062 | { |
---|
1063 | G4cout<<"Error:G4QuasiFreeRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; |
---|
1064 | G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain"); |
---|
1065 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
1066 | } |
---|
1067 | if(sPDG) mS=mNeut; |
---|
1068 | else |
---|
1069 | { |
---|
1070 | G4cout<<"Error:G4QuasiFreeRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl; |
---|
1071 | G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain"); |
---|
1072 | //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception |
---|
1073 | } |
---|
1074 | G4double mT2=mT*mT; |
---|
1075 | G4double mS2=mS*mS; |
---|
1076 | G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT); |
---|
1077 | G4double E2=E*E; |
---|
1078 | if(E<0. || E2<mS2) |
---|
1079 | { |
---|
1080 | #ifdef pdebug |
---|
1081 | G4cerr<<"-Warning-G4QFR::ChEx:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mS2<<G4endl; |
---|
1082 | #endif |
---|
1083 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
1084 | } |
---|
1085 | G4double P=std::sqrt(E2-mS2); // Momentum in pseudo laboratory system |
---|
1086 | G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer(); |
---|
1087 | G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer(); |
---|
1088 | #ifdef debug |
---|
1089 | G4cout<<"G4QFR::ChExer: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl; |
---|
1090 | #endif |
---|
1091 | // @@ Temporary NN t-dependence for all hadrons |
---|
1092 | G4int PDG=2212; // *TMP* instead of pPDG |
---|
1093 | if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG |
---|
1094 | if(!Z && N==1) // Change for Quasi-Elastic on neutron |
---|
1095 | { |
---|
1096 | Z=1; |
---|
1097 | N=0; |
---|
1098 | if (PDG==2212) PDG=2112; |
---|
1099 | else if(PDG==2112) PDG=2212; |
---|
1100 | } |
---|
1101 | G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP* |
---|
1102 | if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP* |
---|
1103 | else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP* |
---|
1104 | #ifdef debug |
---|
1105 | G4cout<<"G4QElast::ChExer:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl; |
---|
1106 | #endif |
---|
1107 | #ifdef nandebug |
---|
1108 | if(xSec>0. || xSec<0. || xSec==0); |
---|
1109 | else G4cout<<"*Warning*G4QElast::ChExer: xSec="<<xSec/millibarn<<G4endl; |
---|
1110 | #endif |
---|
1111 | // @@ check a possibility to separate p, n, or alpha (!) |
---|
1112 | if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing |
---|
1113 | { |
---|
1114 | #ifdef pdebug |
---|
1115 | G4cerr<<"-Warning-G4QFR::ChEx:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl; |
---|
1116 | #endif |
---|
1117 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action |
---|
1118 | } |
---|
1119 | G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP* |
---|
1120 | if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP* |
---|
1121 | else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP* |
---|
1122 | #ifdef pdebug |
---|
1123 | G4cout<<"G4QFR::ChEx:PDG="<<pPDG<<", P="<<P<<", CS="<<xSec<<", -t="<<mint<<G4endl; |
---|
1124 | #endif |
---|
1125 | #ifdef nandebug |
---|
1126 | if(mint>-.0000001); |
---|
1127 | else G4cout<<"*Warning*G4QFR::ChExer: -t="<<mint<<G4endl; |
---|
1128 | #endif |
---|
1129 | G4double maxt=0.; // Prototype of max possible -t |
---|
1130 | if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t |
---|
1131 | else maxt=NCSmanager->GetHMaxT(); // max possible -t |
---|
1132 | G4double cost=1.-mint/maxt; // cos(theta) in CMS |
---|
1133 | #ifdef pdebug |
---|
1134 | G4cout<<"G4QuasiFfreeRatio::ChExer: -t="<<mint<<", maxt="<<maxt<<", cost="<<cost<<G4endl; |
---|
1135 | #endif |
---|
1136 | if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) |
---|
1137 | { |
---|
1138 | if (cost>1.) cost=1.; |
---|
1139 | else if(cost<-1.) cost=-1.; |
---|
1140 | else |
---|
1141 | { |
---|
1142 | G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl; |
---|
1143 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
1144 | } |
---|
1145 | } |
---|
1146 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon |
---|
1147 | pr4M=G4LorentzVector(0.,0.,0.,mS); // 4mom of the scattered hadron |
---|
1148 | G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); |
---|
1149 | if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost)) |
---|
1150 | { |
---|
1151 | G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl; |
---|
1152 | //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp"); |
---|
1153 | return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action |
---|
1154 | } |
---|
1155 | #ifdef debug |
---|
1156 | G4cout<<"G4QFR::ChEx:p4M="<<p4M<<"+r4M="<<reco4M<<"="<<p4M+reco4M<<"="<<tot4M<<G4endl; |
---|
1157 | #endif |
---|
1158 | return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result |
---|
1159 | } // End of ChExer |
---|
1160 | |
---|
1161 | // Calculate ChEx/El ratio (p is in independent units, (Z,N) is target, pPDG is projectile) |
---|
1162 | G4double G4QuasiFreeRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG) |
---|
1163 | { |
---|
1164 | p/=MeV; // Converted from independent units |
---|
1165 | G4double A=Z+N; |
---|
1166 | if(A<1.5) return 0.; |
---|
1167 | G4double C=0.; |
---|
1168 | if (pPDG==2212) C=N/(A+Z); |
---|
1169 | else if(pPDG==2112) C=Z/(A+N); |
---|
1170 | else G4cout<<"*Warning*G4QCohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl; |
---|
1171 | C*=C; // Coherent processes squares the amplitude |
---|
1172 | // @@ This is true only for nucleons: other projectiles must be treated differently |
---|
1173 | G4double sp=std::sqrt(p); |
---|
1174 | G4double p2=p*p; |
---|
1175 | G4double p4=p2*p2; |
---|
1176 | G4double dl1=std::log(p)-5.; |
---|
1177 | G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013); |
---|
1178 | G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p; |
---|
1179 | G4double R=U/T; |
---|
1180 | return C*R*R; |
---|
1181 | } |
---|