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: G4QANuENuclearCrossSection.cc,v 1.1 2009/11/16 18:15:43 mkossov Exp $ |
---|
28 | // GEANT4 tag $Name: hadr-chips-V09-03-08 $ |
---|
29 | // |
---|
30 | // |
---|
31 | // G4 Physics class: G4QANuENuclearCrossSection for (anti_nu_e,e+) cross sections |
---|
32 | // Created: M.V. Kossov, CERN/ITEP(Moscow), 24-SEP-2007 |
---|
33 | // The last update: M.V. Kossov, CERN/ITEP (Moscow) 24-SEP-2007 |
---|
34 | // |
---|
35 | // **************************************************************************************** |
---|
36 | // ********** This CLASS is temporary moved from the photolepton_hadron directory ********* |
---|
37 | // ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ****** |
---|
38 | // **************************************************************************************** |
---|
39 | //========================================================================================= |
---|
40 | |
---|
41 | //#define debug |
---|
42 | //#define edebug |
---|
43 | //#define pdebug |
---|
44 | //#define ppdebug |
---|
45 | //#define tdebug |
---|
46 | //#define sdebug |
---|
47 | |
---|
48 | #include "G4QANuENuclearCrossSection.hh" |
---|
49 | |
---|
50 | // Initialization of the |
---|
51 | G4bool G4QANuENuclearCrossSection::onlyCS=true;//Flag to calculate only CS (not QE) |
---|
52 | G4double G4QANuENuclearCrossSection::lastSig=0.;//Last calculated total cross section |
---|
53 | G4double G4QANuENuclearCrossSection::lastQEL=0.;//Last calculated quasi-el. cross section |
---|
54 | G4int G4QANuENuclearCrossSection::lastL=0; //Last used in cross section TheLastBin |
---|
55 | G4double G4QANuENuclearCrossSection::lastE=0.; //Last used in cross section TheEnergy |
---|
56 | G4double* G4QANuENuclearCrossSection::lastEN=0; //Pointer to the Energy Scale of TX & QE |
---|
57 | G4double* G4QANuENuclearCrossSection::lastTX=0; //Pointer to the LastArray of TX function |
---|
58 | G4double* G4QANuENuclearCrossSection::lastQE=0; //Pointer to the LastArray of QE function |
---|
59 | G4int G4QANuENuclearCrossSection::lastPDG=0; // The last PDG code of the projectile |
---|
60 | G4int G4QANuENuclearCrossSection::lastN=0; // The last N of calculated nucleus |
---|
61 | G4int G4QANuENuclearCrossSection::lastZ=0; // The last Z of calculated nucleus |
---|
62 | G4double G4QANuENuclearCrossSection::lastP=0.; // Last used in cross section Momentum |
---|
63 | G4double G4QANuENuclearCrossSection::lastTH=0.; // Last threshold momentum |
---|
64 | G4double G4QANuENuclearCrossSection::lastCS=0.; // Last value of the Cross Section |
---|
65 | G4int G4QANuENuclearCrossSection::lastI=0; // The last position in the DAMDB |
---|
66 | std::vector<G4double*>* G4QANuENuclearCrossSection::TX = new std::vector<G4double*>; |
---|
67 | std::vector<G4double*>* G4QANuENuclearCrossSection::QE = new std::vector<G4double*>; |
---|
68 | |
---|
69 | // Returns Pointer to the G4VQCrossSection class |
---|
70 | G4VQCrossSection* G4QANuENuclearCrossSection::GetPointer() |
---|
71 | { |
---|
72 | static G4QANuENuclearCrossSection theCrossSection;//**Static body of the Cross Section** |
---|
73 | return &theCrossSection; |
---|
74 | } |
---|
75 | |
---|
76 | G4QANuENuclearCrossSection::~G4QANuENuclearCrossSection() |
---|
77 | { |
---|
78 | G4int lens=TX->size(); |
---|
79 | for(G4int i=0; i<lens; ++i) delete[] (*TX)[i]; |
---|
80 | delete TX; |
---|
81 | G4int hens=QE->size(); |
---|
82 | for(G4int i=0; i<hens; ++i) delete[] (*QE)[i]; |
---|
83 | delete QE; |
---|
84 | } |
---|
85 | |
---|
86 | // The main member function giving the collision cross section (P is in IU, CS is in mb) |
---|
87 | // Make pMom in independent units ! (Now it is MeV) |
---|
88 | G4double G4QANuENuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom, |
---|
89 | G4int tgZ, G4int tgN, G4int pPDG) |
---|
90 | { |
---|
91 | static G4int j; // A#0f records found in DB for this projectile |
---|
92 | static std::vector <G4int> colPDG;// Vector of the projectile PDG code |
---|
93 | static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops) |
---|
94 | static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops) |
---|
95 | static std::vector <G4double> colP; // Vector of last momenta for the reaction |
---|
96 | static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction |
---|
97 | static std::vector <G4double> colCS; // Vector of last cross sections for the reaction |
---|
98 | // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** |
---|
99 | G4double pEn=pMom; |
---|
100 | #ifdef debug |
---|
101 | G4cout<<"G4QAENCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN |
---|
102 | <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz=" |
---|
103 | <<colN.size()<<G4endl; |
---|
104 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
105 | #endif |
---|
106 | if(pPDG!=-12) |
---|
107 | { |
---|
108 | #ifdef debug |
---|
109 | G4cout<<"G4QAENCS::GetCS: *** Found pPDG="<<pPDG<<" ====> CS=0"<<G4endl; |
---|
110 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
111 | #endif |
---|
112 | return 0.; // projectile PDG=0 is a mistake (?!) @@ |
---|
113 | } |
---|
114 | G4bool in=false; // By default the isotope must be found in the AMDB |
---|
115 | if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope |
---|
116 | { |
---|
117 | in = false; // By default the isotope haven't be found in AMDB |
---|
118 | lastP = 0.; // New momentum history (nothing to compare with) |
---|
119 | lastPDG = pPDG; // The last PDG of the projectile |
---|
120 | lastN = tgN; // The last N of the calculated nucleus |
---|
121 | lastZ = tgZ; // The last Z of the calculated nucleus |
---|
122 | lastI = colN.size(); // Size of the Associative Memory DB in the heap |
---|
123 | j = 0; // A#0f records found in DB for this projectile |
---|
124 | if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found |
---|
125 | { // The nucleus with projPDG is found in AMDB |
---|
126 | if(colN[i]==tgN && colZ[i]==tgZ) |
---|
127 | { |
---|
128 | lastI=i; |
---|
129 | lastTH =colTH[i]; // Last THreshold (A-dependent) |
---|
130 | #ifdef pdebug |
---|
131 | G4cout<<"G4QAENCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl; |
---|
132 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
133 | #endif |
---|
134 | if(pEn<=lastTH) |
---|
135 | { |
---|
136 | #ifdef pdebug |
---|
137 | G4cout<<"G4QAENCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl; |
---|
138 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
139 | #endif |
---|
140 | return 0.; // Energy is below the Threshold value |
---|
141 | } |
---|
142 | lastP =colP [i]; // Last Momentum (A-dependent) |
---|
143 | lastCS =colCS[i]; // Last CrossSect (A-dependent) |
---|
144 | if(std::fabs(lastP/pMom-1.)<tolerance) |
---|
145 | { |
---|
146 | #ifdef pdebug |
---|
147 | G4cout<<"G4QAENCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; |
---|
148 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
149 | #endif |
---|
150 | return lastCS*millibarn; // Use theLastCS |
---|
151 | } |
---|
152 | in = true; // This is the case when the isotop is found in DB |
---|
153 | // Momentum pMom is in IU ! @@ Units |
---|
154 | #ifdef pdebug |
---|
155 | G4cout<<"G4QAENCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl; |
---|
156 | #endif |
---|
157 | lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update |
---|
158 | #ifdef pdebug |
---|
159 | G4cout<<"G4QAENCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl; |
---|
160 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
161 | #endif |
---|
162 | if(lastCS<=0. && pEn>lastTH) // Correct the threshold |
---|
163 | { |
---|
164 | #ifdef pdebug |
---|
165 | G4cout<<"G4QAENCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; |
---|
166 | #endif |
---|
167 | lastTH=pEn; |
---|
168 | } |
---|
169 | break; // Go out of the LOOP |
---|
170 | } |
---|
171 | #ifdef pdebug |
---|
172 | G4cout<<"---G4QAENCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i] |
---|
173 | <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl; |
---|
174 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
175 | #endif |
---|
176 | j++; // Increment a#0f records found in DB for this pPDG |
---|
177 | } |
---|
178 | if(!in) // This nucleus has not been calculated previously |
---|
179 | { |
---|
180 | #ifdef pdebug |
---|
181 | G4cout<<"G4QAENCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<lastI<<G4endl; |
---|
182 | #endif |
---|
183 | //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU) |
---|
184 | lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create |
---|
185 | if(lastCS<=0.) |
---|
186 | { |
---|
187 | lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last |
---|
188 | #ifdef pdebug |
---|
189 | G4cout<<"G4QAENCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl; |
---|
190 | #endif |
---|
191 | if(pEn>lastTH) |
---|
192 | { |
---|
193 | #ifdef pdebug |
---|
194 | G4cout<<"G4QAENCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl; |
---|
195 | #endif |
---|
196 | lastTH=pEn; |
---|
197 | } |
---|
198 | } |
---|
199 | #ifdef pdebug |
---|
200 | G4cout<<"G4QAENCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl; |
---|
201 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
202 | #endif |
---|
203 | colN.push_back(tgN); |
---|
204 | colZ.push_back(tgZ); |
---|
205 | colPDG.push_back(pPDG); |
---|
206 | colP.push_back(pMom); |
---|
207 | colTH.push_back(lastTH); |
---|
208 | colCS.push_back(lastCS); |
---|
209 | #ifdef pdebug |
---|
210 | G4cout<<"G4QAENCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl; |
---|
211 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
212 | #endif |
---|
213 | return lastCS*millibarn; |
---|
214 | } // End of creation of the new set of parameters |
---|
215 | else |
---|
216 | { |
---|
217 | #ifdef pdebug |
---|
218 | G4cout<<"G4QAENCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl; |
---|
219 | #endif |
---|
220 | colP[lastI]=pMom; |
---|
221 | colPDG[lastI]=pPDG; |
---|
222 | colCS[lastI]=lastCS; |
---|
223 | } |
---|
224 | } // End of parameters udate |
---|
225 | else if(pEn<=lastTH) |
---|
226 | { |
---|
227 | #ifdef pdebug |
---|
228 | G4cout<<"G4QAENCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl; |
---|
229 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
230 | #endif |
---|
231 | return 0.; // Momentum is below the Threshold Value -> CS=0 |
---|
232 | } |
---|
233 | else if(std::fabs(lastP/pMom-1.)<tolerance) |
---|
234 | { |
---|
235 | #ifdef pdebug |
---|
236 | G4cout<<"G4QAENCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl; |
---|
237 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
238 | #endif |
---|
239 | return lastCS*millibarn; // Use theLastCS |
---|
240 | } |
---|
241 | else |
---|
242 | { |
---|
243 | #ifdef pdebug |
---|
244 | G4cout<<"G4QAENCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl; |
---|
245 | #endif |
---|
246 | lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB |
---|
247 | lastP=pMom; |
---|
248 | } |
---|
249 | #ifdef pdebug |
---|
250 | G4cout<<"G4QAENCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl; |
---|
251 | //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST |
---|
252 | #endif |
---|
253 | return lastCS*millibarn; |
---|
254 | } |
---|
255 | |
---|
256 | // Gives the threshold energy = the same for all nuclei (@@ can be reduced for hevy nuclei) |
---|
257 | G4double G4QANuENuclearCrossSection::ThresholdEnergy(G4int Z, G4int N, G4int) |
---|
258 | { |
---|
259 | //static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/GeV; |
---|
260 | //static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/GeV; |
---|
261 | //static const G4double mDeut = G4NucleiProperties::GetNuclearMass(2,1)/GeV/2.; |
---|
262 | static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV) |
---|
263 | static const G4double dmN=mN+mN; // Doubled nucleon mass (2*AtomicMassUnit, GeV) |
---|
264 | static const G4double me=.00051099892; // electron mass in GeV |
---|
265 | static const G4double me2=me*me; // Squared mass of an electron in GeV^2 |
---|
266 | static const G4double thresh=me+me2/dmN; // Universal threshold in GeV |
---|
267 | // --------- |
---|
268 | //static const G4double infEn = 9.e27; |
---|
269 | G4double dN=0.; |
---|
270 | if(Z>0||N>0) dN=thresh*GeV; // @@ if upgraded, change it in a total cross section |
---|
271 | //@@ "dN=me+me2/G4NucleiProperties::GetNuclearMass(<G4double>(Z+N),<G4double>(Z)/GeV" |
---|
272 | return dN; |
---|
273 | } |
---|
274 | |
---|
275 | // The main member function giving the gamma-A cross section (E_kin in MeV, CS in mb) |
---|
276 | G4double G4QANuENuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I, |
---|
277 | G4int , G4int targZ, G4int targN, G4double Momentum) |
---|
278 | { |
---|
279 | static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2 |
---|
280 | static const G4int nE=65; // !! If change this, change it in GetFunctions() !! |
---|
281 | static const G4int mL=nE-1; |
---|
282 | static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV) |
---|
283 | static const G4double dmN=mN+mN; // Doubled nucleon mass (2*AtomicMassUnit, GeV) |
---|
284 | static const G4double me=.00051099892;// electron mass in GeV |
---|
285 | static const G4double me2=me*me; // Squared mass of an electron in GeV^2 |
---|
286 | static const G4double EMi=me+me2/dmN; // Universal threshold of the reaction in GeV |
---|
287 | static const G4double EMa=300.; // Maximum tabulated Energy of nu_e in GeV |
---|
288 | // *** Begin of the Associative memory for acceleration of the cross section calculations |
---|
289 | static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional) |
---|
290 | static G4bool first=true; // Flag of initialization of the energy axis |
---|
291 | // *** End of Static Definitions (Associative Memory) *** |
---|
292 | //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Electron |
---|
293 | //G4double TotEnergy2=Momentum; |
---|
294 | onlyCS=CS; // Flag to calculate only CS (not TX & QE) |
---|
295 | lastE=Momentum/GeV; // Kinetic energy of the electron neutrino (in GeV) |
---|
296 | if (lastE<=EMi) // Energy is below the minimum energy in the table |
---|
297 | { |
---|
298 | lastE=0.; |
---|
299 | lastSig=0.; |
---|
300 | return 0.; |
---|
301 | } |
---|
302 | G4int Z=targZ; // New Z, which can change the sign |
---|
303 | if(F<=0) // This isotope was not the last used isotop |
---|
304 | { |
---|
305 | if(F<0) // This isotope was found in DAMDB =========> RETRIEVE |
---|
306 | { |
---|
307 | lastTX =(*TX)[I]; // Pointer to the prepared TX function (same isotope) |
---|
308 | lastQE =(*QE)[I]; // Pointer to the prepared QE function (same isotope) |
---|
309 | } |
---|
310 | else // This isotope wasn't calculated previously => CREATE |
---|
311 | { |
---|
312 | if(first) |
---|
313 | { |
---|
314 | lastEN = new G4double[nE]; // This must be done only once! |
---|
315 | Z=-Z; // To explain GetFunctions that E-axis must be filled |
---|
316 | first=false; // To make it only once |
---|
317 | } |
---|
318 | lastTX = new G4double[nE]; // Allocate memory for the new TX function |
---|
319 | lastQE = new G4double[nE]; // Allocate memory for the new QE function |
---|
320 | G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK) |
---|
321 | if(res<0) G4cerr<<"*W*G4NuENuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl; |
---|
322 | // *** The synchronization check *** |
---|
323 | G4int sync=TX->size(); |
---|
324 | if(sync!=I) G4cerr<<"***G4NuENuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl; |
---|
325 | TX->push_back(lastTX); |
---|
326 | QE->push_back(lastQE); |
---|
327 | } // End of creation of the new set of parameters |
---|
328 | } // End of parameters udate |
---|
329 | // ============================== NOW Calculate the Cross Section ===================== |
---|
330 | if (lastE<=EMi) // Check that antiNuEnergy is higher than ThreshE |
---|
331 | { |
---|
332 | lastE=0.; |
---|
333 | lastSig=0.; |
---|
334 | return 0.; |
---|
335 | } |
---|
336 | if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization |
---|
337 | { |
---|
338 | G4int chk=1; |
---|
339 | G4int ran=mL/2; |
---|
340 | G4int sep=ran; // as a result = an index of the left edge of the interval |
---|
341 | while(ran>=2) |
---|
342 | { |
---|
343 | G4int newran=ran/2; |
---|
344 | G4double oldL=lastEN[sep]; |
---|
345 | if(lastE<=oldL) sep-=newran; |
---|
346 | else sep+=newran; |
---|
347 | #ifdef pdebug |
---|
348 | G4cout<<"G4ANuE::CCS:n="<<newran<<",s="<<sep<<",E="<<lastE<<",oE="<<oldL<<G4endl; |
---|
349 | #endif |
---|
350 | ran=newran; |
---|
351 | chk=chk+chk; |
---|
352 | } |
---|
353 | if(chk+chk!=mL) G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Table! mL="<<mL<<G4endl; |
---|
354 | G4double lowE=lastEN[sep]; |
---|
355 | G4double highE=lastEN[sep+1]; |
---|
356 | G4double lowTX=lastTX[sep]; |
---|
357 | if(lastE<lowE||sep>=mL||lastE>highE) |
---|
358 | G4cerr<<"*Warn*G4ANuENuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE |
---|
359 | <<", sep="<<sep<<", mL="<<mL<<G4endl; |
---|
360 | lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E |
---|
361 | if(!onlyCS) // Skip the differential cross-section parameters |
---|
362 | { |
---|
363 | G4double lowQE=lastQE[sep]; |
---|
364 | lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE; |
---|
365 | #ifdef pdebug |
---|
366 | G4cout<<"G4ANuENuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl; |
---|
367 | #endif |
---|
368 | } |
---|
369 | } |
---|
370 | else |
---|
371 | { |
---|
372 | lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking... |
---|
373 | lastQEL=lastQE[mL]; |
---|
374 | } |
---|
375 | if(lastQEL<0.) lastQEL = 0.; |
---|
376 | if(lastSig<0.) lastSig = 0.; |
---|
377 | // The cross-sections are expected to be in mb |
---|
378 | lastSig*=mb38; |
---|
379 | if(!onlyCS) lastQEL*=mb38; |
---|
380 | return lastSig; |
---|
381 | } |
---|
382 | |
---|
383 | // Calculate the cros-section functions |
---|
384 | // **************************************************************************************** |
---|
385 | // *** This tables are the same for all lepto-nuclear reactions, only mass is different *** |
---|
386 | // ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ *** |
---|
387 | // **************************************************************************************** |
---|
388 | G4int G4QANuENuclearCrossSection::GetFunctions (G4int z, G4int n, |
---|
389 | G4double* t, G4double* q, G4double* e) |
---|
390 | { |
---|
391 | static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV) |
---|
392 | static const G4double dmN=mN+mN; // Doubled nucleon mass (2*AtomicMassUnit, GeV) |
---|
393 | static const G4double me=.00051099892; // electron mass in GeV |
---|
394 | static const G4double me2=me*me; // Squared mass of an electron in GeV^2 |
---|
395 | static const G4double thresh=me+me2/dmN; // Universal threshold in GeV |
---|
396 | static const G4int nE=65; // !! If change this, change it in CalculateCrossSection() !! |
---|
397 | static const G4double nuEn[nE]={thresh, |
---|
398 | .00051331,.00053602,.00056078,.00058783,.00061743,.00064990,.00068559,.00072492, |
---|
399 | .00076834,.00081641,.00086975,.00092912,.00099536,.00106950,.00115273,.00124646, |
---|
400 | .00135235,.00147241,.00160901,.00176503,.00194392,.00214986,.00238797,.00266448, |
---|
401 | .00298709,.00336531,.00381094,.00433879,.00496745,.00572047,.00662785,.00772806, |
---|
402 | .00907075,.01072050,.01276190,.01530660,.01850330,.02255110,.02771990,.03437780, |
---|
403 | .04303240,.05438970,.06944210,.08959920,.11688400,.15423600,.20597200,.27851200, |
---|
404 | .38153100,.52979600,.74616300,1.0665200,1.5480900,2.2834800,3.4251100,5.2281000, |
---|
405 | 8.1270200,12.875900,20.808500,34.331200,57.877800,99.796200,176.16300,318.68200}; |
---|
406 | static const G4double TOTX[nE]={0., |
---|
407 | .00046923,.00160693,.00229560,.00288772,.00344344,.00398831,.00453725,.00510087, |
---|
408 | .00568797,.00630663,.00696489,.00767121,.00843477,.00926586,.01017610,.01117910, |
---|
409 | .01229040,.01352820,.01491430,.01647420,.01823820,.02024280,.02253130,.02515600, |
---|
410 | .02818010,.03167950,.03574640,.04049260,.04605330,.05739330,.07028190,.08396290, |
---|
411 | .09969400,.11756200,.13812500,.16207000,.19099600,.22447700,.26434900,.31051900, |
---|
412 | .36357400,.42242900,.48500000,.54692800,.60140500,.63996800,.65519300,.64339200, |
---|
413 | .60784300,.55734400,.50212600,.45116800,.40962600,.37854800,.35702100,.34330900, |
---|
414 | .33576200,.33300100,.33387100,.33737700,.34235200,.34880100,.35507500,.35961200}; |
---|
415 | static const G4double QELX[nE]={0., |
---|
416 | 2.40856e-7,8.61338e-7,1.28732e-6,1.69747e-6,2.12608e-6,2.59201e-6,3.11071e-6,3.69770e-6, |
---|
417 | 4.37028e-6,5.14877e-6,6.05773e-6,7.12746e-6,8.39567e-6,9.90987e-6,1.17304e-5,1.39343e-5, |
---|
418 | 1.66209e-5,1.99191e-5,2.39974e-5,2.90775e-5,3.54536e-5,4.35191e-5,5.38039e-5,6.70278e-5, |
---|
419 | 8.41765e-5,.000106611,.000136228,.000175689,.000228768,.000328317,.000465818,.000648870, |
---|
420 | .000904299,.001260330,.001762740,.002480740,.003534040,.005062210,.007327720,.010675000, |
---|
421 | .015645500,.022975800,.033679400,.049004300,.070294800,.098706100,.134951000,.179193000, |
---|
422 | .231297000,.287287000,.344760000,.404397000,.465915000,.527060000,.584085000,.632748000, |
---|
423 | .671346000,.699744000,.719355000,.732541000,.740996000,.746249000,.749265000,.751160000}; |
---|
424 | // -------------------------------- |
---|
425 | G4int first=0; |
---|
426 | if(z<0.) |
---|
427 | { |
---|
428 | first=1; |
---|
429 | z=-z; |
---|
430 | } |
---|
431 | if(z<1 || z>92) // neutron & plutonium are forbidden |
---|
432 | { |
---|
433 | G4cout<<"***G4QANuENuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl; |
---|
434 | return -1; |
---|
435 | } |
---|
436 | for(G4int k=0; k<nE; k++) |
---|
437 | { |
---|
438 | G4double a=n+z; |
---|
439 | G4double za=z+a; |
---|
440 | G4double dz=z+z; |
---|
441 | G4double da=a+a; |
---|
442 | G4double ta=da+a; |
---|
443 | if(first) e[k]=nuEn[k]; // Energy of neutrino E (first bin k=0 can be modified) |
---|
444 | t[k]=TOTX[k]*nuEn[k]*(za+za)/ta+QELX[k]*(dz+dz-da)/ta; // TotalCrossSection |
---|
445 | q[k]=QELX[k]*dz/a; // QuasiElasticCrossSection |
---|
446 | } |
---|
447 | return first; |
---|
448 | } |
---|
449 | |
---|
450 | // Randomize Q2 from neutrino to the scattered electron when scattering is quasi-elastic |
---|
451 | G4double G4QANuENuclearCrossSection::GetQEL_ExchangeQ2() |
---|
452 | { |
---|
453 | static const G4double me=.00051099892; // electron mass in GeV |
---|
454 | static const G4double me2=me*me; // Squared mass of an electron in GeV^2 |
---|
455 | static const G4double hme2=me2/2; // .5*m_e^2 in GeV^2 |
---|
456 | static const double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV) |
---|
457 | static const double MN2=MN*MN; // M_N^2 in GeV^2 |
---|
458 | static const G4double power=-3.5; // direct power for the magic variable |
---|
459 | static const G4double pconv=1./power;// conversion power for the magic variable |
---|
460 | static const G4int nQ2=101; // #Of point in the Q2l table (in GeV^2) |
---|
461 | static const G4int lQ2=nQ2-1; // index of the last in the Q2l table |
---|
462 | static const G4int bQ2=lQ2-1; // index of the before last in the Q2 ltable |
---|
463 | // Reversed table |
---|
464 | static const G4double Xl[nQ2]={5.20224e-16, |
---|
465 | .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898, |
---|
466 | .093878, .103332, .112844, .122410, .132023, .141680, .151376, .161109, .170875, .180672, |
---|
467 | .190499, .200352, .210230, .220131, .230055, .239999, .249963, .259945, .269944, .279960, |
---|
468 | .289992, .300039, .310099, .320173, .330260, .340359, .350470, .360592, .370724, .380867, |
---|
469 | .391019, .401181, .411352, .421531, .431719, .441915, .452118, .462329, .472547, .482771, |
---|
470 | .493003, .503240, .513484, .523734, .533989, .544250, .554517, .564788, .575065, .585346, |
---|
471 | .595632, .605923, .616218, .626517, .636820, .647127, .657438, .667753, .678072, .688394, |
---|
472 | .698719, .709048, .719380, .729715, .740053, .750394, .760738, .771085, .781434, .791786, |
---|
473 | .802140, .812497, .822857, .833219, .843582, .853949, .864317, .874687, .885060, .895434, |
---|
474 | .905810, .916188, .926568, .936950, .947333, .957719, .968105, .978493, .988883, .999275}; |
---|
475 | // Direct table |
---|
476 | static const G4double Xmax=Xl[lQ2]; |
---|
477 | static const G4double Xmin=Xl[0]; |
---|
478 | static const G4double dX=(Xmax-Xmin)/lQ2; // step in X(Q2, GeV^2) |
---|
479 | static const G4double inl[nQ2]={0, |
---|
480 | 1.52225, 2.77846, 3.96651, 5.11612, 6.23990, 7.34467, 8.43466, 9.51272, 10.5809, 11.6406, |
---|
481 | 12.6932, 13.7394, 14.7801, 15.8158, 16.8471, 17.8743, 18.8979, 19.9181, 20.9353, 21.9496, |
---|
482 | 22.9614, 23.9707, 24.9777, 25.9826, 26.9855, 27.9866, 28.9860, 29.9837, 30.9798, 31.9745, |
---|
483 | 32.9678, 33.9598, 34.9505, 35.9400, 36.9284, 37.9158, 38.9021, 39.8874, 40.8718, 41.8553, |
---|
484 | 42.8379, 43.8197, 44.8007, 45.7810, 46.7605, 47.7393, 48.7174, 49.6950, 50.6718, 51.6481, |
---|
485 | 52.6238, 53.5990, 54.5736, 55.5476, 56.5212, 57.4943, 58.4670, 59.4391, 60.4109, 61.3822, |
---|
486 | 62.3531, 63.3236, 64.2937, 65.2635, 66.2329, 67.2019, 68.1707, 69.1390, 70.1071, 71.0748, |
---|
487 | 72.0423, 73.0095, 73.9763, 74.9429, 75.9093, 76.8754, 77.8412, 78.8068, 79.7721, 80.7373, |
---|
488 | 81.7022, 82.6668, 83.6313, 84.5956, 85.5596, 86.5235, 87.4872, 88.4507, 89.4140, 90.3771, |
---|
489 | 91.3401, 92.3029, 93.2656, 94.2281, 95.1904, 96.1526, 97.1147, 98.0766, 99.0384, 100.000}; |
---|
490 | G4double Enu=lastE; // Get energy of the last calculated cross-section |
---|
491 | G4double dEnu=Enu+Enu; // doubled energy of nu/anu |
---|
492 | G4double Enu2=Enu*Enu; // squared energy of nu/anu |
---|
493 | G4double ME=Enu*MN; // M*E |
---|
494 | G4double dME=ME+ME; // 2*M*E |
---|
495 | G4double dEMN=(dEnu+MN)*ME; |
---|
496 | G4double MEm=ME-hme2; |
---|
497 | G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2); |
---|
498 | G4double E2M=MN*Enu2-(Enu+MN)*hme2; |
---|
499 | G4double ymax=(E2M+sqE)/dEMN; |
---|
500 | G4double ymin=(E2M-sqE)/dEMN; |
---|
501 | G4double rmin=1.-ymin; |
---|
502 | G4double rhm2E=hme2/Enu2; |
---|
503 | G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu) |
---|
504 | G4double Q2ma=dME*ymax; // Q2_max(E_nu) |
---|
505 | G4double Xma=std::pow((1.+Q2mi),power); // X_max(E_nu) |
---|
506 | G4double Xmi=std::pow((1.+Q2ma),power); // X_min(E_nu) |
---|
507 | // Find the integral values integ(Xmi) & integ(Xma) using the direct table |
---|
508 | G4double rXi=(Xmi-Xmin)/dX; |
---|
509 | G4int iXi=static_cast<int>(rXi); |
---|
510 | if(iXi<0) iXi=0; |
---|
511 | if(iXi>bQ2) iXi=bQ2; |
---|
512 | G4double dXi=rXi-iXi; |
---|
513 | G4double bnti=inl[iXi]; |
---|
514 | G4double inti=bnti+dXi*(inl[iXi+1]-bnti); |
---|
515 | // |
---|
516 | G4double rXa=(Xma-Xmin)/dX; |
---|
517 | G4int iXa=static_cast<int>(rXa); |
---|
518 | if(iXa<0) iXa=0; |
---|
519 | if(iXa>bQ2) iXa=bQ2; |
---|
520 | G4double dXa=rXa-iXa; |
---|
521 | G4double bnta=inl[iXa]; |
---|
522 | G4double inta=bnta+dXa*(inl[iXa+1]-bnta); |
---|
523 | // *** Find X using the reversed table *** |
---|
524 | G4double intx=inti+(inta-inti)*G4UniformRand(); |
---|
525 | G4int intc=static_cast<int>(intx); |
---|
526 | if(intc<0) intc=0; |
---|
527 | if(intc>bQ2) intc=bQ2; // If it is more than max, then the BAD extrapolation |
---|
528 | G4double dint=intx-intc; |
---|
529 | G4double mX=Xl[intc]; |
---|
530 | G4double X=mX+dint*(Xl[intc+1]-mX); |
---|
531 | G4double Q2=std::pow(X,pconv)-1.; |
---|
532 | return Q2*GeV*GeV; |
---|
533 | } |
---|
534 | |
---|
535 | // Randomize Q2 from neutrino to the scattered electron when scattering is not quasiElastic |
---|
536 | G4double G4QANuENuclearCrossSection::GetNQE_ExchangeQ2() |
---|
537 | { |
---|
538 | static const double mpi=.13957018; // charged pi meson mass in GeV |
---|
539 | static const G4double me=.00051099892;// electron mass in GeV |
---|
540 | static const G4double me2=me*me; // Squared mass of an electron in GeV^2 |
---|
541 | static const G4double hme2=me2/2; // .5*m_e^2 in GeV^2 |
---|
542 | static const double MN=.931494043; // Nucleon mass (inside nucleus,atomicMassUnit,GeV) |
---|
543 | static const double MN2=MN*MN; // M_N^2 in GeV^2 |
---|
544 | static const double dMN=MN+MN; // 2*M_N in GeV |
---|
545 | static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic |
---|
546 | static const G4int power=7; // direct power for the magic variable |
---|
547 | static const G4double pconv=1./power; // conversion power for the magic variable |
---|
548 | static const G4int nX=21; // #Of point in the Xl table (in GeV^2) |
---|
549 | static const G4int lX=nX-1; // index of the last in the Xl table |
---|
550 | static const G4int bX=lX-1; // @@ index of the before last in the Xl table |
---|
551 | static const G4int nE=20; // #Of point in the El table (in GeV^2) |
---|
552 | static const G4int bE=nE-1; // index of the last in the El table |
---|
553 | static const G4int pE=bE-1; // index of the before last in the El table |
---|
554 | // Reversed table |
---|
555 | static const G4double X0[nX]={5.21412e-05, |
---|
556 | .437860, .681908, .891529, 1.08434, 1.26751, 1.44494, 1.61915, 1.79198, 1.96493, 2.13937, |
---|
557 | 2.31664, 2.49816, 2.68559, 2.88097, 3.08705, 3.30774, 3.54917, 3.82233, 4.15131, 4.62182}; |
---|
558 | static const G4double X1[nX]={.00102591, |
---|
559 | 1.00443, 1.55828, 2.03126, 2.46406, 2.87311, 3.26723, 3.65199, 4.03134, 4.40835, 4.78561, |
---|
560 | 5.16549, 5.55031, 5.94252, 6.34484, 6.76049, 7.19349, 7.64917, 8.13502, 8.66246, 9.25086}; |
---|
561 | static const G4double X2[nX]={.0120304, |
---|
562 | 2.59903, 3.98637, 5.15131, 6.20159, 7.18024, 8.10986, 9.00426, 9.87265, 10.7217, 11.5564, |
---|
563 | 12.3808, 13.1983, 14.0116, 14.8234, 15.6359, 16.4515, 17.2723, 18.1006, 18.9386, 19.7892}; |
---|
564 | static const G4double X3[nX]={.060124, |
---|
565 | 5.73857, 8.62595, 10.9849, 13.0644, 14.9636, 16.7340, 18.4066, 20.0019, 21.5342, 23.0142, |
---|
566 | 24.4497, 25.8471, 27.2114, 28.5467, 29.8564, 31.1434, 32.4102, 33.6589, 34.8912, 36.1095}; |
---|
567 | static const G4double X4[nX]={.0992363, |
---|
568 | 8.23746, 12.1036, 15.1740, 17.8231, 20.1992, 22.3792, 24.4092, 26.3198, 28.1320, 29.8615, |
---|
569 | 31.5200, 33.1169, 34.6594, 36.1536, 37.6044, 39.0160, 40.3920, 41.7353, 43.0485, 44.3354}; |
---|
570 | static const G4double X5[nX]={.0561127, |
---|
571 | 7.33661, 10.5694, 13.0778, 15.2061, 17.0893, 18.7973, 20.3717, 21.8400, 23.2211, 24.5291, |
---|
572 | 25.7745, 26.9655, 28.1087, 29.2094, 30.2721, 31.3003, 32.2972, 33.2656, 34.2076, 35.1265}; |
---|
573 | static const G4double X6[nX]={.0145859, |
---|
574 | 4.81774, 6.83565, 8.37399, 9.66291, 10.7920, 11.8075, 12.7366, 13.5975, 14.4025, 15.1608, |
---|
575 | 15.8791, 16.5628, 17.2162, 17.8427, 18.4451, 19.0259, 19.5869, 20.1300, 20.6566, 21.1706}; |
---|
576 | static const G4double X7[nX]={.00241155, |
---|
577 | 2.87095, 4.02492, 4.89243, 5.61207, 6.23747, 6.79613, 7.30433, 7.77270, 8.20858, 8.61732, |
---|
578 | 9.00296, 9.36863, 9.71682, 10.0495, 10.3684, 10.6749, 10.9701, 11.2550, 11.5306, 11.7982}; |
---|
579 | static const G4double X8[nX]={.000316863, |
---|
580 | 1.76189, 2.44632, 2.95477, 3.37292, 3.73378, 4.05420, 4.34415, 4.61009, 4.85651, 5.08666, |
---|
581 | 5.30299, 5.50738, 5.70134, 5.88609, 6.06262, 6.23178, 6.39425, 6.55065, 6.70149, 6.84742}; |
---|
582 | static const G4double X9[nX]={3.73544e-05, |
---|
583 | 1.17106, 1.61289, 1.93763, 2.20259, 2.42976, 2.63034, 2.81094, 2.97582, 3.12796, 3.26949, |
---|
584 | 3.40202, 3.52680, 3.64482, 3.75687, 3.86360, 3.96557, 4.06323, 4.15697, 4.24713, 4.33413}; |
---|
585 | static const G4double XA[nX]={4.19131e-06, |
---|
586 | .849573, 1.16208, 1.38955, 1.57379, 1.73079, 1.86867, 1.99221, 2.10451, 2.20770, 2.30332, |
---|
587 | 2.39252, 2.47622, 2.55511, 2.62977, 2.70066, 2.76818, 2.83265, 2.89437, 2.95355, 3.01051}; |
---|
588 | static const G4double XB[nX]={4.59981e-07, |
---|
589 | .666131, .905836, 1.07880, 1.21796, 1.33587, 1.43890, 1.53080, 1.61399, 1.69011, 1.76040, |
---|
590 | 1.82573, 1.88682, 1.94421, 1.99834, 2.04959, 2.09824, 2.14457, 2.18878, 2.23107, 2.27162}; |
---|
591 | static const G4double XC[nX]={4.99861e-08, |
---|
592 | .556280, .752730, .893387, 1.00587, 1.10070, 1.18317, 1.25643, 1.32247, 1.38269, 1.43809, |
---|
593 | 1.48941, 1.53724, 1.58203, 1.62416, 1.66391, 1.70155, 1.73728, 1.77128, 1.80371, 1.83473}; |
---|
594 | static const G4double XD[nX]={5.40832e-09, |
---|
595 | .488069, .657650, .778236, .874148, .954621, 1.02432, 1.08599, 1.14138, 1.19172, 1.23787, |
---|
596 | 1.28049, 1.32008, 1.35705, 1.39172, 1.42434, 1.45514, 1.48429, 1.51197, 1.53829, 1.56339}; |
---|
597 | static const G4double XE[nX]={5.84029e-10, |
---|
598 | .445057, .597434, .705099, .790298, .861468, .922865, .976982, 1.02542, 1.06930, 1.10939, |
---|
599 | 1.14630, 1.18050, 1.21233, 1.24208, 1.27001, 1.29630, 1.32113, 1.34462, 1.36691, 1.38812}; |
---|
600 | static const G4double XF[nX]={6.30137e-11, |
---|
601 | .418735, .560003, .659168, .737230, .802138, .857898, .906854, .950515, .989915, 1.02580, |
---|
602 | 1.05873, 1.08913, 1.11734, 1.14364, 1.16824, 1.19133, 1.21306, 1.23358, 1.25298, 1.27139}; |
---|
603 | static const G4double XG[nX]={6.79627e-12, |
---|
604 | .405286, .539651, .633227, .706417, .766929, .818642, .863824, .903931, .939963, .972639, |
---|
605 | 1.00250, 1.02995, 1.05532, 1.07887, 1.10082, 1.12134, 1.14058, 1.15867, 1.17572, 1.19183}; |
---|
606 | static const G4double XH[nX]={7.32882e-13, |
---|
607 | .404391, .535199, .625259, .695036, .752243, .800752, .842823, .879906, .912994, .942802, |
---|
608 | .969862, .994583, 1.01729, 1.03823, 1.05763, 1.07566, 1.09246, 1.10816, 1.12286, 1.13667}; |
---|
609 | static const G4double XI[nX]={7.90251e-14, |
---|
610 | .418084, .548382, .636489, .703728, .758106, .803630, .842633, .876608, .906576, .933269, |
---|
611 | .957233, .978886, .998556, 1.01651, 1.03295, 1.04807, 1.06201, 1.07489, 1.08683, 1.09792}; |
---|
612 | static const G4double XJ[nX]={8.52083e-15, |
---|
613 | .447299, .579635, .666780, .731788, .783268, .825512, .861013, .891356, .917626, .940597, |
---|
614 | .960842, .978802, .994820, 1.00917, 1.02208, 1.03373, 1.04427, 1.05383, 1.06253, 1.07046}; |
---|
615 | // Direct table |
---|
616 | static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0], |
---|
617 | X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]}; |
---|
618 | static const G4double dX[nE]={ |
---|
619 | (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX, |
---|
620 | (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX, |
---|
621 | (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX, |
---|
622 | (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX, |
---|
623 | (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX}; |
---|
624 | static const G4double* Xl[nE]= |
---|
625 | {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ}; |
---|
626 | static const G4double I0[nX]={0, |
---|
627 | .354631, 1.08972, 2.05138, 3.16564, 4.38343, 5.66828, 6.99127, 8.32858, 9.65998, 10.9680, |
---|
628 | 12.2371, 13.4536, 14.6050, 15.6802, 16.6686, 17.5609, 18.3482, 19.0221, 19.5752, 20.0000}; |
---|
629 | static const G4double I1[nX]={0, |
---|
630 | .281625, .877354, 1.67084, 2.60566, 3.64420, 4.75838, 5.92589, 7.12829, 8.34989, 9.57708, |
---|
631 | 10.7978, 12.0014, 13.1781, 14.3190, 15.4162, 16.4620, 17.4496, 18.3724, 19.2245, 20.0000}; |
---|
632 | static const G4double I2[nX]={0, |
---|
633 | .201909, .642991, 1.24946, 1.98463, 2.82370, 3.74802, 4.74263, 5.79509, 6.89474, 8.03228, |
---|
634 | 9.19947, 10.3889, 11.5938, 12.8082, 14.0262, 15.2427, 16.4527, 17.6518, 18.8356, 20.0000}; |
---|
635 | static const G4double I3[nX]={0, |
---|
636 | .140937, .461189, .920216, 1.49706, 2.17728, 2.94985, 3.80580, 4.73758, 5.73867, 6.80331, |
---|
637 | 7.92637, 9.10316, 10.3294, 11.6013, 12.9150, 14.2672, 15.6548, 17.0746, 18.5239, 20.0000}; |
---|
638 | static const G4double I4[nX]={0, |
---|
639 | .099161, .337358, .694560, 1.16037, 1.72761, 2.39078, 3.14540, 3.98768, 4.91433, 5.92245, |
---|
640 | 7.00942, 8.17287, 9.41060, 10.7206, 12.1010, 13.5500, 15.0659, 16.6472, 18.2924, 20.0000}; |
---|
641 | static const G4double I5[nX]={0, |
---|
642 | .071131, .255084, .543312, .932025, 1.41892, 2.00243, 2.68144, 3.45512, 4.32283, 5.28411, |
---|
643 | 6.33859, 7.48602, 8.72621, 10.0590, 11.4844, 13.0023, 14.6128, 16.3158, 18.1115, 20.0000}; |
---|
644 | static const G4double I6[nX]={0, |
---|
645 | .053692, .202354, .443946, .778765, 1.20774, 1.73208, 2.35319, 3.07256, 3.89177, 4.81249, |
---|
646 | 5.83641, 6.96528, 8.20092, 9.54516, 10.9999, 12.5670, 14.2486, 16.0466, 17.9630, 20.0000}; |
---|
647 | static const G4double I7[nX]={0, |
---|
648 | .043065, .168099, .376879, .672273, 1.05738, 1.53543, 2.10973, 2.78364, 3.56065, 4.44429, |
---|
649 | 5.43819, 6.54610, 7.77186, 9.11940, 10.5928, 12.1963, 13.9342, 15.8110, 17.8313, 20.0000}; |
---|
650 | static const G4double I8[nX]={0, |
---|
651 | .036051, .143997, .327877, .592202, .941572, 1.38068, 1.91433, 2.54746, 3.28517, 4.13277, |
---|
652 | 5.09574, 6.17984, 7.39106, 8.73568, 10.2203, 11.8519, 13.6377, 15.5854, 17.7033, 20.0000}; |
---|
653 | static const G4double I9[nX]={0, |
---|
654 | .030977, .125727, .289605, .528146, .846967, 1.25183, 1.74871, 2.34384, 3.04376, 3.85535, |
---|
655 | 4.78594, 5.84329, 7.03567, 8.37194, 9.86163, 11.5150, 13.3430, 15.3576, 17.5719, 20.0000}; |
---|
656 | static const G4double IA[nX]={0, |
---|
657 | .027129, .111420, .258935, .475812, .768320, 1.14297, 1.60661, 2.16648, 2.83034, 3.60650, |
---|
658 | 4.50394, 5.53238, 6.70244, 8.02569, 9.51488, 11.1841, 13.0488, 15.1264, 17.4362, 20.0000}; |
---|
659 | static const G4double IB[nX]={0, |
---|
660 | .024170, .100153, .234345, .433198, .703363, 1.05184, 1.48607, 2.01409, 2.64459, 3.38708, |
---|
661 | 4.25198, 5.25084, 6.39647, 7.70319, 9.18708, 10.8663, 12.7617, 14.8968, 17.2990, 20.0000}; |
---|
662 | static const G4double IC[nX]={0, |
---|
663 | .021877, .091263, .214670, .398677, .650133, .976322, 1.38510, 1.88504, 2.48555, 3.19709, |
---|
664 | 4.03129, 5.00127, 6.12184, 7.40989, 8.88482, 10.5690, 12.4888, 14.6748, 17.1638, 20.0000}; |
---|
665 | static const G4double ID[nX]={0, |
---|
666 | .020062, .084127, .198702, .370384, .606100, .913288, 1.30006, 1.77535, 2.34912, 3.03253, |
---|
667 | 3.83822, 4.78063, 5.87634, 7.14459, 8.60791, 10.2929, 12.2315, 14.4621, 17.0320, 20.0000}; |
---|
668 | static const G4double IE[nX]={0, |
---|
669 | .018547, .078104, .185102, .346090, .567998, .858331, 1.22535, 1.67824, 2.22735, 2.88443, |
---|
670 | 3.66294, 4.57845, 5.64911, 6.89637, 8.34578, 10.0282, 11.9812, 14.2519, 16.8993, 20.0000}; |
---|
671 | static const G4double IF[nX]={0, |
---|
672 | .017143, .072466, .172271, .323007, .531545, .805393, 1.15288, 1.58338, 2.10754, 2.73758, |
---|
673 | 3.48769, 4.37450, 5.41770, 6.64092, 8.07288, 9.74894, 11.7135, 14.0232, 16.7522, 20.0000}; |
---|
674 | static const G4double IG[nX]={0, |
---|
675 | .015618, .066285, .158094, .297316, .490692, .745653, 1.07053, 1.47479, 1.96931, 2.56677, |
---|
676 | 3.28205, 4.13289, 5.14068, 6.33158, 7.73808, 9.40133, 11.3745, 13.7279, 16.5577, 20.0000}; |
---|
677 | static const G4double IH[nX]={0, |
---|
678 | .013702, .058434, .139923, .264115, .437466, .667179, .961433, 1.32965, 1.78283, 2.33399, |
---|
679 | 2.99871, 3.79596, 4.74916, 5.88771, 7.24937, 8.88367, 10.8576, 13.2646, 16.2417, 20.0000}; |
---|
680 | static const G4double II[nX]={0, |
---|
681 | .011264, .048311, .116235, .220381, .366634, .561656, .813132, 1.13008, 1.52322, 2.00554, |
---|
682 | 2.59296, 3.30542, 4.16834, 5.21490, 6.48964, 8.05434, 9.99835, 12.4580, 15.6567, 20.0000}; |
---|
683 | static const G4double IJ[nX]={0, |
---|
684 | .008628, .037206, .089928, .171242, .286114, .440251, .640343, .894382, 1.21208, 1.60544, |
---|
685 | 2.08962, 2.68414, 3.41486, 4.31700, 5.44048, 6.85936, 8.69067, 11.1358, 14.5885, 20.0000}; |
---|
686 | static const G4double* Il[nE]= |
---|
687 | {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ}; |
---|
688 | static const G4double lE[nE]={ |
---|
689 | -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292, |
---|
690 | 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218}; |
---|
691 | static const G4double lEmi=lE[0]; |
---|
692 | static const G4double lEma=lE[nE-1]; |
---|
693 | static const G4double dlE=(lEma-lEmi)/bE; |
---|
694 | //*************************************************************************************** |
---|
695 | G4double Enu=lastE; // Get energy of the last calculated cross-section |
---|
696 | G4double lEn=std::log(Enu); // log(E) for interpolation |
---|
697 | G4double rE=(lEn-lEmi)/dlE; // Position of the energy |
---|
698 | G4int fE=static_cast<int>(rE); // Left bin for interpolation |
---|
699 | if(fE<0) fE=0; |
---|
700 | if(fE>pE)fE=pE; |
---|
701 | G4int sE=fE+1; // Right bin for interpolation |
---|
702 | G4double dE=rE-fE; // relative log shift from the left bin |
---|
703 | G4double dEnu=Enu+Enu; // doubled energy of nu/anu |
---|
704 | G4double Enu2=Enu*Enu; // squared energy of nu/anu |
---|
705 | G4double Ee=Enu-me; // Free Energy of neutrino/anti-neutrino |
---|
706 | G4double ME=Enu*MN; // M*E |
---|
707 | G4double dME=ME+ME; // 2*M*E |
---|
708 | G4double dEMN=(dEnu+MN)*ME; |
---|
709 | G4double MEm=ME-hme2; |
---|
710 | G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2); |
---|
711 | G4double E2M=MN*Enu2-(Enu+MN)*hme2; |
---|
712 | G4double ymax=(E2M+sqE)/dEMN; |
---|
713 | G4double ymin=(E2M-sqE)/dEMN; |
---|
714 | G4double rmin=1.-ymin; |
---|
715 | G4double rhm2E=hme2/Enu2; |
---|
716 | G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu) |
---|
717 | G4double Q2ma=dME*ymax; // Q2_max(E_nu) |
---|
718 | G4double Q2nq=Ee*dMN-mcV; |
---|
719 | if(Q2ma>Q2nq) Q2ma=Q2nq; // Correction for Non Quasi Elastic |
---|
720 | // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma --- |
---|
721 | G4double Rmi=Q2mi/Q2ma; |
---|
722 | G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78); |
---|
723 | // --- E-interpolation must be done in a log scale --- |
---|
724 | G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu) |
---|
725 | G4double Xma=std::pow((shift-1.),power); // X_max(E_nu) |
---|
726 | // Find the integral values integ(Xmi) & integ(Xma) using the direct table |
---|
727 | G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step |
---|
728 | G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum |
---|
729 | G4double rXi=(Xmi-iXmi)/idX; |
---|
730 | G4int iXi=static_cast<int>(rXi); |
---|
731 | if(iXi<0) iXi=0; |
---|
732 | if(iXi>bX) iXi=bX; |
---|
733 | G4double dXi=rXi-iXi; |
---|
734 | G4double bntil=Il[fE][iXi]; |
---|
735 | G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil); |
---|
736 | G4double bntir=Il[sE][iXi]; |
---|
737 | G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir); |
---|
738 | G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral |
---|
739 | // |
---|
740 | G4double rXa=(Xma-iXmi)/idX; |
---|
741 | G4int iXa=static_cast<int>(rXa); |
---|
742 | if(iXa<0) iXa=0; |
---|
743 | if(iXa>bX) iXa=bX; |
---|
744 | G4double dXa=rXa-iXa; |
---|
745 | G4double bntal=Il[fE][iXa]; |
---|
746 | G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal); |
---|
747 | G4double bntar=Il[sE][iXa]; |
---|
748 | G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar); |
---|
749 | G4double inta=intal+dE*(intar-intal);// interpolated end of the integral |
---|
750 | // |
---|
751 | // *** Find X using the reversed table *** |
---|
752 | G4double intx=inti+(inta-inti)*G4UniformRand(); |
---|
753 | G4int intc=static_cast<int>(intx); |
---|
754 | if(intc<0) intc=0; |
---|
755 | if(intc>bX) intc=bX; |
---|
756 | G4double dint=intx-intc; |
---|
757 | G4double mXl=Xl[fE][intc]; |
---|
758 | G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl); |
---|
759 | G4double mXr=Xl[sE][intc]; |
---|
760 | G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr); |
---|
761 | G4double X=Xlb+dE*(Xrb-Xlb); // interpolated X value |
---|
762 | G4double R=shift-std::pow(X,pconv); |
---|
763 | G4double Q2=R*Q2ma; |
---|
764 | return Q2*GeV*GeV; |
---|
765 | } |
---|
766 | |
---|
767 | // It returns a fraction of the direct interaction of the neutrino with quark-partons |
---|
768 | G4double G4QANuENuclearCrossSection::GetDirectPart(G4double Q2) |
---|
769 | { |
---|
770 | G4double f=Q2/4.62; |
---|
771 | G4double ff=f*f; |
---|
772 | G4double r=ff*ff; |
---|
773 | G4double s=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3))); |
---|
774 | //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par) |
---|
775 | return 1.-s*(1.-s/2); |
---|
776 | } |
---|
777 | |
---|
778 | // #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut |
---|
779 | G4double G4QANuENuclearCrossSection::GetNPartons(G4double Q2) |
---|
780 | { |
---|
781 | return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space |
---|
782 | } |
---|
783 | |
---|
784 | // This class can provide only virtual exchange pi- (a substitute for W- boson) |
---|
785 | G4int G4QANuENuclearCrossSection::GetExchangePDGCode() {return -211;} |
---|