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

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

update ti head

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