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