source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/cross_sections/src/G4QANuENuclearCrossSection.cc @ 1353

Last change on this file since 1353 was 1340, checked in by garnier, 14 years ago

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