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

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

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