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 | // G4 Tools program: Glauber Hadron-Nuclear Cross Sections |
---|
28 | // ..................................................... |
---|
29 | // Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Jan-2005 |
---|
30 | // |
---|
31 | //===================================================================== |
---|
32 | |
---|
33 | //#define debug |
---|
34 | //#define bestest |
---|
35 | //#define integrcs |
---|
36 | #define eldiffcs |
---|
37 | |
---|
38 | #include "globals.hh" |
---|
39 | #include <iostream> |
---|
40 | #include <fstream> |
---|
41 | #include <vector> |
---|
42 | #include "G4ios.hh" |
---|
43 | #include "G4QBesIKJY.hh" |
---|
44 | #include "G4QPDGCode.hh" |
---|
45 | #include "G4NucleiPropertiesTable.hh" |
---|
46 | #include "G4NucleiProperties.hh" |
---|
47 | //#include <CLHEP/GenericFunctions/Bessel.hh> |
---|
48 | |
---|
49 | // Integrated cross-sections |
---|
50 | G4double TotalCrSec, InelCrSec, ProdCrSec, ElasticCrSec, QuasyElasticCrSec; |
---|
51 | // Parameters of hadron-nucleon interaction |
---|
52 | G4double HadrTot, HadrSlope, HadrReIm, DDSect2, DDSect3, Tot00; |
---|
53 | // CHIPS parameters of hadron-proton & hadron-neutron interaction (Re/Im is common) |
---|
54 | G4double HadrTElN, HadrTElP, HadrSlpN, HadrSlpP, HadrPexN, HadrPexP; |
---|
55 | |
---|
56 | // CHIPS Parameters (only for protons, nn is necessary for neutrons np+nn) |
---|
57 | void CHIPSParameters(G4int iPDG, G4double HadrMoment) // HadrMoment is in MeV |
---|
58 | { |
---|
59 | G4double p=HadrMoment/GeV; // Momentum in GeV |
---|
60 | G4double p2=p*p; // Squared momentum in GeV^2 |
---|
61 | G4double sp=std::sqrt(p); |
---|
62 | G4double p2s=p2*sp; |
---|
63 | G4double p3=p2*p; |
---|
64 | G4double p4=p2*p2; |
---|
65 | G4double p5=p4*p; |
---|
66 | G4double lp=std::log(p); |
---|
67 | G4double d2=(lp-5.)*(lp-5.); |
---|
68 | G4double p8=p4*p4; |
---|
69 | if (iPDG==2212) |
---|
70 | { |
---|
71 | HadrTElN=12./(p2s+.05*p+.0001/sp)+.35/p+(6.75+.14*d2+19./p)/(1.+.6/p4); // pn (mb) |
---|
72 | HadrTElP=2.91/(p2s+.0024)+(6.75+.14*d2+23./p)/(1.+.4/p4); // pp (mb) |
---|
73 | HadrSlpN=(7.2+4.32/(p8+.012*p3))/(1.+2.5/p4); // pn (GeV^-2) |
---|
74 | HadrSlpP=8.*std::pow(p,.055)/(1.+3.64/p4); // pp (GeV^-2) |
---|
75 | HadrPexN=(6.75+.14*d2+13./p)/(1.+.14/p4)+.6/(p4+.00013); // pn (mb/sr) |
---|
76 | HadrPexP=(74.+3.*d2)/(1.+3.4/p5)+(.2/p2+17.*p)/(p4+.001*sp); // pn (mb/sr) |
---|
77 | HadrReIm=-.54+.0954*lp+.975/(p2+.09865/p); // no unit (the same for pp & pn) |
---|
78 | //G4cout<<"CHIPSPar:p="<<p<<",N="<<HadrTotN<<",P="<<HadrTotN<<",RI="<<HadrReIm<<G4endl; |
---|
79 | } |
---|
80 | else G4cout<<"CHIPSPar: Only proton projectile (2212) is implemented PDG="<<iPDG<<G4endl; |
---|
81 | } |
---|
82 | |
---|
83 | // Dubna Parameters |
---|
84 | G4int CalculateParameters(G4int iPDG, G4double HadrMoment) // HadrMoment is in MeV |
---|
85 | { |
---|
86 | static const G4double mN=.938; // Dubna's mass of the nucleon in nuclei in GeV |
---|
87 | //static const G4double mN=.9315; // AtomicMassUnit (mass of the nucleon in nuclei) in GeV |
---|
88 | static const G4double dmN=mN+mN; // Doubled mass of the nucleon in nuclei in GeV |
---|
89 | static const G4double mN2=mN*mN; // Squared mass of the nucleon in nuclei in GeV^2 |
---|
90 | G4int iHadron=-1; |
---|
91 | if(iPDG==2212||iPDG==2112||iPDG==3122||iPDG==3222||iPDG==3112|| |
---|
92 | iPDG==3212||iPDG==3312||iPDG==3322||iPDG==3334) iHadron = 0; |
---|
93 | else if(iPDG==-2212||iPDG==-2112||iPDG==-3122||iPDG==-3222|| |
---|
94 | iPDG==-3112||iPDG==-3212||iPDG==-3312||iPDG==-3322|| |
---|
95 | iPDG==-3334) iHadron = 1; // Anti-nucleons, Anti-hyperons |
---|
96 | else if(iPDG== 211) iHadron = 2; // Pi+ |
---|
97 | else if(iPDG==-211) iHadron = 3; // Pi- |
---|
98 | else if(iPDG== 321) iHadron = 4; // K+ |
---|
99 | else if(iPDG==-321) iHadron = 5; // K- @@ What about K0? |
---|
100 | else |
---|
101 | { |
---|
102 | G4cout<<"CalculateParameters: PDG= "<<iPDG<<" is not supported"<<G4endl; |
---|
103 | return iHadron; |
---|
104 | } |
---|
105 | G4double mHadr = G4QPDGCode(iPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
106 | G4double mHadr2 = mHadr*mHadr; // Squared Hadron Mass in GeV |
---|
107 | HadrMoment/= 1000.; // Momentum in GeV (input parameter) |
---|
108 | G4double HadrMoment2= HadrMoment*HadrMoment; // Squared Momentum in GeV^2 |
---|
109 | G4double HadrEnergy2= mHadr2+HadrMoment2; // Tot energy in GeV |
---|
110 | G4double HadrEnergy = std::sqrt(HadrEnergy2); // Tot energy in GeV (global parameter) |
---|
111 | G4double sHadr = HadrEnergy*dmN+mN2+mHadr2; // s in GeV^2 |
---|
112 | G4double sqrS = std::sqrt(sHadr); // W in GeV |
---|
113 | //G4cout<<"GHAD: E="<<HadrEnergy<<",dM="<<dmN<<",sN="<<mN2<<",sH="<<mHadr2<<G4endl; |
---|
114 | switch (iHadron) |
---|
115 | { |
---|
116 | case 0: // =====>>> proton or hyperons (the same?) |
---|
117 | //G4double Delta=1; // LowEnergy corr |
---|
118 | //if(HadrEnergy<40.) Delta = 0.916+0.0021*HadrEnergy; |
---|
119 | HadrTot = 5.2+5.2*std::log(HadrEnergy)+51*std::pow(HadrEnergy,-0.35); // mb |
---|
120 | HadrSlope = 6.44+0.88*std::log(sHadr)-1; // GeV^-2 |
---|
121 | HadrReIm = 0.13*std::log(sHadr/350)*std::pow(sHadr,-0.18); // no unit |
---|
122 | //G4cout<<"GHAD:S="<<sHadr<<",T="<<HadrTot<<",R="<<HadrReIm<<",B="<<HadrSlope<<G4endl; |
---|
123 | DDSect2 = 11.; // mb*GeV^-2 |
---|
124 | DDSect3 = 3.; // mb*GeV^-2 |
---|
125 | // Hyperon corrections |
---|
126 | if(iPDG==3122||iPDG==3222||iPDG==3112||iPDG==3212) //Lambda,Sig0+- |
---|
127 | { |
---|
128 | HadrTot *= 0.80; |
---|
129 | HadrSlope *= 0.85; |
---|
130 | } |
---|
131 | else if(iPDG==3312||iPDG==3322) // ---> Xi0- |
---|
132 | { |
---|
133 | HadrTot *= 0.70; |
---|
134 | HadrSlope *= 0.75; |
---|
135 | } |
---|
136 | else if(iPDG==3334) // ---> Omega- |
---|
137 | { |
---|
138 | HadrTot *= 0.60; |
---|
139 | HadrSlope *= 0.65; |
---|
140 | } |
---|
141 | break; |
---|
142 | case 1: // =====>>> antiproton and antihyperons |
---|
143 | HadrTot = 5.2+5.2*std::log(HadrEnergy)+123.2*std::pow(HadrEnergy,-0.5); // mb |
---|
144 | HadrSlope = 8.32+0.57*std::log(sHadr); // GeV^-2 |
---|
145 | if(HadrEnergy<1000.) HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8); |
---|
146 | else HadrReIm = 0.6*std::log(sHadr/350)*std::pow(sHadr,-0.25); |
---|
147 | DDSect2 = 11.; //mb*GeV-2 @@ the same as in case 0 |
---|
148 | DDSect3 = 3.; //mb*GeV-2 @@ the same as in case 0 |
---|
149 | // Hyperon corrections |
---|
150 | if(iPDG==-3122||iPDG==-3222||iPDG==-3112||iPDG==-3212)//AntiLam,Sig |
---|
151 | { |
---|
152 | HadrTot *=0.75; |
---|
153 | HadrSlope *=0.85; |
---|
154 | } |
---|
155 | else if(iPDG==-3312||iPDG==-3322) // =====>>> AntiXi0- |
---|
156 | { |
---|
157 | HadrTot *=0.65; |
---|
158 | HadrSlope *=0.75; |
---|
159 | } |
---|
160 | if(iPDG==-3334) // ======>>>AntiOmega- |
---|
161 | { |
---|
162 | HadrTot *=0.55; |
---|
163 | HadrSlope *=0.65; |
---|
164 | } |
---|
165 | break; |
---|
166 | case 2: // =====>>> pi plus |
---|
167 | if(HadrMoment>2.) HadrTot = 10.6+2.*std::log(HadrEnergy)+25*std::pow(HadrEnergy,-0.43);// mb |
---|
168 | else HadrTot = 40-50*(HadrMoment-1.5)*(HadrMoment-1.7); |
---|
169 | HadrSlope = 7.28+0.245*std::log(sHadr); // GeV^-2 |
---|
170 | HadrReIm = 0.2*std::log(sHadr/100)*std::pow(sHadr,-0.15); // no dim |
---|
171 | DDSect2 = 4.6; // mb*GeV^-2 |
---|
172 | DDSect3 = 1.33; // mb*GeV^-2 |
---|
173 | break; |
---|
174 | case 3: // =====>>> pi minus |
---|
175 | HadrTot = 10.6+2*std::log(HadrEnergy)+30*std::pow(HadrEnergy,-0.43); // mb @@ E->inf? |
---|
176 | if(HadrMoment<1.399) HadrTot = HadrTot+21.0/0.4*(1.4-HadrMoment); // A huge jump ? |
---|
177 | HadrSlope = 7.28+0.245*std::log(sHadr); // GeV-2 |
---|
178 | HadrReIm = 0.2*std::log(sHadr/100)*std::pow(sHadr,-0.15); |
---|
179 | DDSect2 = 4.6; //mb*GeV-2 @@ the same as in case 2 |
---|
180 | DDSect3 = 1.33; //mb*GeV-2 |
---|
181 | break; |
---|
182 | case 4: // =====>>> K plus |
---|
183 | HadrTot = 10.6+1.8*std::log(HadrEnergy)+9.*std::pow(HadrEnergy,-0.55); // mb |
---|
184 | if(HadrEnergy>100) HadrSlope = 15.; // Jump at 100 GeV (?) |
---|
185 | else HadrSlope = 1.0+1.76*std::log(sHadr)-2.84*std::pow(sHadr,-0.5); // GeV-2 |
---|
186 | //else HadrSlope = 5.28+1.76*std::log(sHadr)-2.84*std::pow(sHadr,-0.5); // GeV-2 |
---|
187 | HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1); |
---|
188 | DDSect2 = 3.5; //mb*GeV-2 |
---|
189 | DDSect3 = 1.03; //mb*GeV-2 |
---|
190 | break; |
---|
191 | case 5: // ======>>> K minus |
---|
192 | HadrTot = 10+1.8*std::log(HadrEnergy)+25*std::pow(HadrEnergy,-0.5); // mb |
---|
193 | HadrSlope = 6.98+0.127*std::log(sHadr); // GeV-2 |
---|
194 | //if(HadrEnergy<8) HadrReIm = 0.7; |
---|
195 | HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1); |
---|
196 | DDSect2 = 3.5; //mb*GeV-2 |
---|
197 | DDSect3 = 1.03; //mb*GeV-2 |
---|
198 | break; |
---|
199 | } |
---|
200 | return iHadron; |
---|
201 | } |
---|
202 | |
---|
203 | void CalculateIntegralCS(G4int Anuc, G4double HadrEnergy) |
---|
204 | { |
---|
205 | static const G4double eps = 0.0001; // Accuracy of calculations |
---|
206 | static const G4double mb2G2 = 2.568; // Transformation from mb to GeV^-2 |
---|
207 | static const G4double incrf = 10; // Increase factor of radius for intergation |
---|
208 | static const G4double m2G10 = incrf*mb2G2; // Big conversion factor |
---|
209 | static const G4double dR0 = 1.06+1.06; // Doubled nominal R0 |
---|
210 | G4double An = Anuc; // A float |
---|
211 | G4double Stot = HadrTot*mb2G2; // Total hadron-nucleon CS in GeV-2 |
---|
212 | G4double Bhad = HadrSlope; // Diffraction cone nH In GeV-2 |
---|
213 | G4double Asq = 1+HadrReIm*HadrReIm; // |M|^2/(ImM)^2 |
---|
214 | |
---|
215 | G4double R0 = std::sqrt(0.99); // in fermi (strange value? 1.07?) |
---|
216 | |
---|
217 | if (Anuc == 58) R0 = std::sqrt(0.6); // Personal correction for Fe |
---|
218 | if (Anuc >20) R0 = std::sqrt((35.34+0.5*Anuc)/(40.97+Anuc)); |
---|
219 | if (Anuc == 16) R0 = std::sqrt(0.75); // Personal correction for O16 |
---|
220 | if (Anuc >10) R0 = std::sqrt(0.84); // The same R0 for A=11-20 (except O16) |
---|
221 | if (Anuc == 4) R0 *= 1.2; // Personal correction for Helium |
---|
222 | |
---|
223 | G4double Rnucl = R0*std::pow(static_cast<double>(Anuc),.3333333); // in fermi |
---|
224 | G4double Rnuc2 = Rnucl*Rnucl*m2G10; // in GeV^-2 (10 -> to make it >> R) |
---|
225 | G4double RB = Rnuc2+Bhad; // Effective increase of R2 (Convolution?) |
---|
226 | G4double R2B = RB+Bhad; // Over-convolution (?) |
---|
227 | G4double Delta = Stot/R2B/twopi; // Integration step for Total & Inelastic CS |
---|
228 | G4double Delt2 = Delta+Delta; // Two delta |
---|
229 | G4double Delt3 = Stot*Asq*R2B/RB/Bhad/16/pi; // Integration step for Production CS |
---|
230 | |
---|
231 | G4double Tot0 = 0.; // =SUM[(-Delta)^(i-1)*(C^A_i)/i] |
---|
232 | G4double Inel0 = 0.; // =SUM[(-Delta)^(i-1)*(2A)!/(2A-i)!/i!/i] |
---|
233 | G4double N = -1./Delta; |
---|
234 | G4double N1 = N; |
---|
235 | G4double N3 = -1./Delt2; |
---|
236 | G4double Prod0 = 0.; |
---|
237 | G4int Ai = Anuc+1; // Working value to avoid multiplication |
---|
238 | G4double RBi = 0.; // Working value to avoid multiplication |
---|
239 | for (G4int i=1; i<= Anuc; i++) // ==> Loop over nucleons |
---|
240 | { |
---|
241 | Ai--; // Working value to avoid multiplication |
---|
242 | RBi += RB; // Working value to avoid multiplication |
---|
243 | G4double Di=Delta/i; // Working value to avoid repetition |
---|
244 | N *= -Di*Ai; // (-Delta)^(i-1)*A!/(A-i)!/i! (C^A_i) |
---|
245 | N1 *= -Di*(Anuc+Ai); // (-Delta)^(i-1)*(2A)!/(2A-i)!/i! |
---|
246 | N3 *= -Delt2*Ai/i; // (-2*Delt)^(i-1)*A!/(A-i)!/i! (C^A_i) |
---|
247 | Tot0 += N/i; // =SUM[(-Delta)^(i-1)*(C^A_i)/i] |
---|
248 | Inel0 += N1/i; // =SUM[(-Delta)^(i-1)*(2A)!/(2A-i)!/i!/i] |
---|
249 | //G4double N2 = 1.; |
---|
250 | G4double N4 = 1.; |
---|
251 | //G4double Inel1 = 0.; |
---|
252 | //G4double Inel2 = 0.; |
---|
253 | G4double Prod1 = 0.; |
---|
254 | G4double Bhl = 0.; |
---|
255 | for (G4int l=0; l<= i; l++) // ==> Loop over nucleon-nucleon pairs (0?) |
---|
256 | { |
---|
257 | //Inel1 += N2/(i+l); // SUM[(-Delta)^(l-1)*i!(i-1)!/(l-1)!/l!/(i+l)!] |
---|
258 | Prod1 += N4*RB/(RBi+Bhl); // |
---|
259 | //N2 *= -Delt*(i-l)/(l+1); // (-Delta)^l*i!/l!/(l+1)! |
---|
260 | N4 *= -Delt3*(i-l)/(l+1); |
---|
261 | Bhl += Bhad; |
---|
262 | } |
---|
263 | //Inel2 += Inel1*N3; |
---|
264 | Prod0 += Prod1*N3; |
---|
265 | if(std::fabs(N1/i/Inel0) < eps) break; |
---|
266 | } |
---|
267 | |
---|
268 | Tot0 *= HadrTot; // Multiply by the hN cross section |
---|
269 | Inel0 *= HadrTot/2; // Multiply by a half of hN cross section |
---|
270 | //Inel2 *= HadrTot; // Multiply by the hN cross section |
---|
271 | Prod0 *= HadrTot; // Multiply by the hN cross section |
---|
272 | Tot00 = Tot0; // Remember because it will be changed |
---|
273 | G4double ak = Rnuc2*twopi/Stot; // cross-section ratio for integration |
---|
274 | G4double Ank = An/ak; // Working value to avoid repetition |
---|
275 | G4double rak = 1./ak; // Working value to avoid repetition |
---|
276 | G4double ak1 = 1.-rak/4; // Working value to avoid repetition |
---|
277 | G4double Ank1 = Ank*ak1; // Working value to avoid repetition |
---|
278 | G4double DDSect1 = (DDSect2+DDSect3*std::log(dR0*HadrEnergy/Rnucl/std::sqrt(m2G10)/4)); |
---|
279 | G4double Dtot = 8*pi*ak/HadrTot*(1.-(1.+Ank)*std::exp(-Ank))*DDSect1/mb2G2;//TotCSCor |
---|
280 | G4double bk = (1.-rak)/Stot/ak1; // Working value to avoid power(,2) |
---|
281 | G4double bd = bk*bk*DDSect1*(1.-(1.+Ank1)*std::exp(-Ank1))*Rnuc2; // Working |
---|
282 | G4double Dprod = bd*4*pi2*mb2G2; // Correction for the Production & Inel CS |
---|
283 | TotalCrSec = Tot0-Dtot; // Total Cross Section (primary) |
---|
284 | InelCrSec = Inel0-Dprod; // Inelastic Cross Section (primary) |
---|
285 | //InelCrSec1 = Inel2; // Direct calculation without correction |
---|
286 | ProdCrSec = Prod0-Dprod; // Production Cross Section (primary) |
---|
287 | ElasticCrSec = TotalCrSec-InelCrSec; // Elastic Cross Section (secondary) |
---|
288 | QuasyElasticCrSec = InelCrSec-ProdCrSec; // Quasy-Elastic Cross Section (secondary) |
---|
289 | } |
---|
290 | |
---|
291 | #ifdef eldiffcs |
---|
292 | // Functions for calculation of differential elastic cross-sections |
---|
293 | // Static globals |
---|
294 | static const G4double rAmax = 2.5; // Factor of maximum radius |
---|
295 | static const G4int NptB = 500; // A#of steps in the impact parameter integration |
---|
296 | static const G4int nFact = 253; // Maximum A (for Factorials and Factors) |
---|
297 | |
---|
298 | // Globals |
---|
299 | G4double r0 = 1.1; // Wood-Saxon's factor |
---|
300 | G4double Mnoj[nFact], Factorials[nFact]; // Factors and factorials |
---|
301 | G4double R1, R2, Pnucl, Aeff; // Nuclear Parameters |
---|
302 | G4double InCoh; // Incoherent quasi-elastic cross section |
---|
303 | |
---|
304 | G4double binom(G4int N, G4int M) |
---|
305 | { |
---|
306 | if(N>170 && N>M && M) |
---|
307 | { |
---|
308 | G4double fN=N; |
---|
309 | G4double fM=M; |
---|
310 | G4double fD=N-M; |
---|
311 | G4double man=N*std::log(fN)-M*std::log(fM)-fD*std::log(fD)-1.; |
---|
312 | //G4cout<<"G4GCS::binom: N="<<N<<", M="<<M<<", C="<<man<<", E="<<std::exp(man)<<G4endl; |
---|
313 | return std::exp(man); |
---|
314 | } |
---|
315 | else if(N>1 && N>M && M) return Factorials[N]/Factorials[M]/Factorials[N-M]; |
---|
316 | return 1.; |
---|
317 | } |
---|
318 | |
---|
319 | // Initialize factorials and combinatoric coefficients |
---|
320 | void Initialize() |
---|
321 | { |
---|
322 | //static const G4double sat = 3.03; // Stirling saturation of Dubna |
---|
323 | G4double Val = Factorials[0] = Mnoj[0] = 1.; |
---|
324 | for (G4int i=1; i<nFact; i++) |
---|
325 | { |
---|
326 | G4double Si = 0.; |
---|
327 | Val *= i; // Calculation of the factorial values |
---|
328 | Factorials[i] = Val; |
---|
329 | for (G4int l = 0; l<=i; l++) |
---|
330 | { |
---|
331 | G4double Sm = 1.; // sum of reversed binom coefficients (?) |
---|
332 | for (G4int m=1; m<=i-l; m++) Sm += 1./binom(i-l,m); |
---|
333 | //if(i>169)G4cout<<"G4GCS::Init:S="<<Sm<<", b("<<i<<","<<l<<")="<<binom(i,l)<<G4endl; |
---|
334 | Si += Sm/binom(i,l); |
---|
335 | } // End of l LOOP |
---|
336 | Mnoj[i] = Si; |
---|
337 | //G4double d=Si-sat; |
---|
338 | //G4double r=d/sat; |
---|
339 | //if(i>100)G4cout<<"G4GCS::Init:i="<<i<<":"<<Si<<"-"<<sat<<" = "<<d<<", r="<<r<<G4endl; |
---|
340 | } // End of i LOOP |
---|
341 | } // End of Initialization |
---|
342 | |
---|
343 | // Dubna paranmeters of nuclei |
---|
344 | void CalculateNuclearParameters(G4int Anuc) // R1,R2,Pnuc,Aeff |
---|
345 | {// ============================================== |
---|
346 | G4double A=Anuc; |
---|
347 | // Personal corrections fore some nuclei |
---|
348 | if(Anuc == 4) // Personal correction for He (!) |
---|
349 | { |
---|
350 | R1 = 5.5; |
---|
351 | R2 = 3.7; |
---|
352 | Pnucl = 0.4; |
---|
353 | Aeff = 0.87; |
---|
354 | } |
---|
355 | else if(Anuc == 9) // Personal correction for Be (!) |
---|
356 | { |
---|
357 | R1 = 9.0; |
---|
358 | R2 = 7.0; |
---|
359 | Pnucl = 0.190; |
---|
360 | Aeff = 0.9; |
---|
361 | } |
---|
362 | //else if(Anuc == 11) // Personal correction for B (!) |
---|
363 | //{ |
---|
364 | // R1 = 10.8; |
---|
365 | // R2 = 7.5; |
---|
366 | // Pnucl = 0.85; |
---|
367 | // Aeff = 1.2; |
---|
368 | //} |
---|
369 | //else if(Anuc == 12) // Personal correction for C (!) |
---|
370 | //{ |
---|
371 | // R1 = 9.336; |
---|
372 | // R2 = 5.63; |
---|
373 | // Pnucl = 0.197; |
---|
374 | // Aeff = 1.0; |
---|
375 | //} |
---|
376 | else if(Anuc == 16) // Personal correction for O (!) |
---|
377 | { |
---|
378 | R1 = 10.50; |
---|
379 | R2 = 5.5; |
---|
380 | Pnucl = 0.7; |
---|
381 | Aeff = 0.98; |
---|
382 | //R1 = 11.3; |
---|
383 | //R2 = 2.5; |
---|
384 | //Pnucl = 0.75; |
---|
385 | //Aeff = 0.9; |
---|
386 | } |
---|
387 | else if(Anuc == 58) // Personal correction for Ni (!) |
---|
388 | { |
---|
389 | R1 = 15.0; |
---|
390 | R2 = 9.9; |
---|
391 | Pnucl = 0.45; |
---|
392 | Aeff = 0.85; |
---|
393 | } |
---|
394 | else if(Anuc == 90) // Personal correction for Zr (!) |
---|
395 | { |
---|
396 | //R1 = 16.5; |
---|
397 | //R2 = 11.62; |
---|
398 | //Pnucl = 0.4; |
---|
399 | //Aeff = 0.9; |
---|
400 | R1 = 16.5; |
---|
401 | R2 = 11.62; |
---|
402 | Pnucl = 0.4; |
---|
403 | Aeff = 0.7; |
---|
404 | } |
---|
405 | else if(Anuc == 208) // Personal correction for Pb (!) |
---|
406 | { |
---|
407 | // R1 = 20.73; R2 = 15.74. |
---|
408 | //R1 = 4.1408*std::pow(A,0.3018); |
---|
409 | //R2 = 3.806*std::pow(A-10.068,0.2685); |
---|
410 | //Pnucl = 0.9; |
---|
411 | //Aeff = 1.1; |
---|
412 | R1 = 19.5; |
---|
413 | R2 = 15.74; |
---|
414 | Pnucl = 0.4; |
---|
415 | Aeff = 0.7; |
---|
416 | } |
---|
417 | else |
---|
418 | { |
---|
419 | R1 = 4.45*std::pow(A-1.,.309); // First diffractional radius |
---|
420 | if(Anuc == 28) R1 *= .955; // Personal correction for Si (?!) |
---|
421 | R2 = 2.3*std::pow(A,.36); // Second diffraction radius |
---|
422 | Pnucl = 0.176+(.00167+.00000869*Anuc)*Anuc; // Momentum of mean Fermi motion |
---|
423 | Aeff = 0.9; // Effective screaning |
---|
424 | } |
---|
425 | } |
---|
426 | |
---|
427 | // Independent transition method calculation. Has problems for A>95 (short cut). |
---|
428 | G4double DiffElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) // All in MeV |
---|
429 | {// ================================================================= |
---|
430 | static const G4double eps = 0.000001; // Accuracy of calculations |
---|
431 | static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
432 | static const G4double piMG = pi/mb2G2; // PiTrans from GeV^-2 to mb |
---|
433 | static const G4double dR0 = 1.06+1.06; // Doubled nominal R0 |
---|
434 | G4double hMom = HadrMom/1000.; // Momentum (GeV, inputParam) |
---|
435 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
436 | G4double MassH2 = MassH*MassH; |
---|
437 | G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2); // Tot energy in GeV |
---|
438 | //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl; |
---|
439 | if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3"); |
---|
440 | if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!"); |
---|
441 | //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)"); |
---|
442 | CalculateParameters(hPDG, HadrMom); |
---|
443 | CalculateIntegralCS(A, HadrEnergy); |
---|
444 | CalculateNuclearParameters(A); |
---|
445 | G4double Q2 = aQ2/1000/1000; // in GeV^2 |
---|
446 | //MassN = A*0.938; // @@ This is bad! |
---|
447 | G4double MassN = A*0.9315; // @@ Atomic Unit is bad too |
---|
448 | G4double MassN2 = MassN*MassN; |
---|
449 | G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s |
---|
450 | G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // CM energy of a hadron |
---|
451 | G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum |
---|
452 | |
---|
453 | G4double Stot = HadrTot*mb2G2; // in GeV^-2 |
---|
454 | G4double Bhad = HadrSlope; // in GeV^-2 |
---|
455 | G4double Asq = 1+HadrReIm*HadrReIm; // |M|^2/(ImM)^2 |
---|
456 | G4double Rho2 = std::sqrt(Asq); // M/ImM |
---|
457 | G4double R12 = R1*R1; |
---|
458 | G4double R22 = R2*R2; |
---|
459 | G4double R12B = R12+Bhad+Bhad; |
---|
460 | G4double R22B = R22+Bhad+Bhad; |
---|
461 | G4double R12Ap = R12+20; |
---|
462 | G4double R22Ap = R22+20; |
---|
463 | G4double R13Ap = R12*R1/R12Ap; |
---|
464 | G4double R23Ap = R22*R2/R22Ap*Pnucl; |
---|
465 | G4double R23dR13 = R23Ap/R13Ap; |
---|
466 | G4double R12Apd = 2./R12Ap; |
---|
467 | G4double R22Apd = 2./R22Ap; |
---|
468 | G4double R122Apd = (R12Apd+R22Apd)/2; |
---|
469 | G4double dR0H = dR0*HadrEnergy/4; |
---|
470 | G4double DDSec1p = DDSect2+DDSect3*std::log(dR0H/R1); |
---|
471 | G4double DDSec2p = DDSect2+DDSect3*std::log(dR0H/std::sqrt((R12+R22)/2)); |
---|
472 | G4double DDSec3p = DDSect2+DDSect3*std::log(dR0H/R2); |
---|
473 | G4double Norm = (R12*R1-Pnucl*R22*R2)*Aeff; // Some questionable norming |
---|
474 | G4double R13 = R12*R1/R12B; |
---|
475 | G4double R23 = Pnucl*R22*R2/R22B; |
---|
476 | G4double norFac = Stot/twopi/Norm; // totCS (in GeV^-2) is used |
---|
477 | G4double Unucl = norFac*R13; |
---|
478 | G4double UnuclScr = norFac*R13Ap; |
---|
479 | G4double SinFi = HadrReIm/Rho2; |
---|
480 | G4double FiH = std::asin(SinFi); |
---|
481 | G4double N = -1.; |
---|
482 | G4double N2 = R23/R13; |
---|
483 | G4double ImElA0 = 0.; |
---|
484 | G4double ReElA0 = 0.; |
---|
485 | G4double Tot1 = 0.; |
---|
486 | for(G4int i=1; i<=A; i++) // @@ Make separately for n and p |
---|
487 | { |
---|
488 | N *= -Unucl*(A-i+1)*Rho2/i; // Includes total cross-section |
---|
489 | G4double N4 = 1.; |
---|
490 | G4double Prod1 = std::exp(-Q2*R12B/i/4)*R12B/i; // Includes the slope |
---|
491 | G4double medTot = R12B/i; |
---|
492 | for(G4int l=1; l<=i; l++) |
---|
493 | { |
---|
494 | G4double exp1 = l/R22B+(i-l)/R12B; |
---|
495 | N4 *= -(i-l+1)*N2/l; |
---|
496 | G4double Nexp = N4/exp1; |
---|
497 | Prod1 += Nexp*std::exp(-Q2/exp1/4); |
---|
498 | medTot += Nexp; |
---|
499 | } |
---|
500 | G4double FiHi = FiH*i; |
---|
501 | G4double nCos = N*std::cos(FiHi); |
---|
502 | ReElA0 += Prod1*N*std::sin(FiHi); |
---|
503 | ImElA0 += Prod1*nCos; |
---|
504 | Tot1 += medTot*nCos; |
---|
505 | if(std::abs(Prod1*N/ImElA0) < eps) break; |
---|
506 | } |
---|
507 | ImElA0 *= piMG; // The amplitude in mb |
---|
508 | ReElA0 *= piMG; // The amplitude in mb |
---|
509 | Tot1 *= piMG+piMG; |
---|
510 | G4double N1p = 1.; |
---|
511 | G4double R13Ap1 = DDSec1p*R13Ap*R13Ap/2; |
---|
512 | G4double R22Ap1 = DDSec2p*R13Ap*R23Ap; |
---|
513 | G4double R23Ap1 = DDSec3p*R23Ap*R23Ap*R22Ap/2; |
---|
514 | G4double R13Ap2 = R13Ap1*R12Ap/4; |
---|
515 | G4double R22Ap2 = R22Ap1/R122Apd/2; |
---|
516 | G4double R23Ap2 = R23Ap1*R22Ap/4; |
---|
517 | G4double Q28 = Q2/8; |
---|
518 | G4double Q24 = Q28+Q28; |
---|
519 | G4double Din1 = R13Ap2*std::exp(-Q28*R12Ap)-R22Ap2*std::exp(-(Q28+Q28)/R122Apd)+ |
---|
520 | R23Ap2*std::exp(-Q28*R22Ap); // i=0 start value |
---|
521 | G4double DTot1 = R13Ap2-R22Ap2+R23Ap2; // i=0 start value |
---|
522 | if (A<96) |
---|
523 | { |
---|
524 | for(G4int i=1; i<A-1; i++) |
---|
525 | { |
---|
526 | N1p *= -UnuclScr*(A-i-1)/i*Rho2; |
---|
527 | G4double N2p = 1.; |
---|
528 | G4double Din2 = 0.; |
---|
529 | G4double DmedTot = 0.; |
---|
530 | G4double BinCoeff = 1.; // Start for Binomial Coefficient |
---|
531 | for(G4int l = 0; l<=i; l++) |
---|
532 | { |
---|
533 | if(l) BinCoeff *= (i-l+1)/l; |
---|
534 | G4double exp1 = l/R22B+(i-l)/R12B; |
---|
535 | G4double exp1p = exp1+R12Apd; |
---|
536 | G4double exp2p = exp1+R122Apd; |
---|
537 | G4double exp3p = exp1+R22Apd; |
---|
538 | G4double NBinC = N2p*BinCoeff; |
---|
539 | Din2 += NBinC*(R13Ap1/exp1p*std::exp(-Q24/exp1p)- |
---|
540 | R22Ap1/exp2p*std::exp(-Q24/exp2p)+ |
---|
541 | R23Ap1/exp3p*std::exp(-Q24/exp3p)); |
---|
542 | DmedTot += NBinC * (R13Ap1/exp1p - R22Ap1/exp2p + R23Ap1/exp3p); |
---|
543 | N2p *= -R23dR13; |
---|
544 | } |
---|
545 | G4double comFac = N1p*Mnoj[i]/(i+2)/(i+1)*std::cos(FiH*i); |
---|
546 | Din1 += Din2*comFac; |
---|
547 | DTot1 += DmedTot*comFac; |
---|
548 | if(std::abs(Din2*N1p/Din1) < eps) break; |
---|
549 | } |
---|
550 | G4double cFac = A*(A-1)*4/Norm/Norm; |
---|
551 | Din1 *= -cFac; |
---|
552 | DTot1 *= cFac; |
---|
553 | } |
---|
554 | //else Din1=0.; |
---|
555 | |
---|
556 | G4double Corr0 = Tot00/Tot1; |
---|
557 | ImElA0 *= Corr0; |
---|
558 | |
---|
559 | //return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*CMMom*CMMom*mb2G2/4/pi2; // ds/do |
---|
560 | return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*mb2G2/(twopi+twopi); |
---|
561 | } // End of DiffElasticCS |
---|
562 | |
---|
563 | G4double CHIPSDiffElCS(G4int hPDG, G4double HadrMom, G4int Z, G4int A, G4double aQ2) // MeV |
---|
564 | {// ================================================================= |
---|
565 | static const G4double eps = 0.000001; // Accuracy of calculations |
---|
566 | static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
567 | static const G4double piMG = pi/mb2G2; // PiTrans from GeV^-2 to mb |
---|
568 | G4double p = HadrMom/1000.; // Momentum (GeV, inputParam) |
---|
569 | G4double p2 = p*p; |
---|
570 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
571 | G4double MassH2 = MassH*MassH; |
---|
572 | G4double HadrEnergy = std::sqrt(p2+MassH2); // Tot energy in GeV |
---|
573 | //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl; |
---|
574 | //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)"); |
---|
575 | CalculateParameters(hPDG, HadrMom); // ? |
---|
576 | CalculateIntegralCS(A, HadrEnergy); // ? |
---|
577 | CalculateNuclearParameters(A); // ? |
---|
578 | G4double Q2 = aQ2/1000/1000; // in GeV^2 |
---|
579 | //G4doubleMassN = A*0.938; // @@ This is bad! |
---|
580 | G4double MassN = A*0.9315; // @@ Atomic Unit is bad too |
---|
581 | if(G4NucleiPropertiesTable::IsInTable(Z,A)) |
---|
582 | MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.; // Geant4 NuclearMass in GeV |
---|
583 | G4double MassN2 = MassN*MassN; |
---|
584 | G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s |
---|
585 | G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // CM energy of a hadron |
---|
586 | G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum |
---|
587 | //G4cout<<"P="<<CMMom<<",E="<<HadrEnergy<<",N="<<MassN2<<",h="<<MassH2<<",p="<<p<<G4endl; |
---|
588 | G4double p3 = p2*p; |
---|
589 | G4double p4 = p2*p2; |
---|
590 | G4double sp = std::sqrt(p); |
---|
591 | G4double p2s = p2*sp; |
---|
592 | G4double ap = std::log(p); |
---|
593 | G4double dl = ap-3.; |
---|
594 | G4double dl2 = dl*dl; |
---|
595 | G4double shPPTot = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb |
---|
596 | G4double shPNTot = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p |
---|
597 | +(38.+.3*dl2+8./p)/(1.+1.2/p3); // SigPP in mb |
---|
598 | G4double hNTot = (Z*shPPTot+(A-Z)*shPNTot)/A; // SighN in mb |
---|
599 | G4double Stot = hNTot*mb2G2; // in GeV^-2 |
---|
600 | G4double shPPSl = 8.*std::pow(p,.055)/(1.+3.64/p4); // PPslope in GeV^-2 |
---|
601 | G4double shPNSl = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2 |
---|
602 | G4double Bhad = (Z*shPPSl+(A-Z)*shPNSl)/A; // B-slope in GeV^-2 |
---|
603 | G4double hNReIm = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit |
---|
604 | G4double Asq = 1+hNReIm*hNReIm; // |M|^2/(ImM)^2 |
---|
605 | G4double Rho2 = std::sqrt(Asq); // M/ImM |
---|
606 | G4double r1 = 3.9*std::pow(A-1.,.309); // Positive diffractionalRadius |
---|
607 | G4double r2 = 2.*std::pow(A,.36); // Negative diffraction radius |
---|
608 | G4double pN = Pnucl; // Dubna value |
---|
609 | //G4double pN = .4; // Screaning factor |
---|
610 | G4double Ae = Aeff; // Dubna value |
---|
611 | //G4double Ae = .75; // Normalization |
---|
612 | G4double R12 = r1*r1; |
---|
613 | G4double R22 = r2*r2; |
---|
614 | G4double R1C = R12*r1; |
---|
615 | G4double R2C = R22*r2; |
---|
616 | G4double R12B = R12+Bhad+Bhad; // Slope is used |
---|
617 | G4double R22B = R22+Bhad+Bhad; // Slope is used |
---|
618 | G4double Norm = (R1C-pN*R2C)*Ae; // ScreanFac & NormFac are used |
---|
619 | G4double R13 = R1C/R12B; // Slope is used |
---|
620 | G4double R23 = pN*R2C/R22B; // Slope & ScreanFact are used |
---|
621 | G4double norFac = Stot/twopi/Norm; // totCS (in GeV^-2) is used |
---|
622 | G4double Unucl = norFac*R13; |
---|
623 | G4double SinFi = hNReIm/Rho2; // Real part |
---|
624 | G4double FiH = std::asin(SinFi); |
---|
625 | G4double N = -1.; |
---|
626 | G4double N2 = R23/R13; // Slope is used |
---|
627 | G4double ImElA0 = 0.; |
---|
628 | G4double ReElA0 = 0.; |
---|
629 | G4double Tot1 = 0.; |
---|
630 | for(G4int i=1; i<=A; i++) |
---|
631 | { |
---|
632 | N *= -Unucl*(A-i+1)*Rho2/i; // TotCS is used |
---|
633 | G4double N4 = 1.; |
---|
634 | G4double Prod1 = std::exp(-Q2*R12B/i/4)*R12B/i; // Slope is used |
---|
635 | G4double medTot = R12B/i; // Slope is used |
---|
636 | for(G4int l=1; l<=i; l++) |
---|
637 | { |
---|
638 | G4double exp1 = l/R22B+(i-l)/R12B; // Slope is used |
---|
639 | N4 *= -(i-l+1)*N2/l; // Slope is used |
---|
640 | G4double Nexp = N4/exp1; |
---|
641 | Prod1 += Nexp*std::exp(-Q2/exp1/4); |
---|
642 | medTot += Nexp; |
---|
643 | } |
---|
644 | G4double FiHi = FiH*i; |
---|
645 | G4double nCos = N*std::cos(FiHi); |
---|
646 | ReElA0 += Prod1*N*std::sin(FiHi); |
---|
647 | ImElA0 += Prod1*nCos; |
---|
648 | Tot1 += medTot*nCos; |
---|
649 | if(std::abs(Prod1*N/ImElA0) < eps) break; |
---|
650 | } |
---|
651 | ImElA0 *= piMG; // The amplitude in mb |
---|
652 | ReElA0 *= piMG; // The amplitude in mb |
---|
653 | Tot1 *= piMG+piMG; |
---|
654 | G4double Corr0 = Tot00/Tot1; |
---|
655 | ImElA0 *= Corr0; |
---|
656 | |
---|
657 | //return (ReElA0*ReElA0+ImElA0*ImElA0)*CMMom*CMMom*mb2G2/4/pi2; // ds/do |
---|
658 | return (ReElA0*ReElA0+ImElA0*ImElA0)*mb2G2/(twopi+twopi); |
---|
659 | } // End of DiffElasticCS |
---|
660 | |
---|
661 | G4double CHIPS_Tb(G4int A, G4double b) // T(b) in fm-2 |
---|
662 | { |
---|
663 | static G4int Am=0; // Associative memory value A |
---|
664 | static G4double B=0.; // Effective edge |
---|
665 | static G4double C=0.; // Normalization factor |
---|
666 | static G4double D=0.; // Slope parameter |
---|
667 | if(A!=Am) |
---|
668 | { |
---|
669 | B=.0008*A*A; // no units |
---|
670 | D=.42*std::pow(A,-.26); // fm^-2 |
---|
671 | C=A*D/pi/std::log(1.+B); // fm^-2 |
---|
672 | } |
---|
673 | G4double E=B*std::exp(-D*b*b); |
---|
674 | return C*E/(1.+E); // T(b) in fm^-2 |
---|
675 | } |
---|
676 | |
---|
677 | G4double Thickness(G4int A, G4double b, G4double R) // T(b) in fm^-2 |
---|
678 | {// ========================================== // rAmax=2.5 is a GlobStatConstPar |
---|
679 | static const G4double bTh = .545; // Surface thicknes |
---|
680 | static const G4double bTh2 = bTh*bTh; // SquaredSurface thicknes |
---|
681 | G4double b2 = b*b; // Squared impact parameter |
---|
682 | G4double dr = rAmax*R/(NptB-1); // Step of integration in z |
---|
683 | G4double R2 = R*R; // Working parameter |
---|
684 | G4double Norm = .75/pi/R2/R/(1.+bTh2*pi2/R2); // Corrected Volume (?) |
---|
685 | G4double r = -dr; // Pre value for the radius |
---|
686 | G4double SumZ = 0.; // Prototype of the result |
---|
687 | //G4double SumN = 0.; // is not used |
---|
688 | for(G4int k=0; k<NptB; k++) |
---|
689 | { |
---|
690 | r += dr; |
---|
691 | SumZ += 1./(1.+std::exp((std::sqrt(b2+r*r)-R)/bTh)); // Woods-Saxon density |
---|
692 | //SumN += r*r/(1.+std::exp((r-rAfm)/bTh)); // is not used |
---|
693 | } |
---|
694 | //SumN *= (twopi+twopi)*dr*Norm; // is not used |
---|
695 | return SumZ*(A+A)*dr*Norm; // T(b) in fm^-2 |
---|
696 | } // End of Thickness |
---|
697 | |
---|
698 | G4double CoherentDifElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) //in MeV |
---|
699 | {// ========================================================================= |
---|
700 | //static const G4double eps = 0.000001; // CalculationAccuracy@@NotUsed |
---|
701 | //static const G4double mN = .9315; // Atomic Unit GeV |
---|
702 | static const G4double mN = .938; // Atomic Unit GeV |
---|
703 | static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
704 | static const G4double incrf = 10; // Increase factor of radius for intergation |
---|
705 | static const G4double m2G10 = incrf*mb2G2; // Big conversion factor |
---|
706 | static const G4double s2G10 = std::sqrt(m2G10); // sqrt Big conversion factor |
---|
707 | G4double ReIntegrand[NptB],ImIntegrand[NptB],Thick[NptB]; // Calculated arrays |
---|
708 | G4QBesIKJY QI0(BessI0); // I0 Bessel function |
---|
709 | G4QBesIKJY QJ0(BessJ0); // J0 Bessel function |
---|
710 | G4double hMom = HadrMom/1000.; // Momentum (GeV, inputParam) |
---|
711 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
712 | G4double MassH2 = MassH*MassH; |
---|
713 | G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2); // Tot energy in GeV |
---|
714 | if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3"); |
---|
715 | if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!"); |
---|
716 | //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)"); |
---|
717 | G4double Re = std::pow(A, 0.33333333); |
---|
718 | G4double Re2 = Re*Re; |
---|
719 | CalculateParameters(hPDG, HadrMom); |
---|
720 | G4double Q2 = aQ2/1000/1000; // GeV |
---|
721 | // Individual corrections for nuclei which had data |
---|
722 | if (A==208) r0 = 1.125; |
---|
723 | else if(A==90) r0 = 1.12; |
---|
724 | else if(A==64) r0 = 1.1; |
---|
725 | else if(A==58) r0 = 1.09; |
---|
726 | else if(A==48) r0 = 1.07; |
---|
727 | else if(A==40) r0 = 1.15; |
---|
728 | else if(A==28) r0 = 0.93; |
---|
729 | else if(A==16) r0 = 0.92; |
---|
730 | else if(A==12) r0 = 0.80; |
---|
731 | else r0 = 1.16*(1.-1.16/Re2); // For other nuclei which have not been tested |
---|
732 | G4double MassN = A*0.9315; // @@ Atomic Unit is bad too |
---|
733 | G4double MassN2 = MassN*MassN; |
---|
734 | G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2; // Mondelstam S |
---|
735 | G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // Hadron CM Energy |
---|
736 | G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum |
---|
737 | G4double rAfm = r0*Re; // Nuclear radius |
---|
738 | G4double rAGeV = rAfm*s2G10; // Big integration radius |
---|
739 | G4double stepB = rAmax*rAGeV/(NptB-1); // dr step of integration |
---|
740 | G4double hTotG2 = HadrTot*mb2G2; |
---|
741 | G4double ReSum = 0.; |
---|
742 | G4double ImSum = 0.; |
---|
743 | G4double ValB = -stepB; |
---|
744 | for(G4int i=0; i<NptB; i++) |
---|
745 | { |
---|
746 | ValB += stepB; // An incident parameter |
---|
747 | G4double ValB2 = ValB*ValB; // A working value |
---|
748 | G4double IPH = ValB/HadrSlope; // A working value |
---|
749 | G4double dHS = HadrSlope+HadrSlope; // A working value |
---|
750 | G4double Integ = 0.; |
---|
751 | G4double ValS = 0.; |
---|
752 | for(G4int j=1; j<NptB; j++) |
---|
753 | { |
---|
754 | ValS += stepB; // Impact parameter GeV^-1 |
---|
755 | if(!i) Thick[j] = Thickness(A,ValS/s2G10,rAfm)/m2G10; // Calculate only once |
---|
756 | G4double FunS = IPH*ValS; |
---|
757 | if(FunS > 320.) break; |
---|
758 | Integ += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Thick[j]; |
---|
759 | } |
---|
760 | G4double InExp = -hTotG2*Integ*stepB/dHS; |
---|
761 | G4double expB = std::exp(InExp); |
---|
762 | G4double HRIE = HadrReIm*InExp; |
---|
763 | ReIntegrand[i] = (1.-expB*std::cos(HRIE)); // Real part of the amplitude |
---|
764 | ImIntegrand[i] = expB*std::sin(HRIE); // Imaginary part of the amplit |
---|
765 | } |
---|
766 | InCoh = 0.; // incohirent (quasi-elastic) |
---|
767 | ValB = -stepB; |
---|
768 | for(G4int k=0; k<NptB; k++) // Third integration (?) |
---|
769 | { |
---|
770 | ValB += stepB; |
---|
771 | InCoh += Thick[k]*ValB*std::exp(-hTotG2*Thick[k]); |
---|
772 | G4double J0qb = QJ0(std::sqrt(Q2)*ValB)*ValB; // Bessel0(Q*b) |
---|
773 | ReSum += J0qb*ReIntegrand[k]; |
---|
774 | ImSum += J0qb*ImIntegrand[k]; |
---|
775 | } |
---|
776 | //G4cout<<"GHAD:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2=" |
---|
777 | // <<Q2<<",m="<<mb2G2<<G4endl; |
---|
778 | InCoh *= stepB*hTotG2*hTotG2*(1.+HadrReIm*HadrReIm)*std::exp(-HadrSlope*Q2)/8/mb2G2; |
---|
779 | //G4cout<<"GHAD:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl; |
---|
780 | //return (ReSum*ReSum+ImSum*ImSum)*m2G10*CMMom*CMMom*stepB*stepB/twopi; // ds/do |
---|
781 | return (ReSum*ReSum+ImSum*ImSum)*m2G10*stepB*stepB/12; // ds/dt |
---|
782 | } |
---|
783 | |
---|
784 | G4double CHIPSDifElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2) |
---|
785 | {// ============================================================================ |
---|
786 | static const G4int Npb = 500; // A#of intergation points |
---|
787 | //static const G4double mN = .9315; // Atomic Unit GeV |
---|
788 | static const G4double mN = .938; // Atomic Unit GeV |
---|
789 | static const G4double hc2 = .3893793; // Transform from GeV^-2 to mb |
---|
790 | static const G4double mb2G2 = 1./hc2; // Transform from mb to GeV^-2 |
---|
791 | static const G4double f22mb = 10; // Transform from fermi^2 to mb |
---|
792 | static const G4double f22G2 = f22mb*mb2G2; // Transform from fm2 to GeV^-2 |
---|
793 | static const G4double f2Gm1 = std::sqrt(f22G2); // Transform from fm to GeV^-1 |
---|
794 | G4double Re = std::pow(A,.33333333); // A-dep coefficient |
---|
795 | G4double Lim = 100*Re; // Integration accuracy limit |
---|
796 | G4double Tb[Npb]; // Calculated T(b) array |
---|
797 | G4QBesIKJY QI0(BessI0); // I0 Bessel function |
---|
798 | G4QBesIKJY QJ0(BessJ0); // J0 Bessel function |
---|
799 | G4double p = hMom/1000.; // Momentum (GeV, inputParam) |
---|
800 | G4double p2 = p*p; |
---|
801 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
802 | G4double MassH2 = MassH*MassH; // Squared mass of the hadron |
---|
803 | G4double hEnergy = std::sqrt(p2+MassH2); // Tot energy in GeV |
---|
804 | G4double Q2 = aQ2/1000000.; // -t in GeV |
---|
805 | G4double MassN = A*0.9315; // @@ Atomic Unit is bad too |
---|
806 | if(G4NucleiPropertiesTable::IsInTable(Z,A)) |
---|
807 | MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.; // Geant4 NuclearMass in GeV |
---|
808 | G4double MassN2 = MassN*MassN; // Squared mass of the target |
---|
809 | G4double S = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s |
---|
810 | G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // Hadron CM Energy |
---|
811 | G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum (to norm CS) |
---|
812 | //G4cout<<"2: P="<<CMMom<<",E="<<hEnergy<<",N="<<MassN2<<",h="<<MassH2<<",p="<<p<<G4endl; |
---|
813 | // The mean value of the total can be used |
---|
814 | G4double p3 = p2*p; |
---|
815 | G4double p4 = p2*p2; |
---|
816 | G4double sp = std::sqrt(p); |
---|
817 | G4double p2s = p2*sp; |
---|
818 | G4double ap = std::log(p); |
---|
819 | G4double dl = ap-3.; |
---|
820 | G4double dl2 = dl*dl; |
---|
821 | G4double shPPTot = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb |
---|
822 | G4double shPNTot = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p |
---|
823 | +(38.+.3*dl2+8./p)/(1.+1.2/p3); // SigPP in mb |
---|
824 | G4double shNTot = (Z*shPPTot+(A-Z)*shPNTot)/A; // SighN in mb |
---|
825 | #ifdef debug |
---|
826 | G4cout<<"CHIPS:SI,p="<<p<<",n="<<shNTot<<",P="<<shPPTot<<",N="<<shPNTot |
---|
827 | <<",Z="<<Z<<",A="<<A<<G4endl; |
---|
828 | #endif |
---|
829 | G4double shPPSl = 8.*std::pow(p,.055)/(1.+3.64/p4); // PPslope in GeV^-2 |
---|
830 | G4double shPNSl = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2 |
---|
831 | G4double shNSl = (Z*shPPSl+(A-Z)*shPNSl)/A; // B-slope in GeV^-2 |
---|
832 | #ifdef debug |
---|
833 | G4cout<<"CHIPS:SL,n="<<shNSl<<",P="<<shPPSl<<",N="<<shPNSl<<G4endl; |
---|
834 | #endif |
---|
835 | G4double shNReIm = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit |
---|
836 | //G4cout<<"CHPS: s="<<S<<",T="<<shNTot<<",R="<<shNReIm<<",B="<<shNSl<<G4endl; |
---|
837 | // @@ End of temporary ^^^^^^^ |
---|
838 | G4double dHS = shNSl+shNSl; // Working: doubled B-slope |
---|
839 | G4double stepB = (Re+Re+2.7)*f2Gm1/(Npb-1); // in GeV^-1, step of integral |
---|
840 | G4double hTotG2 = shNTot*mb2G2; // sigma_hN in GeV^-2 |
---|
841 | G4double ReSum = 0.; // Integration of RePart of Amp |
---|
842 | G4double ImSum = 0.; // Integration of ImPart of Amp |
---|
843 | G4double ValB = -stepB; |
---|
844 | for(G4int i=0; i<Npb; i++) // First integration over b |
---|
845 | { |
---|
846 | ValB += stepB; // An incident parameter |
---|
847 | G4double ValB2 = ValB*ValB; // A working value b^2 |
---|
848 | G4double IPH = ValB/shNSl; // A working value slope |
---|
849 | G4double Integ = 0.; // Integral over ImpactParam. |
---|
850 | G4double ValS = 0.; // Prototype of ImpactParameter |
---|
851 | for(G4int j=1; j<Npb; j++) // Second integration over b |
---|
852 | { |
---|
853 | ValS += stepB; // back to fm // Impact parameter GeV^-1 |
---|
854 | if(!i) Tb[j] = CHIPS_Tb(A,ValS/f2Gm1)/f22G2; // GeV^2, calculate only once |
---|
855 | //if(!i) Tb[j] = Thickness(A,ValS/f2Gm1,rAfm)/f22G2; // Calculate T(b) only once |
---|
856 | G4double FunS = IPH*ValS; // b1*b2/slope |
---|
857 | if(FunS > Lim) break; // To avoid NAN |
---|
858 | Integ += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Tb[j]; // BessI0 |
---|
859 | } |
---|
860 | G4double InExp = -hTotG2*Integ*stepB/dHS; // Integrated absorption |
---|
861 | G4double expB = std::exp(InExp); // Exponential absorption |
---|
862 | G4double HRIE = shNReIm*InExp; // Phase shift |
---|
863 | G4double J0qb = QJ0(std::sqrt(Q2)*ValB)*ValB; // Bessel0(Q*b) |
---|
864 | ReSum += J0qb*(1.-expB*std::cos(HRIE)); |
---|
865 | ImSum += J0qb*expB*std::sin(HRIE); |
---|
866 | } |
---|
867 | //G4cout<<"CHPS:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl; |
---|
868 | //return (ReSum*ReSum+ImSum*ImSum)*mb2G2*CMMom*CMMom*stepB*stepB/twopi; |
---|
869 | //return (ReSum*ReSum+ImSum*ImSum)*f22G2*CMMom*CMMom*stepB*stepB/twopi; // ds/do |
---|
870 | return (ReSum*ReSum+ImSum*ImSum)*f22G2*stepB*stepB/12; // ds/dt |
---|
871 | } |
---|
872 | |
---|
873 | // Separate quasielastic calculation |
---|
874 | G4double CHIPSDifQuasiElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2) |
---|
875 | {// ================================================================================= |
---|
876 | static const G4int Npb = 500; // A#of intergation points |
---|
877 | static const G4double hc2 = .3893793; // Transform from GeV^-2 to mb |
---|
878 | static const G4double mb2G2 = 1./hc2; // Transform from mb to GeV^-2 |
---|
879 | //static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
880 | static const G4double f22mb = 10; // Transform from fermi^2 to mb |
---|
881 | static const G4double f22G2 = f22mb*mb2G2; // Transform from fm2 to GeV^-2 |
---|
882 | static const G4double f2Gm1 = std::sqrt(f22G2); // Transform from fm to GeV^-1 |
---|
883 | G4double p = hMom/1000.; // Momentum (GeV, inputParam) |
---|
884 | G4double p2 = p*p; |
---|
885 | G4double Q2 = aQ2/1000000.; // -t in GeV |
---|
886 | // @@ Temporary only for nucleons |
---|
887 | G4double p3 = p2*p; |
---|
888 | G4double p4 = p2*p2; |
---|
889 | G4double sp = std::sqrt(p); |
---|
890 | G4double p2s = p2*sp; |
---|
891 | G4double ap = std::log(p); |
---|
892 | G4double dl = ap-3.; |
---|
893 | G4double dl2 = dl*dl; |
---|
894 | G4double shPPTot = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb |
---|
895 | G4double shPNTot = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p |
---|
896 | +(38.+.3*dl2+8./p)/(1.+1.2/p3); // SigPP in mb |
---|
897 | G4double shNTot = (Z*shPPTot+(A-Z)*shPNTot)/A; // SighN in mb |
---|
898 | G4double shPPSl = 8.*std::pow(p,.055)/(1.+3.64/p4); // PPslope in GeV^-2 |
---|
899 | G4double shPNSl = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2 |
---|
900 | G4double shNSl = (Z*shPPSl+(A-Z)*shPNSl)/A; // B-slope in GeV^-2 |
---|
901 | G4double shNReIm = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit |
---|
902 | // @@ End of temporary ^^^^^^^ |
---|
903 | G4double Re = std::pow(A,.33333333); // A-dep coefficient |
---|
904 | G4double stepB = (Re+Re+2.7)*f2Gm1/(Npb-1); // in GeV^-1, step of integral |
---|
905 | G4double InQE = 0.; // Quasielastic integral |
---|
906 | G4double ValB = -stepB; // Impact parameter |
---|
907 | G4double hTotG2 = shNTot*mb2G2; // sigma_hN in GeV^-2 |
---|
908 | for(G4int k=0; k<Npb; k++) |
---|
909 | { |
---|
910 | ValB += stepB; |
---|
911 | //G4double Tb = Thickness(A,ValB/f2Gm1,rAfm)/f22G2; // Calculate T(b) only once |
---|
912 | G4double Tb = CHIPS_Tb(A,ValB/f2Gm1)/f22G2; // GeV^2, calculate only once |
---|
913 | InQE += Tb*ValB*std::exp(-hTotG2*Tb); |
---|
914 | } |
---|
915 | //G4cout<<"CHPS:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2=" |
---|
916 | // <<Q2<<",m="<<mb2G2<<G4endl; |
---|
917 | return InQE*stepB*hTotG2*hTotG2*(1.+shNReIm*shNReIm)*std::exp(-shNSl*Q2)/8/mb2G2; |
---|
918 | } |
---|
919 | #endif |
---|
920 | |
---|
921 | int main() |
---|
922 | { |
---|
923 | #ifdef bestest |
---|
924 | // *** Test of CHIPS implementation of Bessel functions (accuracy 1.E-8 and better) *** |
---|
925 | const G4int nj=21; |
---|
926 | const G4int ny=16; |
---|
927 | const G4int ni=20; |
---|
928 | const G4int nk=20; |
---|
929 | const G4double jX[nj]= |
---|
930 | {-5.,-4.,-3.,-2.,-1.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.}; |
---|
931 | const G4double yX[ny]={.1,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.}; |
---|
932 | const G4double iX[ni]= |
---|
933 | {0.,.2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,6.,8.,10.}; |
---|
934 | const G4double kX[nk]= |
---|
935 | {.1,.2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,6.,8.,10.}; |
---|
936 | const G4double vBJ0[nj]= |
---|
937 | {-.1775968,-.3971498,-.2600520,0.2238908,0.7651976,1.0000000,0.7651977, |
---|
938 | 0.2238908,-.2600520,-.3971498,-.1775968,0.1506453,0.3000793,0.1716508, |
---|
939 | -.0903336,-.2459358,-.1711903,0.0476893,0.2069261,0.1710735,-.0142245}; |
---|
940 | const G4double vBJ1[nj]= |
---|
941 | {0.3275791,0.0660433,-.3390590,-.5767248,-.4400506,0.0000000,0.4400506, |
---|
942 | 0.5767248,0.3390590,-.0660433,-.3275791,-.2766839,-.0046828,0.2346364, |
---|
943 | 0.2453118,0.0434728,-.1767853,-.2234471,-.0703181,0.1333752,0.2051040}; |
---|
944 | const G4double vBY0[ny]= |
---|
945 | {-1.534239,0.0882570,.51037567,.37685001,-.0169407,-.3085176,-.2881947,-.0259497, |
---|
946 | 0.2235215,0.2499367,0.0556712,-.1688473,-.2252373,-.0782079,0.1271926,0.2054743}; |
---|
947 | const G4double vBY1[ny]= |
---|
948 | {-6.458951,-.7812128,-.1070324,0.3246744,0.3979257,0.1478631,-.1750103,-.3026672, |
---|
949 | -.1580605,0.1043146,0.2490154,0.1637055,-.0570992,-.2100814,-.1666448,0.0210736}; |
---|
950 | const G4double vBI0[ni]={1.0000000,1.0100250,1.0404018,1.0920453,1.1665149, |
---|
951 | 1.2660658,1.3937256,1.5533951,1.7499807,1.9895593, |
---|
952 | 2.2795852,3.2898391,4.8807925,7.3782035,11.301922, |
---|
953 | 17.481172,27.239871,67.234406,427.56411,2815.7167}; |
---|
954 | const G4double vBI1[ni]={.00000000,.10050083,.20402675,.31370403,.43286480, |
---|
955 | .56515912,.71467794,.88609197,1.0848107,1.3171674, |
---|
956 | 1.5906369,2.5167163,3.9533700,6.2058350,9.7594652, |
---|
957 | 15.389221,24.335643,61.341937,399.87313,2670.9883}; |
---|
958 | const G4double vBK0[nk]={2.4270690,1.7527038,1.1145291,.77752208,.56534710, |
---|
959 | .42102445,.31850821,.24365506,.18795475,.14593140, |
---|
960 | .11389387,.062347553,.0347395,.019598897,.011159676, |
---|
961 | .0063998572,.0036910983,.0012439943,.00014647071,.000017780062}; |
---|
962 | const G4double vBK1[nk]={9.8538451,4.7759725,2.1843544,1.3028349,.86178163, |
---|
963 | .60190724,.43459241,.32083589,.24063392,.18262309, |
---|
964 | .13986588,.073890816,.040156431,.022239393,.012483499, |
---|
965 | .0070780949,.0040446134,.0013439197,.00015536921,.000018648773}; |
---|
966 | //Bessel J0('J',0); // CLHEP J/Y is possible insead of 'J' |
---|
967 | //Bessel J1('J',1); |
---|
968 | G4QBesIKJY QI0(BessI0); // CHIPS I/K/J/Y 0/1 Bessel functions |
---|
969 | G4QBesIKJY QI1(BessI1); |
---|
970 | G4QBesIKJY QJ0(BessJ0); |
---|
971 | G4QBesIKJY QJ1(BessJ1); |
---|
972 | G4QBesIKJY QK0(BessK0); |
---|
973 | G4QBesIKJY QK1(BessK1); |
---|
974 | G4QBesIKJY QY0(BessY0); |
---|
975 | G4QBesIKJY QY1(BessY1); |
---|
976 | G4cout<<"G4GCS: ***J0*** Test: ================================================"<<G4endl; |
---|
977 | for(G4int j0=0; j0<nj; j0++) |
---|
978 | { |
---|
979 | G4double F=QJ0(jX[j0]); |
---|
980 | G4double d=F-vBJ0[j0]; |
---|
981 | G4double r=std::abs(d/F); |
---|
982 | G4cout<<"G4GCS: x="<<jX[j0]<<", J0="<<F<<" - "<<vBJ0[j0]<<" = "<<d<<", r="<<r<<G4endl; |
---|
983 | } |
---|
984 | G4cout<<"G4GCS: ***J1*** Test: ================================================"<<G4endl; |
---|
985 | for(G4int j1=0; j1<nj; j1++) |
---|
986 | { |
---|
987 | G4double F=QJ1(jX[j1]); |
---|
988 | G4double d=F-vBJ1[j1]; |
---|
989 | G4double r=std::abs(d/F); |
---|
990 | G4cout<<"G4GCS: x="<<jX[j1]<<", J1="<<F<<" - "<<vBJ1[j1]<<" = "<<d<<", r="<<r<<G4endl; |
---|
991 | } |
---|
992 | G4cout<<"G4GCS: ***Y0*** Test: ================================================"<<G4endl; |
---|
993 | for(G4int y0=0; y0<ny; y0++) |
---|
994 | { |
---|
995 | G4double F=QY0(yX[y0]); |
---|
996 | G4double d=F-vBY0[y0]; |
---|
997 | G4double r=std::abs(d/F); |
---|
998 | G4cout<<"G4GCS: x="<<yX[y0]<<", Y0="<<F<<" - "<<vBY0[y0]<<" = "<<d<<", r="<<r<<G4endl; |
---|
999 | } |
---|
1000 | G4cout<<"G4GCS: ***Y1*** Test: ================================================"<<G4endl; |
---|
1001 | for(G4int y1=0; y1<ny; y1++) |
---|
1002 | { |
---|
1003 | G4double F=QY1(yX[y1]); |
---|
1004 | G4double d=F-vBY1[y1]; |
---|
1005 | G4double r=std::abs(d/F); |
---|
1006 | G4cout<<"G4GCS: x="<<yX[y1]<<", Y1="<<F<<" - "<<vBY1[y1]<<" = "<<d<<", r="<<r<<G4endl; |
---|
1007 | } |
---|
1008 | G4cout<<"G4GCS: ***I0*** Test: ================================================"<<G4endl; |
---|
1009 | for(G4int i0=0; i0<ni; i0++) |
---|
1010 | { |
---|
1011 | G4double F=QI0(iX[i0]); |
---|
1012 | G4double d=F-vBI0[i0]; |
---|
1013 | G4double r=std::abs(d/F); |
---|
1014 | G4cout<<"G4GCS: x="<<iX[i0]<<", I0="<<F<<" - "<<vBI0[i0]<<" = "<<d<<", r="<<r<<G4endl; |
---|
1015 | } |
---|
1016 | G4cout<<"G4GCS: ***I1*** Test: ================================================"<<G4endl; |
---|
1017 | for(G4int i1=0; i1<ni; i1++) |
---|
1018 | { |
---|
1019 | G4double F=QI1(iX[i1]); |
---|
1020 | G4double d=F-vBI1[i1]; |
---|
1021 | G4double r=std::abs(d/F); |
---|
1022 | G4cout<<"G4GCS: x="<<iX[i1]<<", I1="<<F<<" - "<<vBI1[i1]<<" = "<<d<<", r="<<r<<G4endl; |
---|
1023 | } |
---|
1024 | G4cout<<"G4GCS: ***K0*** Test: ================================================"<<G4endl; |
---|
1025 | for(G4int k0=0; k0<nk; k0++) |
---|
1026 | { |
---|
1027 | G4double F=QK0(kX[k0]); |
---|
1028 | G4double d=F-vBK0[k0]; |
---|
1029 | G4double r=std::abs(d/F); |
---|
1030 | G4cout<<"G4GCS: x="<<kX[k0]<<", I1="<<F<<" - "<<vBK0[k0]<<" = "<<d<<", r="<<r<<G4endl; |
---|
1031 | } |
---|
1032 | G4cout<<"G4GCS: ***K1*** Test: ================================================"<<G4endl; |
---|
1033 | for(G4int k1=0; k1<nk; k1++) |
---|
1034 | { |
---|
1035 | G4double F=QK1(kX[k1]); |
---|
1036 | G4double d=F-vBK1[k1]; |
---|
1037 | G4double r=std::abs(d/F); |
---|
1038 | G4cout<<"G4GCS: x="<<kX[k1]<<", I1="<<F<<" - "<<vBK1[k1]<<" = "<<d<<", r="<<r<<G4endl; |
---|
1039 | } |
---|
1040 | G4cout<<"End .................................................................."<<G4endl; |
---|
1041 | #endif |
---|
1042 | |
---|
1043 | // Test of integrated cross sections |
---|
1044 | //const G4int na=12; |
---|
1045 | const G4int na=1; |
---|
1046 | // |
---|
1047 | const G4int np=7; |
---|
1048 | const G4int nm=1; |
---|
1049 | //// He Be C O Al Ti Ni Cu Sn Ta Pb U |
---|
1050 | //const G4int A[na]={4,9,12,16,27,48,58,64,120,181,207,238}; // A's of target nuclei |
---|
1051 | // He Al Pb |
---|
1052 | const G4int A[na]={119}; // A's of target nuclei |
---|
1053 | // p n pi+ pi- K+ K- antip |
---|
1054 | const G4int pdg[np]={2212,2112,211,-211,321,-321,-2212}; // projectiles |
---|
1055 | const G4double mom[nm]={120.}; // momentum in MeV/c |
---|
1056 | #ifdef integrc |
---|
1057 | for(G4int ip=0; ip<np; ip++) |
---|
1058 | { |
---|
1059 | G4int PDG=pdg[ip]; |
---|
1060 | for(G4int im=0; im<nm; im++) |
---|
1061 | { |
---|
1062 | CalculateParameters(PDG, mom[im]); |
---|
1063 | G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<":A,Tot,Inel,Prod,El,QEl"<<G4endl; |
---|
1064 | for(G4int ia=0; ia<na; ia++) |
---|
1065 | { |
---|
1066 | CalculateIntegralCrossSections(A[ia]); |
---|
1067 | G4cout<<A[ia]<<" "<<TotalCrSec<<" "<<InelCrSec<<" "<<ProdCrSec |
---|
1068 | <<" "<<ElasticCrSec<<" "<<QuasyElasticCrSec<<G4endl; |
---|
1069 | } |
---|
1070 | } |
---|
1071 | } |
---|
1072 | #endif |
---|
1073 | |
---|
1074 | #ifdef eldiffcs |
---|
1075 | const G4double ms=.000001; |
---|
1076 | const G4int nt=200; |
---|
1077 | G4double t[nt]; // Q2=-t in Mev^2 |
---|
1078 | const G4double tMin=200.; |
---|
1079 | const G4double tMax=2000000.; |
---|
1080 | const G4double ltMi=std::log(tMin); |
---|
1081 | const G4double ltMa=std::log(tMax); |
---|
1082 | const G4double dlt=(ltMa-ltMi)/(nt-1); |
---|
1083 | G4double lt=ltMi-dlt; |
---|
1084 | for(G4int it=0; it<nt; it++) |
---|
1085 | { |
---|
1086 | lt+=dlt; |
---|
1087 | t[it]=std::exp(lt); |
---|
1088 | } |
---|
1089 | //// He Be C O Al Ti Ni Cu Sn Ta Pb U |
---|
1090 | //const G4int Z[na]={2,4, 6, 8,13,22,28,29, 50, 73, 82, 92}; // Z's of target nuclei |
---|
1091 | // He Al Pb |
---|
1092 | const G4int Z[na]={50}; // Z's of target nuclei |
---|
1093 | // Test of differential ellastic cross sections |
---|
1094 | Initialize(); |
---|
1095 | //for(G4int ip=0; ip<np; ip++) |
---|
1096 | for(G4int ip=0; ip<1; ip++) |
---|
1097 | { |
---|
1098 | G4int PDG=pdg[ip]; |
---|
1099 | for(G4int im=0; im<nm; im++) |
---|
1100 | { |
---|
1101 | CalculateParameters(PDG, mom[im]); |
---|
1102 | G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<G4endl; |
---|
1103 | for(G4int ia=0; ia<na; ia++) |
---|
1104 | { |
---|
1105 | G4cout<<"G4GCS:-------------------A="<<A[ia]<<G4endl; |
---|
1106 | for(G4int it=0; it<nt; it++) |
---|
1107 | { |
---|
1108 | G4double Sig1 = CHIPSDiffElCS(PDG, mom[im], Z[ia], A[ia], t[it]); |
---|
1109 | //G4double Sig1 = DiffElasticCS(PDG, mom[im], A[ia], t[it]); //Doesn't work |
---|
1110 | G4double Sig2 = CoherentDifElasticCS(PDG, mom[im], A[ia], t[it]); |
---|
1111 | G4double Sig3 = CHIPSDifElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]); |
---|
1112 | G4double CQEl = CHIPSDifQuasiElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]); |
---|
1113 | G4cout<<std::setw(12)<<ms*t[it]<<std::setw(12)<<Sig1<<std::setw(12)<<Sig2 |
---|
1114 | <<std::setw(12)<<InCoh<<std::setw(12)<<Sig3<<std::setw(12)<<CQEl<<G4endl; |
---|
1115 | } |
---|
1116 | } |
---|
1117 | } |
---|
1118 | } |
---|
1119 | #endif |
---|
1120 | return 0; |
---|
1121 | } |
---|