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