source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QANuENuclearCrossSection.cc @ 1055

Last change on this file since 1055 was 1055, checked in by garnier, 15 years ago

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