Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QNuENuclearCrossSection.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: G4QNuENuclearCrossSection for A(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: nu_e -> e nuclear XS
40 // --------------------------------------------------------------------------------------
41 
42 //#define debug
43 //#define edebug
44 //#define pdebug
45 //#define ppdebug
46 //#define tdebug
47 //#define sdebug
48 
49 #include "G4QNuENuclearCrossSection.hh"
50 #include "G4SystemOfUnits.hh"
51 
52 // Initialization of the
53 G4bool G4QNuENuclearCrossSection::onlyCS=true;// Flag to calculate only CS (not QE)
54 G4double G4QNuENuclearCrossSection::lastSig=0.;// Last calculated total cross section
55 G4double G4QNuENuclearCrossSection::lastQEL=0.;// Last calculated quasi-el. cross section
56 G4int G4QNuENuclearCrossSection::lastL=0; // Last used in cross section TheLastBin
57 G4double G4QNuENuclearCrossSection::lastE=0.; // Last used in cross section TheEnergy
58 G4double* G4QNuENuclearCrossSection::lastEN=0; // Pointer to the Energy Scale of TX & QE
59 G4double* G4QNuENuclearCrossSection::lastTX=0; // Pointer to the LastArray of TX function
60 G4double* G4QNuENuclearCrossSection::lastQE=0; // Pointer to the LastArray of QE function
61 G4int G4QNuENuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
62 G4int G4QNuENuclearCrossSection::lastN=0; // The last N of calculated nucleus
63 G4int G4QNuENuclearCrossSection::lastZ=0; // The last Z of calculated nucleus
64 G4double G4QNuENuclearCrossSection::lastP=0.; // Last used in cross section Momentum
65 G4double G4QNuENuclearCrossSection::lastTH=0.; // Last threshold momentum
66 G4double G4QNuENuclearCrossSection::lastCS=0.; // Last value of the Cross Section
67 G4int G4QNuENuclearCrossSection::lastI=0; // The last position in the DAMDB
68 std::vector<G4double*>* G4QNuENuclearCrossSection::TX = new std::vector<G4double*>;
69 std::vector<G4double*>* G4QNuENuclearCrossSection::QE = new std::vector<G4double*>;
70 
71 // Returns Pointer to the G4VQCrossSection class
72 G4VQCrossSection* G4QNuENuclearCrossSection::GetPointer()
73 {
74  static G4QNuENuclearCrossSection theCrossSection; //**Static body of the Cross Section**
75  return &theCrossSection;
76 }
77 
78 G4QNuENuclearCrossSection::~G4QNuENuclearCrossSection()
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)
90 G4double G4QNuENuclearCrossSection::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"---G4QNENCrossSec::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<<"G4QNENCS::GetCrSec: CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
184 #endif
185  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<<"G4QNENCrossSection::GetCrossSect:NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
192 #endif
193  if(pEn>lastTH)
194  {
195 #ifdef pdebug
196  G4cout<<"G4QNENCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
197 #endif
198  lastTH=pEn;
199  }
200  }
201 #ifdef pdebug
202  G4cout<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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)
259 G4double G4QNuENuclearCrossSection::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)
278 G4double G4QNuENuclearCrossSection::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()33/65!!
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 the neutrinoEnergy 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  if(lastE<=lastEN[sep]) sep-=newran;
347  else sep+=newran;
348  ran=newran;
349  chk=chk+chk;
350  }
351  if(chk+chk!=mL) G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
352  G4double lowE=lastEN[sep];
353  G4double highE=lastEN[sep+1];
354  G4double lowTX=lastTX[sep];
355  if(lastE<lowE||sep>=mL||lastE>highE)
356  G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
357  <<", sep="<<sep<<", mL="<<mL<<G4endl;
358  lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
359  if(!onlyCS) // Skip the differential cross-section parameters
360  {
361  G4double lowQE=lastQE[sep];
362  lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
363 #ifdef pdebug
364  G4cout<<"G4NuENuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
365 #endif
366  }
367  }
368  else
369  {
370  lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
371  lastQEL=lastQE[mL];
372  }
373  if(lastQEL<0.) lastQEL = 0.;
374  if(lastSig<0.) lastSig = 0.;
375  // The cross-sections are expected to be in mb
376  lastSig*=mb38;
377  if(!onlyCS) lastQEL*=mb38;
378  return lastSig;
379 }
380 
381 // Calculate the cros-section functions
382 // ****************************************************************************************
383 // *** This tables are the same for all lepto-nuclear reactions, only mass is different ***
384 // ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ ***
385 // ****************************************************************************************
386 G4int G4QNuENuclearCrossSection::GetFunctions(G4int z, G4int n,
387  G4double* t, G4double* q, G4double* e)
388 {
389  static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
390  static const G4double dmN=mN+mN; // Doubled nucleon mass (2*AtomicMassUnit, GeV)
391  static const G4double me=.00051099892; // electron mass in GeV
392  static const G4double me2=me*me; // Squared mass of an electron in GeV^2
393  static const G4double thresh=me+me2/dmN; // Universal threshold in GeV
394  static const G4int nE=65; // !! If change this, change it in CalculateCrossSection() !!
395  static const G4double nuEn[nE]={thresh,
396  .00051331,.00053602,.00056078,.00058783,.00061743,.00064990,.00068559,.00072492,
397  .00076834,.00081641,.00086975,.00092912,.00099536,.00106950,.00115273,.00124646,
398  .00135235,.00147241,.00160901,.00176503,.00194392,.00214986,.00238797,.00266448,
399  .00298709,.00336531,.00381094,.00433879,.00496745,.00572047,.00662785,.00772806,
400  .00907075,.01072050,.01276190,.01530660,.01850330,.02255110,.02771990,.03437780,
401  .04303240,.05438970,.06944210,.08959920,.11688400,.15423600,.20597200,.27851200,
402  .38153100,.52979600,.74616300,1.0665200,1.5480900,2.2834800,3.4251100,5.2281000,
403  8.1270200,12.875900,20.808500,34.331200,57.877800,99.796200,176.16300,318.68200};
404  static const G4double TOTX[nE]={0.,
405  .00047551,.00162896,.00232785,.00292938,.00349456,.00404939,.00460908,.00518455,
406  .00578488,.00641848,.00709376,.00781964,.00860585,.00946334,.01040460,.01144420,
407  .01259910,.01388920,.01533830,.01697480,.01883260,.02095280,.02338500,.02618960,
408  .02944060,.03322870,.03766580,.04289050,.04907540,.06123530,.07521120,.09034730,
409  .10803800,.12856800,.15277600,.18174900,.21764300,.26083100,.31424900,.37935600,
410  .45871900,.55375100,.66506600,.79039400,.92276600,1.0489000,1.1500300,1.2071700,
411  1.2096800,1.1612200,1.0782900,.98251100,.89137400,.81472700,.75500100,.71061200,
412  .67873900,.65619000,.64098400,.63085100,.62389900,.61664900,.61261000,.60635700};
413  static const G4double QELX[nE]={0.,
414  2.44084e-7,8.73147e-7,1.30540e-6,1.72196e-6,2.15765e-6,2.63171e-6,3.15996e-6,3.75836e-6,
415  4.44474e-6,5.24008e-6,6.16982e-6,7.26537e-6,8.56595e-6,1.01211e-5,1.19937e-5,1.42647e-5,
416  1.70384e-5,2.04506e-5,2.46796e-5,2.99611e-5,3.66090e-5,4.50456e-5,5.58425e-5,6.97817e-5,
417  8.79417e-5,.000111825,.000143542,.000186093,.000243780,.000350295,.000498489,.000698209,
418  .000979989,.001378320,.001949710,.002781960,.004027110,.005882030,.008710940,.013041400,
419  .019739800,.030118300,.046183600,.070818700,.107857000,.161778000,.236873000,.336212000,
420  .455841000,.565128000,.647837000,.701208000,.729735000,.742062000,.746495000,.748182000,
421  .749481000,.750637000,.751471000,.752237000,.752763000,.753103000,.753159000,.753315000};
422 
423  // --------------------------------
424  G4int first=0;
425  if(z<0.)
426  {
427  first=1;
428  z=-z;
429  }
430  if(z<1 || z>92) // neutron and plutonium are forbidden
431  {
432  G4cout<<"***G4QNuENuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
433  return -1;
434  }
435  for(G4int k=0; k<nE; k++)
436  {
437  G4double a=n+z;
438  G4double na=n+a;
439  G4double dn=n+n;
440  G4double da=a+a;
441  G4double ta=da+a;
442  if(first) e[k]=nuEn[k]; // Energy of neutrino E (first bin k=0 can be modified)
443  t[k]=TOTX[k]*nuEn[k]*(na+na)/ta+QELX[k]*(dn+dn-da)/ta; // TotalCrossSection
444  q[k]=QELX[k]*dn/a; // QuasiElasticCrossSection
445  }
446  return first;
447 }
448 
449 // Randomize Q2 from neutrino to the scattered electron when scattering is quasi-elastic
450 G4double G4QNuENuclearCrossSection::GetQEL_ExchangeQ2()
451 {
452  static const G4double me=.00051099892; // electron mass in GeV
453  static const G4double me2=me*me; // Squared mass of an electron in GeV^2
454  static const G4double hme2=me2/2; // .5*m_e^2 in GeV^2
455  static const G4double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
456  static const double MN2=MN*MN; // M_N^2 in GeV^2
457  static const G4double power=-3.5; // direct power for the magic variable
458  static const G4double pconv=1./power;// conversion power for the magic variable
459  static const G4int nQ2=101; // #Of point in the Q2l table (in GeV^2)
460  static const G4int lQ2=nQ2-1; // index of the last in the Q2l table
461  static const G4int bQ2=lQ2-1; // index of the before last in the Q2 ltable
462  // Reversed table
463  static const G4double Xl[nQ2]={1.87905e-10,
464  .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374,
465  .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069,
466  .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718,
467  .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169,
468  .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670,
469  .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720,
470  .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457,
471  .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366,
472  .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128,
473  .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539};
474  // Direct table
475  static const G4double Xmax=Xl[lQ2];
476  static const G4double Xmin=Xl[0];
477  static const G4double dX=(Xmax-Xmin)/lQ2; // step in X(Q2, GeV^2)
478  static const G4double inl[nQ2]={0,
479  1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811,
480  15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039,
481  26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876,
482  36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950,
483  46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377,
484  55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162,
485  65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852,
486  74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767,
487  83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106,
488  92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000};
489  G4double Enu=lastE; // Get energy of the last calculated cross-section
490  G4double dEnu=Enu+Enu; // doubled energy of nu/anu
491  G4double Enu2=Enu*Enu; // squared energy of nu/anu
492  G4double ME=Enu*MN; // M*E
493  G4double dME=ME+ME; // 2*M*E
494  G4double dEMN=(dEnu+MN)*ME;
495  G4double MEm=ME-hme2;
496  G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
497  G4double E2M=MN*Enu2-(Enu+MN)*hme2;
498  G4double ymax=(E2M+sqE)/dEMN;
499  G4double ymin=(E2M-sqE)/dEMN;
500  G4double rmin=1.-ymin;
501  G4double rhm2E=hme2/Enu2;
502  G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
503  G4double Q2ma=dME*ymax; // Q2_max(E_nu)
504  G4double Xma=std::pow((1.+Q2mi),power); // X_max(E_nu)
505  G4double Xmi=std::pow((1.+Q2ma),power); // X_min(E_nu)
506  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
507  G4double rXi=(Xmi-Xmin)/dX;
508  G4int iXi=static_cast<int>(rXi);
509  if(iXi<0) iXi=0;
510  if(iXi>bQ2) iXi=bQ2;
511  G4double dXi=rXi-iXi;
512  G4double bnti=inl[iXi];
513  G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
514  //
515  G4double rXa=(Xma-Xmin)/dX;
516  G4int iXa=static_cast<int>(rXa);
517  if(iXa<0) iXa=0;
518  if(iXa>bQ2) iXa=bQ2;
519  G4double dXa=rXa-iXa;
520  G4double bnta=inl[iXa];
521  G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
522  // *** Find X using the reversed table ***
523  G4double intx=inti+(inta-inti)*G4UniformRand();
524  G4int intc=static_cast<int>(intx);
525  if(intc<0) intc=0;
526  if(intc>bQ2) intc=bQ2; // If it is more than max, then the BAD extrapolation
527  G4double dint=intx-intc;
528  G4double mX=Xl[intc];
529  G4double X=mX+dint*(Xl[intc+1]-mX);
530  G4double Q2=std::pow(X,pconv)-1.;
531  return Q2*GeV*GeV;
532 }
533 
534 // Randomize Q2 from neutrino to the scattered electron when scattering is not quasiElastic
535 G4double G4QNuENuclearCrossSection::GetNQE_ExchangeQ2()
536 {
537  static const double mpi=.13957018; // charged pi meson mass in GeV
538  static const G4double me=.00051099892;// electron mass in GeV
539  static const G4double me2=me*me; // Squared mass of an electron in GeV^2
540  static const G4double hme2=me2/2; // .5*m_e^2 in GeV^2
541  static const double MN=.931494043; // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
542  static const double MN2=MN*MN; // M_N^2 in GeV^2
543  static const double dMN=MN+MN; // 2*M_N in GeV
544  static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
545  static const G4int power=7; // direct power for the magic variable
546  static const G4double pconv=1./power; // conversion power for the magic variable
547  static const G4int nX=21; // #Of point in the Xl table (in GeV^2)
548  static const G4int lX=nX-1; // index of the last in the Xl table
549  static const G4int bX=lX-1; // @@ index of the before last in the Xl table
550  static const G4int nE=20; // #Of point in the El table (in GeV^2)
551  static const G4int bE=nE-1; // index of the last in the El table
552  static const G4int pE=bE-1; // index of the before last in the El table
553  // Reversed table
554  static const G4double X0[nX]={6.14081e-05,
555  .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244,
556  2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419};
557  static const G4double X1[nX]={.00125268,
558  .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049,
559  4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330};
560  static const G4double X2[nX]={.015694,
561  1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728,
562  10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276};
563  static const G4double X3[nX]={.0866877,
564  4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302,
565  21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115};
566  static const G4double X4[nX]={.160483,
567  5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127,
568  29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265};
569  static const G4double X5[nX]={.0999307,
570  5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398,
571  26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649};
572  static const G4double X6[nX]={.0276367,
573  3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519,
574  16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148};
575  static const G4double X7[nX]={.00472383,
576  2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473,
577  9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057};
578  static const G4double X8[nX]={.000630783,
579  1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480,
580  5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037};
581  static const G4double X9[nX]={7.49179e-05,
582  .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823,
583  3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028};
584  static const G4double XA[nX]={8.43437e-06,
585  .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704,
586  2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823};
587  static const G4double XB[nX]={9.27028e-07,
588  .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059,
589  1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341};
590  static const G4double XC[nX]={1.00807e-07,
591  .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668,
592  1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014};
593  static const G4double XD[nX]={1.09102e-08,
594  .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878,
595  1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617};
596  static const G4double XE[nX]={1.17831e-09,
597  .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374,
598  .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571};
599  static const G4double XF[nX]={1.27141e-10,
600  .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046,
601  .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338};
602  static const G4double XG[nX]={1.3713e-11,
603  .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870,
604  .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379};
605  static const G4double XH[nX]={1.47877e-12,
606  .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421,
607  .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197};
608  static const G4double XI[nX]={1.59454e-13,
609  .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806,
610  .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872};
611  static const G4double XJ[nX]={1.71931e-14,
612  .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147,
613  .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813};
614  // Direct table
615  static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
616  X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
617  static const G4double dX[nE]={
618  (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
619  (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
620  (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
621  (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
622  (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
623  static const G4double* Xl[nE]=
624  {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
625  static const G4double I0[nX]={0,
626  .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555,
627  13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000};
628  static const G4double I1[nX]={0,
629  .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365,
630  13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000};
631  static const G4double I2[nX]={0,
632  .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586,
633  12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000};
634  static const G4double I3[nX]={0,
635  .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461,
636  12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000};
637  static const G4double I4[nX]={0,
638  .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515,
639  11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000};
640  static const G4double I5[nX]={0,
641  .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040,
642  10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000};
643  static const G4double I6[nX]={0,
644  .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243,
645  9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000};
646  static const G4double I7[nX]={0,
647  .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717,
648  8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000};
649  static const G4double I8[nX]={0,
650  .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415,
651  8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000};
652  static const G4double I9[nX]={0,
653  .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672,
654  8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000};
655  static const G4double IA[nX]={0,
656  .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265,
657  7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000};
658  static const G4double IB[nX]={0,
659  .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298,
660  7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000};
661  static const G4double IC[nX]={0,
662  .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696,
663  7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000};
664  static const G4double ID[nX]={0,
665  .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212,
666  7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000};
667  static const G4double IE[nX]={0,
668  .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949,
669  7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000};
670  static const G4double IF[nX]={0,
671  .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164,
672  7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000};
673  static const G4double IG[nX]={0,
674  .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249,
675  7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000};
676  static const G4double IH[nX]={0,
677  .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804,
678  7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000};
679  static const G4double II[nX]={0,
680  .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722,
681  7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000};
682  static const G4double IJ[nX]={0,
683  .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196,
684  6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000};
685  static const G4double* Il[nE]=
686  {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
687  static const G4double lE[nE]={
688 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
689  2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
690  static const G4double lEmi=lE[0];
691  static const G4double lEma=lE[nE-1];
692  static const G4double dlE=(lEma-lEmi)/bE;
693  //***************************************************************************************
694  G4double Enu=lastE; // Get energy of the last calculated cross-section
695  G4double lEn=std::log(Enu); // log(E) for interpolation
696  G4double rE=(lEn-lEmi)/dlE; // Position of the energy
697  G4int fE=static_cast<int>(rE); // Left bin for interpolation
698  if(fE<0) fE=0;
699  if(fE>pE)fE=pE;
700  G4int sE=fE+1; // Right bin for interpolation
701  G4double dE=rE-fE; // relative log shift from the left bin
702  G4double dEnu=Enu+Enu; // doubled energy of nu/anu
703  G4double Enu2=Enu*Enu; // squared energy of nu/anu
704  G4double Ee=Enu-me; // Free Energy of neutrino/anti-neutrino
705  G4double ME=Enu*MN; // M*E
706  G4double dME=ME+ME; // 2*M*E
707  G4double dEMN=(dEnu+MN)*ME;
708  G4double MEm=ME-hme2;
709  G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
710  G4double E2M=MN*Enu2-(Enu+MN)*hme2;
711  G4double ymax=(E2M+sqE)/dEMN;
712  G4double ymin=(E2M-sqE)/dEMN;
713  G4double rmin=1.-ymin;
714  G4double rhm2E=hme2/Enu2;
715  G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
716  G4double Q2ma=dME*ymax; // Q2_max(E_nu)
717  G4double Q2nq=Ee*dMN-mcV;
718  if(Q2ma>Q2nq) Q2ma=Q2nq; // Correction for Non Quasi Elastic
719  // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
720  G4double Rmi=Q2mi/Q2ma;
721  G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu
722  // --- E-interpolation must be done in a log scale ---
723  G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
724  G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
725  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
726  G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
727  G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
728  G4double rXi=(Xmi-iXmi)/idX;
729  G4int iXi=static_cast<int>(rXi);
730  if(iXi<0) iXi=0;
731  if(iXi>bX) iXi=bX;
732  G4double dXi=rXi-iXi;
733  G4double bntil=Il[fE][iXi];
734  G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
735  G4double bntir=Il[sE][iXi];
736  G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
737  G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
738  //
739  G4double rXa=(Xma-iXmi)/idX;
740  G4int iXa=static_cast<int>(rXa);
741  if(iXa<0) iXa=0;
742  if(iXa>bX) iXa=bX;
743  G4double dXa=rXa-iXa;
744  G4double bntal=Il[fE][iXa];
745  G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
746  G4double bntar=Il[sE][iXa];
747  G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
748  G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
749  //
750  // *** Find X using the reversed table ***
751  G4double intx=inti+(inta-inti)*G4UniformRand();
752  G4int intc=static_cast<int>(intx);
753  if(intc<0) intc=0;
754  if(intc>bX) intc=bX;
755  G4double dint=intx-intc;
756  G4double mXl=Xl[fE][intc];
757  G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
758  G4double mXr=Xl[sE][intc];
759  G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
760  G4double X=Xlb+dE*(Xrb-Xlb); // interpolated X value
761  G4double R=shift-std::pow(X,pconv);
762  G4double Q2=R*Q2ma;
763  return Q2*GeV*GeV;
764 }
765 
766 // It returns a fraction of the direct interaction of the neutrino with quark-partons
767 G4double G4QNuENuclearCrossSection::GetDirectPart(G4double Q2)
768 {
769  G4double f=Q2/4.62;
770  G4double ff=f*f;
771  G4double r=ff*ff;
772  G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
773  //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
774  return 1.-s_value*(1.-s_value/2);
775 }
776 
777 // #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
778 G4double G4QNuENuclearCrossSection::GetNPartons(G4double Q2)
779 {
780  return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
781 }
782 
783 // This class can provide only virtual exchange pi+ (a substitute for W+ boson)
784 G4int G4QNuENuclearCrossSection::GetExchangePDGCode() {return 211;}
Generated on Sat May 25 2013 14:33:49 for Geant4 by   doxygen 1.8.4