source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QANuANuNuclearCrossSection.cc @ 968

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

update processes
File size: 38.6 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: G4QANuANuNuclearCrossSection.cc,v 1.2 2007/11/02 15:57:16 mkossov Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31// G4 Physics class: G4QANuANuNuclearCrossSection for gamma+A cross sections
32// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
33// The last update: M.V. Kossov, CERN/ITEP (Moscow) 17-Oct-03
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 "G4QANuANuNuclearCrossSection.hh"
49
50// Initialization of the
51G4bool    G4QANuANuNuclearCrossSection::onlyCS=true;//Flag to calculate only CS (not QE)
52G4double  G4QANuANuNuclearCrossSection::lastSig=0.;//Last calculated total cross section
53G4double  G4QANuANuNuclearCrossSection::lastQEL=0.;//Last calculated quasi-el cross section
54G4int     G4QANuANuNuclearCrossSection::lastL=0;   //Last used in cross section TheLastBin
55G4double  G4QANuANuNuclearCrossSection::lastE=0.;  //Last used in cross section TheEnergy
56G4double* G4QANuANuNuclearCrossSection::lastEN=0;  //Pointer to theEnergy Scale of TX & QE
57G4double* G4QANuANuNuclearCrossSection::lastTX=0;  //Pointer to theLastArray of TX function
58G4double* G4QANuANuNuclearCrossSection::lastQE=0;  //Pointer to theLastArray of QE function
59G4int     G4QANuANuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
60G4int     G4QANuANuNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
61G4int     G4QANuANuNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
62G4double  G4QANuANuNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
63G4double  G4QANuANuNuclearCrossSection::lastTH=0.; // Last threshold momentum
64G4double  G4QANuANuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
65G4int     G4QANuANuNuclearCrossSection::lastI=0;   // The last position in the DAMDB
66
67// Returns Pointer to the G4VQCrossSection class
68G4VQCrossSection* G4QANuANuNuclearCrossSection::GetPointer()
69{
70  static G4QANuANuNuclearCrossSection 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 G4QANuANuNuclearCrossSection::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<<"G4QAMNCS::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!=-14)
95  {
96#ifdef debug
97    G4cout<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"---G4QAMNCrossSec::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<<"G4QAMNCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<lastI<<G4endl;
170#endif
171      //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
172      lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
173      if(lastCS<=0.)
174                                                {
175        lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
176#ifdef pdebug
177        G4cout<<"G4QAMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
178#endif
179        if(pEn>lastTH)
180        {
181#ifdef pdebug
182          G4cout<<"G4QAMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
183#endif
184          lastTH=pEn;
185        }
186                                                }
187#ifdef pdebug
188      G4cout<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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 G4QANuANuNuclearCrossSection::ThresholdEnergy(G4int, G4int, G4int) {return 0.;}
246
247// The main member function giving the gamma-A cross section (E_kin in MeV, CS in mb)
248G4double G4QANuANuNuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I,
249                                       G4int , G4int targZ, G4int targN, G4double Momentum)
250{
251  static const G4double mb38=1.E-11;// Conversion 10^-38 cm^2 to mb=10^-27 cm^2
252  static const G4int nE=65;   // !! If change this, change it in GetFunctions() (*.hh) !!
253  static const G4int mL=nE-1;
254  static const G4double EMi=0.; // Universal threshold of the reaction in GeV
255  static const G4double EMa=300.;       // Maximum tabulated Energy of nu_mu in GeV
256  // *** Begin of the Associative memory for acceleration of the cross section calculations
257  static std::vector <G4double> colH;   //?? Vector of HighEnergyCoefficients (functional)
258  static std::vector <G4double*> TX;    // Vector of pointers to the TX tabulated functions
259  static std::vector <G4double*> QE;    // Vector of pointers to the QE tabulated functions
260  static G4bool first=true;             // Flag of initialization of the energy axis
261  // *** End of Static Definitions (Associative Memory) ***
262  //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
263  //G4double TotEnergy2=Momentum;
264  onlyCS=CS;                            // Flag to calculate only CS (not TX & QE)
265  lastE=Momentum/GeV;                   // Kinetic energy of the muon neutrino (in GeV!)
266  if (lastE<=EMi)                       // Energy is below the minimum energy in the table
267  {
268    lastE=0.;
269    lastSig=0.;
270    return 0.;
271  }
272  G4int Z=targZ;                        // New Z, which can change the sign
273  if(F<=0)                              // This isotope was not the last used isotop
274  {
275    if(F<0)                          // This isotope was found in DAMDB =========> RETRIEVE
276                                {
277      lastTX =TX[I];                 // Pointer to the prepared TX function (same isotope)
278      lastQE =QE[I];                 // Pointer to the prepared QE function (same isotope)
279          }
280          else                              // This isotope wasn't calculated previously => CREATE
281          {
282      if(first)
283      {
284        lastEN = new G4double[nE];   // This must be done only once!
285        Z=-Z;                        // To explain GetFunctions that E-axis must be filled
286        first=false;                 // To make it only once
287      }
288      lastTX = new G4double[nE];     // Allocate memory for the new TX function
289      lastQE = new G4double[nE];     // Allocate memory for the new QE function
290      G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
291      if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
292      // *** The synchronization check ***
293      G4int sync=TX.size();
294      if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
295      TX.push_back(lastTX);
296      QE.push_back(lastQE);
297           } // End of creation of the new set of parameters
298  } // End of parameters udate
299  // ============================== NOW Calculate the Cross Section =====================
300  if (lastE<=EMi)                    // Check that antiNuEnergy is higher than ThreshE
301  {
302    lastE=0.;
303    lastSig=0.;
304    return 0.;
305  }
306  if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
307  {
308    G4int chk=1;
309    G4int ran=mL/2;
310    G4int sep=ran;  // as a result = an index of the left edge of the interval
311    while(ran>=2)
312                                {
313      G4int newran=ran/2;
314      if(lastE<=lastEN[sep]) sep-=newran;
315      else                   sep+=newran;
316      ran=newran;
317      chk=chk+chk; 
318    }
319    if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
320    G4double lowE=lastEN[sep];
321    G4double highE=lastEN[sep+1];
322    G4double lowTX=lastTX[sep];
323    if(lastE<lowE||sep>=mL||lastE>highE)
324      G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
325            <<", sep="<<sep<<", mL="<<mL<<G4endl;
326    lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
327    if(!onlyCS)                       // Skip the differential cross-section parameters
328    {
329      G4double lowQE=lastQE[sep];
330      lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
331#ifdef pdebug
332      G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
333#endif
334    }
335  }
336  else
337  {
338    lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
339    lastQEL=lastQE[mL];
340  }
341  if(lastQEL<0.) lastQEL = 0.;
342  if(lastSig<0.) lastSig = 0.;
343  // The cross-sections are expected to be in mb
344  lastSig*=mb38;
345  if(!onlyCS) lastQEL*=mb38;
346  return lastSig;
347}
348
349// Calculate the cros-section functions
350// ****************************************************************************************
351// *** This tables are the same for all lepto-nuclear reactions, only mass is different ***
352// ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ ***
353// ****************************************************************************************
354G4int G4QANuANuNuclearCrossSection::GetFunctions (G4int z, G4int n,
355                                                     G4double* t, G4double* q, G4double* e)
356{
357  static const G4int nE=65; // !! If change this, change it in GetCrossSection() (*.cc) !!
358  static const G4double nuEn[nE]={0.,
3591.00463e-5,1.05336e-5,1.10692e-5,1.16592e-5,1.23109e-5,1.30323e-5,1.38331e-5,1.47245e-5,
3601.57194e-5,1.68335e-5,1.80848e-5,1.94948e-5,2.10894e-5,2.28991e-5,2.49608e-5,2.73189e-5,
3613.00273e-5,3.31516e-5,3.67722e-5,4.09881e-5,4.59217e-5,5.17255e-5,5.85908e-5,6.67583e-5,
3627.65338e-5,8.83078e-5,.000102583,.000120011,.000141441,.000167995,.000201160,.000242926,
363.000295985,.000364008,.000452051,.000567152,.000719210,.000922307,.001196710,.001571930,
364.002091530,.002820590,.003857810,.005354930,.007548840,.010815300,.015760100,.023376900,
365.035325600,.054430800,.085595700,.137508000,.225898000,.379892000,.654712000,1.15767000,
3662.10277000,3.92843000,7.55861000,14.9991000,30.7412000,65.1734000,143.155000,326.326000};
367  static const G4double TOTX[nE]={0.,
3683.63538e-5,3.81165e-5,4.00539e-5,4.21884e-5,4.45454e-5,4.71548e-5,5.00511e-5,5.32747e-5,
3695.68730e-5,6.09016e-5,6.54261e-5,7.05246e-5,7.62894e-5,8.28315e-5,9.02838e-5,9.88066e-5,
370.000108594,.000119884,.000132964,.000148191,.000166007,.000186959,.000211736,.000241202,
371.000276453,.000318890,.000370311,.000433042,.000510116,.000605516,.000724514,.000874144,
372.001063870,.001306520,.001619660,.002027520,.002563780,.003275710,.004230080,.005521890,
373.007286920,.009719890,.013099700,.018695100,.029208400,.042095000,.059253700,.082373900,
374.113071000,.151041000,.191803000,.224208000,.234187000,.217774000,.187139000,.157818000,
375.137494000,.125872000,.120462000,.119148000,.120418000,.123027000,.126408000,.129071000};
376  static const G4double QELX[nE]={0.,
377.365220e-9,.401502e-9,.443364e-9,.491885e-9,.548393e-9,.614536e-9,.692362e-9,.784441e-9,
378.894012e-9,1.02519e-9,1.18322e-9,1.37487e-9,1.60890e-9,1.89677e-9,2.25355e-9,2.69928e-9,
3793.26079e-9,3.97433e-9,4.88937e-9,6.07407e-9,7.62331e-9,9.67058e-9,1.24058e-8,1.61022e-8,
3802.11580e-8,2.81605e-8,3.79876e-8,5.19696e-8,7.21515e-8,1.01724e-7,1.45743e-7,2.12353e-7,
3813.14890e-7,4.75585e-7,7.32171e-7,1.14991e-6,1.84390e-6,3.02121e-6,5.06217e-6,8.68002e-6,
3821.52408e-5,2.74159e-5,5.05363e-5,.000100111,.000220489,.000455269,.000933841,.001925650,
383.003994300,.008221270,.016417600,.030830400,.052902400,.082519200,.115560000,.149598000,
384.184112000,.215102000,.238253000,.252949000,.261267000,.265626000,.267782000,.268791000};
385  // --------------------------------
386  G4int first=0;
387  if(z<0.)
388                {
389    first=1;
390    z=-z;
391  }
392  if(z<1 || z>92)             // neutron & plutonium are forbidden
393  {
394    G4cout<<"*G4QANuANuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
395    return -1;
396  }
397  for(G4int k=0; k<nE; k++)
398  {
399    G4double a=n+z;
400    G4double za=z+a;
401    G4double dz=z+z;
402    G4double da=a+a;
403    G4double ta=da+a;
404    if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
405    t[k]=TOTX[k]*nuEn[k]*(za+za)/ta+QELX[k]*(dz+dz-da)/ta; // TotalCrossSection
406    q[k]=QELX[k]*dz/a;                                     // QuasiElasticCrossSection
407         }
408  return first;
409}
410
411// Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
412G4double G4QANuANuNuclearCrossSection::GetQEL_ExchangeQ2()
413{
414  static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
415  static const G4double power=-3.5;    // direct power for the magic variable
416  static const G4double pconv=1./power;// conversion power for the magic variable
417  static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
418  static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
419  static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
420  // Reversed table
421  static const G4double Xl[nQ2]={5.20224e-16,
422 .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898,
423        .093878,        .103332,        .112844,        .122410,        .132023,        .141680,        .151376,        .161109,        .170875,        .180672,
424        .190499,        .200352,        .210230,        .220131,        .230055,        .239999,        .249963,        .259945,        .269944,        .279960,
425        .289992,        .300039,        .310099,        .320173,        .330260,        .340359,        .350470,        .360592,        .370724,        .380867,
426        .391019,        .401181,        .411352,        .421531,        .431719,        .441915,        .452118,        .462329,        .472547,        .482771,
427        .493003,        .503240,        .513484,        .523734,        .533989,        .544250,        .554517,        .564788,        .575065,        .585346,
428        .595632,        .605923,        .616218,        .626517,        .636820,        .647127,        .657438,        .667753,        .678072,        .688394,
429        .698719,        .709048,        .719380,        .729715,        .740053,        .750394,        .760738,        .771085,        .781434,        .791786,
430        .802140,        .812497,        .822857,        .833219,        .843582,        .853949,        .864317,        .874687,        .885060,        .895434,
431        .905810,        .916188,        .926568,        .936950,        .947333,        .957719,        .968105,        .978493,        .988883,        .999275};
432          // Direct table
433  static const G4double Xmax=Xl[lQ2];
434  static const G4double Xmin=Xl[0];
435  static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
436  static const G4double inl[nQ2]={0,
437        1.52225,        2.77846,        3.96651,        5.11612,        6.23990,        7.34467,        8.43466,        9.51272,        10.5809,        11.6406,
438        12.6932,        13.7394,        14.7801,        15.8158,        16.8471,        17.8743,        18.8979,        19.9181,        20.9353,        21.9496,
439        22.9614,        23.9707,        24.9777,        25.9826,        26.9855,        27.9866,        28.9860,        29.9837,        30.9798,        31.9745,
440        32.9678,        33.9598,        34.9505,        35.9400,        36.9284,        37.9158,        38.9021,        39.8874,        40.8718,        41.8553,
441        42.8379,        43.8197,        44.8007,        45.7810,        46.7605,        47.7393,        48.7174,        49.6950,        50.6718,        51.6481,
442        52.6238,        53.5990,        54.5736,        55.5476,        56.5212,        57.4943,        58.4670,        59.4391,        60.4109,        61.3822,
443        62.3531,        63.3236,        64.2937,        65.2635,        66.2329,        67.2019,        68.1707,        69.1390,        70.1071,        71.0748,
444        72.0423,        73.0095,        73.9763,        74.9429,        75.9093,        76.8754,        77.8412,        78.8068,        79.7721,        80.7373,
445        81.7022,        82.6668,        83.6313,        84.5956,        85.5596,        86.5235,        87.4872,        88.4507,        89.4140,        90.3771,
446        91.3401,        92.3029,        93.2656,        94.2281,        95.1904,        96.1526,        97.1147,        98.0766,        99.0384,        100.000};
447  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
448  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
449  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
450  G4double ME=Enu*MN;                 // M*E
451  G4double dME=ME+ME;                 // 2*M*E
452  G4double dEMN=(dEnu+MN)*ME;
453  G4double sqE=Enu*ME;
454  G4double E2M=MN*Enu2;
455  G4double ymax=(E2M+sqE)/dEMN;
456  G4double Q2mi=0.; // Q2_min(E_nu)
457  G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
458  G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
459  G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
460  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
461  G4double rXi=(Xmi-Xmin)/dX;
462  G4int    iXi=static_cast<int>(rXi);
463  if(iXi<0) iXi=0;
464  if(iXi>bQ2) iXi=bQ2;
465  G4double dXi=rXi-iXi;
466  G4double bnti=inl[iXi];
467  G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
468  //
469  G4double rXa=(Xma-Xmin)/dX;
470  G4int    iXa=static_cast<int>(rXa);
471  if(iXa<0) iXa=0;
472  if(iXa>bQ2) iXa=bQ2;
473  G4double dXa=rXa-iXa;
474  G4double bnta=inl[iXa];
475  G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
476  // *** Find X using the reversed table ***
477  G4double intx=inti+(inta-inti)*G4UniformRand();
478  G4int    intc=static_cast<int>(intx);
479  if(intc<0) intc=0;
480  if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
481  G4double dint=intx-intc;
482  G4double mX=Xl[intc];
483  G4double X=mX+dint*(Xl[intc+1]-mX);
484  G4double Q2=std::pow(X,pconv)-1.;
485  return Q2*GeV*GeV;
486}
487
488// Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
489G4double G4QANuANuNuclearCrossSection::GetNQE_ExchangeQ2()
490{
491  static const double mpi=.13957018;    // charged pi meson mass in GeV
492  static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
493  static const double dMN=MN+MN;        // 2*M_N in GeV
494  static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
495  static const G4int power=7;           // direct power for the magic variable
496  static const G4double pconv=1./power; // conversion power for the magic variable
497  static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
498  static const G4int lX=nX-1;           // index of the last in the Xl table
499  static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
500  static const G4int nE=20;             // #Of point in the El table (in GeV^2)
501  static const G4int bE=nE-1;           // index of the last in the El table
502  static const G4int pE=bE-1;           // index of the before last in the El table
503  // Reversed table
504  static const G4double X0[nX]={5.21412e-05,
505        .437860,        .681908,        .891529,        1.08434,        1.26751,        1.44494,        1.61915,        1.79198,        1.96493,        2.13937,
506        2.31664,        2.49816,        2.68559,        2.88097,        3.08705,        3.30774,        3.54917,        3.82233,        4.15131, 4.62182};
507  static const G4double X1[nX]={.00102591,
508        1.00443,        1.55828,        2.03126,        2.46406,        2.87311,        3.26723,        3.65199,        4.03134,        4.40835,        4.78561,
509        5.16549,        5.55031,        5.94252,        6.34484,        6.76049,        7.19349,        7.64917,        8.13502,        8.66246, 9.25086};
510  static const G4double X2[nX]={.0120304,
511        2.59903,        3.98637,        5.15131,        6.20159,        7.18024,        8.10986,        9.00426,        9.87265,        10.7217,        11.5564,
512        12.3808,        13.1983,        14.0116,        14.8234,        15.6359,        16.4515,        17.2723,        18.1006,        18.9386, 19.7892};
513  static const G4double X3[nX]={.060124,
514        5.73857,        8.62595,        10.9849,        13.0644,        14.9636,        16.7340,        18.4066,        20.0019,        21.5342,        23.0142,
515        24.4497,        25.8471,        27.2114,        28.5467,        29.8564,        31.1434,        32.4102,        33.6589,        34.8912, 36.1095};
516  static const G4double X4[nX]={.0992363,
517        8.23746,        12.1036,        15.1740,        17.8231,        20.1992,        22.3792,        24.4092,        26.3198,        28.1320,        29.8615,
518        31.5200,        33.1169,        34.6594,        36.1536,        37.6044,        39.0160,        40.3920,        41.7353,        43.0485, 44.3354};
519  static const G4double X5[nX]={.0561127,
520        7.33661,        10.5694,        13.0778,        15.2061,        17.0893,        18.7973,        20.3717,        21.8400,        23.2211,        24.5291,
521        25.7745,        26.9655,        28.1087,        29.2094,        30.2721,        31.3003,        32.2972,        33.2656,        34.2076, 35.1265};
522  static const G4double X6[nX]={.0145859,
523        4.81774,        6.83565,        8.37399,        9.66291,        10.7920,        11.8075,        12.7366,        13.5975,        14.4025,        15.1608,
524        15.8791,        16.5628,        17.2162,        17.8427,        18.4451,        19.0259,        19.5869,        20.1300,        20.6566, 21.1706};
525  static const G4double X7[nX]={.00241155,
526        2.87095,        4.02492,        4.89243,        5.61207,        6.23747,        6.79613,        7.30433,        7.77270,        8.20858,        8.61732,
527        9.00296,        9.36863,        9.71682,        10.0495,        10.3684,        10.6749,        10.9701,        11.2550,        11.5306, 11.7982};
528  static const G4double X8[nX]={.000316863,
529        1.76189,        2.44632,        2.95477,        3.37292,        3.73378,        4.05420,        4.34415,        4.61009,        4.85651,        5.08666,
530        5.30299,        5.50738,        5.70134,        5.88609,        6.06262,        6.23178,        6.39425,        6.55065,        6.70149, 6.84742};
531  static const G4double X9[nX]={3.73544e-05,
532        1.17106,        1.61289,        1.93763,        2.20259,        2.42976,        2.63034,        2.81094,        2.97582,        3.12796,        3.26949,
533        3.40202,        3.52680,        3.64482,        3.75687,        3.86360,        3.96557,        4.06323,        4.15697,        4.24713, 4.33413};
534  static const G4double XA[nX]={4.19131e-06,
535        .849573,        1.16208,        1.38955,        1.57379,        1.73079,        1.86867,        1.99221,        2.10451,        2.20770,        2.30332,
536        2.39252,        2.47622,        2.55511,        2.62977,        2.70066,        2.76818,        2.83265,        2.89437,        2.95355, 3.01051};
537  static const G4double XB[nX]={4.59981e-07,
538        .666131,        .905836,        1.07880,        1.21796,        1.33587,        1.43890,        1.53080,        1.61399,        1.69011,        1.76040,
539        1.82573,        1.88682,        1.94421,        1.99834,        2.04959,        2.09824,        2.14457,        2.18878,        2.23107, 2.27162};
540  static const G4double XC[nX]={4.99861e-08,
541        .556280,        .752730,        .893387,        1.00587,        1.10070,        1.18317,        1.25643,        1.32247,        1.38269,        1.43809,
542        1.48941,        1.53724,        1.58203,        1.62416,        1.66391,        1.70155,        1.73728,        1.77128,        1.80371, 1.83473};
543  static const G4double XD[nX]={5.40832e-09,
544        .488069,        .657650,        .778236,        .874148,        .954621,        1.02432,        1.08599,        1.14138,        1.19172,        1.23787,
545        1.28049,        1.32008,        1.35705,        1.39172,        1.42434,        1.45514,        1.48429,        1.51197,        1.53829, 1.56339};
546  static const G4double XE[nX]={5.84029e-10,
547        .445057,        .597434,        .705099,        .790298,        .861468,        .922865,        .976982,        1.02542,        1.06930,        1.10939,
548        1.14630,        1.18050,        1.21233,        1.24208,        1.27001,        1.29630,        1.32113,        1.34462,        1.36691, 1.38812};
549  static const G4double XF[nX]={6.30137e-11,
550        .418735,        .560003,        .659168,        .737230,        .802138,        .857898,        .906854,        .950515,        .989915,        1.02580,
551        1.05873,        1.08913,        1.11734,        1.14364,        1.16824,        1.19133,        1.21306,        1.23358,        1.25298, 1.27139};
552  static const G4double XG[nX]={6.79627e-12,
553        .405286,        .539651,        .633227,        .706417,        .766929,        .818642,        .863824,        .903931,        .939963,        .972639,
554        1.00250,        1.02995,        1.05532,        1.07887,        1.10082,        1.12134,        1.14058,        1.15867,        1.17572, 1.19183};
555  static const G4double XH[nX]={7.32882e-13,
556        .404391,        .535199,        .625259,        .695036,        .752243,        .800752,        .842823,        .879906,        .912994,        .942802,
557        .969862,        .994583,        1.01729,        1.03823,        1.05763,        1.07566,        1.09246,        1.10816,        1.12286, 1.13667};
558  static const G4double XI[nX]={7.90251e-14,
559        .418084,        .548382,        .636489,        .703728,        .758106,        .803630,        .842633,        .876608,        .906576,        .933269,
560        .957233,        .978886,        .998556,        1.01651,        1.03295,        1.04807,        1.06201,        1.07489,        1.08683, 1.09792};
561  static const G4double XJ[nX]={8.52083e-15,
562        .447299,        .579635,        .666780,        .731788,        .783268,        .825512,        .861013,        .891356,        .917626,        .940597,
563        .960842,        .978802,        .994820,        1.00917,        1.02208,        1.03373,        1.04427,        1.05383,        1.06253, 1.07046};
564  // Direct table
565  static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
566                        X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
567  static const G4double dX[nE]={
568    (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
569    (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
570    (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
571    (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
572    (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
573  static const G4double* Xl[nE]=
574                             {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
575  static const G4double I0[nX]={0,
576        .354631,        1.08972,        2.05138,        3.16564,        4.38343,        5.66828,        6.99127,        8.32858,        9.65998,        10.9680,
577        12.2371,        13.4536,        14.6050,        15.6802,        16.6686,        17.5609,        18.3482,        19.0221,        19.5752,        20.0000};
578  static const G4double I1[nX]={0,
579        .281625,        .877354,        1.67084,        2.60566,        3.64420,        4.75838,        5.92589,        7.12829,        8.34989,        9.57708,
580        10.7978,        12.0014,        13.1781,        14.3190,        15.4162,        16.4620,        17.4496,        18.3724,        19.2245,        20.0000};
581  static const G4double I2[nX]={0,
582        .201909,        .642991,        1.24946,        1.98463,        2.82370,        3.74802,        4.74263,        5.79509,        6.89474,        8.03228,
583        9.19947,        10.3889,        11.5938,        12.8082,        14.0262,        15.2427,        16.4527,        17.6518,        18.8356,        20.0000};
584  static const G4double I3[nX]={0,
585        .140937,        .461189,        .920216,        1.49706,        2.17728,        2.94985,        3.80580,        4.73758,        5.73867,        6.80331,
586        7.92637,        9.10316,        10.3294,        11.6013,        12.9150,        14.2672,        15.6548,        17.0746,        18.5239,        20.0000};
587  static const G4double I4[nX]={0,
588        .099161,        .337358,        .694560,        1.16037,        1.72761,        2.39078,        3.14540,        3.98768,        4.91433,        5.92245,
589        7.00942,        8.17287,        9.41060,        10.7206,        12.1010,        13.5500,        15.0659,        16.6472,        18.2924,        20.0000};
590  static const G4double I5[nX]={0,
591        .071131,        .255084,        .543312,        .932025,        1.41892,        2.00243,        2.68144,        3.45512,        4.32283,        5.28411,
592        6.33859,        7.48602,        8.72621,        10.0590,        11.4844,        13.0023,        14.6128,        16.3158,        18.1115,        20.0000};
593  static const G4double I6[nX]={0,
594        .053692,        .202354,        .443946,        .778765,        1.20774,        1.73208,        2.35319,        3.07256,        3.89177,        4.81249,
595        5.83641,        6.96528,        8.20092,        9.54516,        10.9999,        12.5670,        14.2486,        16.0466,        17.9630,        20.0000};
596  static const G4double I7[nX]={0,
597        .043065,        .168099,        .376879,        .672273,        1.05738,        1.53543,        2.10973,        2.78364,        3.56065,        4.44429,
598        5.43819,        6.54610,        7.77186,        9.11940,        10.5928,        12.1963,        13.9342,        15.8110,        17.8313,        20.0000};
599  static const G4double I8[nX]={0,
600        .036051,        .143997,        .327877,        .592202,        .941572,        1.38068,        1.91433,        2.54746,        3.28517,        4.13277,
601        5.09574,        6.17984,        7.39106,        8.73568,        10.2203,        11.8519,        13.6377,        15.5854,        17.7033,        20.0000};
602  static const G4double I9[nX]={0,
603        .030977,        .125727,        .289605,        .528146,        .846967,        1.25183,        1.74871,        2.34384,        3.04376,        3.85535,
604        4.78594,        5.84329,        7.03567,        8.37194,        9.86163,        11.5150,        13.3430,        15.3576,        17.5719,        20.0000};
605  static const G4double IA[nX]={0,
606        .027129,        .111420,        .258935,        .475812,        .768320,        1.14297,        1.60661,        2.16648,        2.83034,        3.60650,
607        4.50394,        5.53238,        6.70244,        8.02569,        9.51488,        11.1841,        13.0488,        15.1264,        17.4362,        20.0000};
608  static const G4double IB[nX]={0,
609        .024170,        .100153,        .234345,        .433198,        .703363,        1.05184,        1.48607,        2.01409,        2.64459,        3.38708,
610        4.25198,        5.25084,        6.39647,        7.70319,        9.18708,        10.8663,        12.7617,        14.8968,        17.2990,        20.0000};
611  static const G4double IC[nX]={0,
612        .021877,        .091263,        .214670,        .398677,        .650133,        .976322,        1.38510,        1.88504,        2.48555,        3.19709,
613        4.03129,        5.00127,        6.12184,        7.40989,        8.88482,        10.5690,        12.4888,        14.6748,        17.1638,        20.0000};
614  static const G4double ID[nX]={0,
615        .020062,        .084127,        .198702,        .370384,        .606100,        .913288,        1.30006,        1.77535,        2.34912,        3.03253,
616        3.83822,        4.78063,        5.87634,        7.14459,        8.60791,        10.2929,        12.2315,        14.4621,        17.0320,        20.0000};
617  static const G4double IE[nX]={0,
618        .018547,        .078104,        .185102,        .346090,        .567998,        .858331,        1.22535,        1.67824,        2.22735,        2.88443,
619        3.66294,        4.57845,        5.64911,        6.89637,        8.34578,        10.0282,        11.9812,        14.2519,        16.8993,        20.0000};
620  static const G4double IF[nX]={0,
621        .017143,        .072466,        .172271,        .323007,        .531545,        .805393,        1.15288,        1.58338,        2.10754,        2.73758,
622        3.48769,        4.37450,        5.41770,        6.64092,        8.07288,        9.74894,        11.7135,        14.0232,        16.7522,        20.0000};
623  static const G4double IG[nX]={0,
624        .015618,        .066285,        .158094,        .297316,        .490692,        .745653,        1.07053,        1.47479,        1.96931,        2.56677,
625        3.28205,        4.13289,        5.14068,        6.33158,        7.73808,        9.40133,        11.3745,        13.7279,        16.5577,        20.0000};
626  static const G4double IH[nX]={0,
627        .013702,        .058434,        .139923,        .264115,        .437466,        .667179,        .961433,        1.32965,        1.78283,        2.33399,
628        2.99871,        3.79596,        4.74916,        5.88771,        7.24937,        8.88367,        10.8576,        13.2646,        16.2417,        20.0000};
629  static const G4double II[nX]={0,
630        .011264,        .048311,        .116235,        .220381,        .366634,        .561656,        .813132,        1.13008,        1.52322,        2.00554,
631        2.59296,        3.30542,        4.16834,        5.21490,        6.48964,        8.05434,        9.99835,        12.4580,        15.6567,        20.0000};
632  static const G4double IJ[nX]={0,
633        .008628,        .037206,        .089928,        .171242,        .286114,        .440251,        .640343,        .894382,        1.21208,        1.60544,
634        2.08962,        2.68414,        3.41486,        4.31700,        5.44048,        6.85936,        8.69067,        11.1358,        14.5885,        20.0000};
635  static const G4double* Il[nE]=
636                             {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
637  static const G4double lE[nE]={
638-1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215,  .459141,        .867068,        1.27499,        1.68292,
639 2.09085,       2.49877,        2.90670,        3.31463,        3.72255,        4.13048,        4.53840,        4.94633,        5.35426,        5.76218};
640  static const G4double lEmi=lE[0];
641  static const G4double lEma=lE[nE-1];
642  static const G4double dlE=(lEma-lEmi)/bE;
643         //***************************************************************************************
644  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
645  G4double lEn=std::log(Enu);         // log(E) for interpolation
646  G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
647  G4int fE=static_cast<int>(rE);      // Left bin for interpolation
648  if(fE<0) fE=0;
649  if(fE>pE)fE=pE;
650  G4int    sE=fE+1;                   // Right bin for interpolation
651  G4double dE=rE-fE;                  // relative log shift from the left bin
652  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
653  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
654  G4double Emu=Enu;                   // Free Energy of neutrino/anti-neutrino
655  G4double ME=Enu*MN;                 // M*E
656  G4double dME=ME+ME;                 // 2*M*E
657  G4double dEMN=(dEnu+MN)*ME;
658  G4double sqE=Enu*ME;
659  G4double E2M=MN*Enu2;
660  G4double ymax=(E2M+sqE)/dEMN;
661  G4double Q2mi=0.;                   // Q2_min(E_nu)
662  G4double Q2ma=dME*ymax;             // Q2_max(E_nu)
663  G4double Q2nq=Emu*dMN-mcV;
664  if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
665  // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
666  G4double Rmi=Q2mi/Q2ma;
667  G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78);
668  // --- E-interpolation must be done in a log scale ---
669  G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
670  G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
671  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
672  G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
673  G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
674  G4double rXi=(Xmi-iXmi)/idX;
675  G4int    iXi=static_cast<int>(rXi);
676  if(iXi<0) iXi=0;
677  if(iXi>bX) iXi=bX;
678  G4double dXi=rXi-iXi;
679  G4double bntil=Il[fE][iXi];
680  G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
681  G4double bntir=Il[sE][iXi];
682  G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
683  G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
684  //
685  G4double rXa=(Xma-iXmi)/idX;
686  G4int    iXa=static_cast<int>(rXa);
687  if(iXa<0) iXa=0;
688  if(iXa>bX) iXa=bX;
689  G4double dXa=rXa-iXa;
690  G4double bntal=Il[fE][iXa];
691  G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
692  G4double bntar=Il[sE][iXa];
693  G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
694  G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
695  //
696  // *** Find X using the reversed table ***
697  G4double intx=inti+(inta-inti)*G4UniformRand(); 
698  G4int    intc=static_cast<int>(intx);
699  if(intc<0) intc=0;
700  if(intc>bX) intc=bX;
701  G4double dint=intx-intc;
702  G4double mXl=Xl[fE][intc];
703  G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
704  G4double mXr=Xl[sE][intc];
705  G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
706  G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
707  G4double R=shift-std::pow(X,pconv);
708  G4double Q2=R*Q2ma;
709  return Q2*GeV*GeV;
710}
711
712// It returns a fraction of the direct interaction of the neutrino with quark-partons
713G4double G4QANuANuNuclearCrossSection::GetDirectPart(G4double Q2)
714{
715  G4double f=Q2/4.62;
716  G4double ff=f*f;
717  G4double r=ff*ff;
718  G4double s=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
719  //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
720  return 1.-s*(1.-s/2);
721}
722
723// #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
724G4double G4QANuANuNuclearCrossSection::GetNPartons(G4double Q2)
725{
726  return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
727}
728
729// This class can provide only virtual exchange pi+ (a substitute for W+ boson)
730G4int G4QANuANuNuclearCrossSection::GetExchangePDGCode() {return 22;}
Note: See TracBrowser for help on using the repository browser.