source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/cross_sections/src/G4QNuENuclearCrossSection.cc @ 1340

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

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