G4QANuANuNuclearCrossSection.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 //
00030 // G4 Physics class: G4QANuANuNuclearCrossSection for gamma+A cross sections
00031 // Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
00032 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 17-Oct-03
00033 // 
00034 // ****************************************************************************************
00035 // ********** This CLASS is temporary moved from the photolepton_hadron directory *********
00036 // ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ******
00037 // ****************************************************************************************
00038 // ----------------------------------------------------------------------------------------
00039 // Short description: Provides the (anti-nu,anti-nu)A cross-section (Created by M. Kossov)
00040 // ----------------------------------------------------------------------------------------
00041 
00042 //#define debug
00043 //#define edebug
00044 //#define pdebug
00045 //#define ppdebug
00046 //#define tdebug
00047 //#define sdebug
00048 
00049 #include "G4QANuANuNuclearCrossSection.hh"
00050 #include "G4SystemOfUnits.hh"
00051 
00052 // Initialization of the
00053 G4bool    G4QANuANuNuclearCrossSection::onlyCS=true;//Flag to calculate only CS (not QE)
00054 G4double  G4QANuANuNuclearCrossSection::lastSig=0.;//Last calculated total cross section
00055 G4double  G4QANuANuNuclearCrossSection::lastQEL=0.;//Last calculated quasi-el cross section
00056 G4int     G4QANuANuNuclearCrossSection::lastL=0;   //Last used in cross section TheLastBin
00057 G4double  G4QANuANuNuclearCrossSection::lastE=0.;  //Last used in cross section TheEnergy
00058 G4double* G4QANuANuNuclearCrossSection::lastEN=0;  //Pointer to theEnergy Scale of TX & QE
00059 G4double* G4QANuANuNuclearCrossSection::lastTX=0;  //Pointer to theLastArray of TX function
00060 G4double* G4QANuANuNuclearCrossSection::lastQE=0;  //Pointer to theLastArray of QE function
00061 G4int     G4QANuANuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
00062 G4int     G4QANuANuNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
00063 G4int     G4QANuANuNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
00064 G4double  G4QANuANuNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
00065 G4double  G4QANuANuNuclearCrossSection::lastTH=0.; // Last threshold momentum
00066 G4double  G4QANuANuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
00067 G4int     G4QANuANuNuclearCrossSection::lastI=0;   // The last position in the DAMDB
00068 std::vector<G4double*>* G4QANuANuNuclearCrossSection::TX = new std::vector<G4double*>;
00069 std::vector<G4double*>* G4QANuANuNuclearCrossSection::QE = new std::vector<G4double*>;
00070 
00071 // Returns Pointer to the G4VQCrossSection class
00072 G4VQCrossSection* G4QANuANuNuclearCrossSection::GetPointer()
00073 {
00074   static G4QANuANuNuclearCrossSection theCrossSection; //*Static body of the Cross Section*
00075   return &theCrossSection;
00076 }
00077 
00078 G4QANuANuNuclearCrossSection::~G4QANuANuNuclearCrossSection()
00079 {
00080   G4int lens=TX->size();
00081   for(G4int i=0; i<lens; ++i) delete[] (*TX)[i];
00082   delete TX;
00083   G4int hens=QE->size();
00084   for(G4int i=0; i<hens; ++i) delete[] (*QE)[i];
00085   delete QE;
00086 }
00087 
00088 // The main member function giving the collision cross section (P is in IU, CS is in mb)
00089 // Make pMom in independent units ! (Now it is MeV)
00090 G4double G4QANuANuNuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
00091                                                       G4int tgZ, G4int tgN, G4int pPDG)
00092 {
00093   static G4int j;                      // A#0f records found in DB for this projectile
00094   static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
00095   static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
00096   static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
00097   static std::vector <G4double> colP;  // Vector of last momenta for the reaction
00098   static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
00099   static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
00100   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
00101   G4double pEn=pMom;
00102 #ifdef debug
00103   G4cout<<"G4QAMNCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
00104         <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
00105         <<colN.size()<<G4endl;
00106   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00107 #endif
00108   if(pPDG!=-14)
00109   {
00110 #ifdef debug
00111     G4cout<<"G4QAMNCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00112     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00113 #endif
00114     return 0.;                         // projectile PDG=0 is a mistake (?!) @@
00115   }
00116   G4bool in=false;                     // By default the isotope must be found in the AMDB
00117   if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
00118   {
00119     in = false;                        // By default the isotope haven't be found in AMDB  
00120     lastP   = 0.;                      // New momentum history (nothing to compare with)
00121     lastPDG = pPDG;                    // The last PDG of the projectile
00122     lastN   = tgN;                     // The last N of the calculated nucleus
00123     lastZ   = tgZ;                     // The last Z of the calculated nucleus
00124     lastI   = colN.size();             // Size of the Associative Memory DB in the heap
00125     j  = 0;                            // A#0f records found in DB for this projectile
00126     if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found
00127     {                                  // The nucleus with projPDG is found in AMDB
00128       if(colN[i]==tgN && colZ[i]==tgZ)
00129       {
00130         lastI=i;
00131         lastTH =colTH[i];                // Last THreshold (A-dependent)
00132 #ifdef pdebug
00133         G4cout<<"G4QAMNCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00134         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00135 #endif
00136         if(pEn<=lastTH)
00137         {
00138 #ifdef pdebug
00139           G4cout<<"G4QAMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl;
00140           //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00141 #endif
00142           return 0.;                     // Energy is below the Threshold value
00143         }
00144         lastP  =colP [i];                // Last Momentum  (A-dependent)
00145         lastCS =colCS[i];                // Last CrossSect (A-dependent)
00146         if(std::fabs(lastP/pMom-1.)<tolerance)
00147         {
00148 #ifdef pdebug
00149           G4cout<<"G4QAMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00150           //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00151 #endif
00152           return lastCS*millibarn;     // Use theLastCS
00153         }
00154         in = true;                       // This is the case when the isotop is found in DB
00155         // Momentum pMom is in IU ! @@ Units
00156 #ifdef pdebug
00157         G4cout<<"G4QAMNCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00158 #endif
00159         lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
00160 #ifdef pdebug
00161         G4cout<<"G4QAMNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00162         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00163 #endif
00164         if(lastCS<=0. && pEn>lastTH)    // Correct the threshold
00165         {
00166 #ifdef pdebug
00167           G4cout<<"G4QAMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00168 #endif
00169           lastTH=pEn;
00170         }
00171         break;                           // Go out of the LOOP
00172       }
00173 #ifdef pdebug
00174       G4cout<<"---G4QAMNCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
00175             <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl;
00176       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00177 #endif
00178       j++;                             // Increment a#0f records found in DB for this pPDG
00179     }
00180     if(!in)                            // This nucleus has not been calculated previously
00181     {
00182 #ifdef pdebug
00183       G4cout<<"G4QAMNCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<lastI<<G4endl;
00184 #endif
00186       lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
00187       if(lastCS<=0.)
00188       {
00189         lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
00190 #ifdef pdebug
00191         G4cout<<"G4QAMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
00192 #endif
00193         if(pEn>lastTH)
00194         {
00195 #ifdef pdebug
00196           G4cout<<"G4QAMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00197 #endif
00198           lastTH=pEn;
00199         }
00200       }
00201 #ifdef pdebug
00202       G4cout<<"G4QAMNCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00203       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00204 #endif
00205       colN.push_back(tgN);
00206       colZ.push_back(tgZ);
00207       colPDG.push_back(pPDG);
00208       colP.push_back(pMom);
00209       colTH.push_back(lastTH);
00210       colCS.push_back(lastCS);
00211 #ifdef pdebug
00212       G4cout<<"G4QAMNCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl;
00213       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00214 #endif
00215       return lastCS*millibarn;
00216     } // End of creation of the new set of parameters
00217     else
00218     {
00219 #ifdef pdebug
00220       G4cout<<"G4QAMNCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00221 #endif
00222       colP[lastI]=pMom;
00223       colPDG[lastI]=pPDG;
00224       colCS[lastI]=lastCS;
00225     }
00226   } // End of parameters udate
00227   else if(pEn<=lastTH)
00228   {
00229 #ifdef pdebug
00230     G4cout<<"G4QAMNCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00231     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00232 #endif
00233     return 0.;                         // Momentum is below the Threshold Value -> CS=0
00234   }
00235   else if(std::fabs(lastP/pMom-1.)<tolerance)
00236   {
00237 #ifdef pdebug
00238     G4cout<<"G4QAMNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00239     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00240 #endif
00241     return lastCS*millibarn;     // Use theLastCS
00242   }
00243   else
00244   {
00245 #ifdef pdebug
00246     G4cout<<"G4QAMNCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00247 #endif
00248     lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00249     lastP=pMom;
00250   }
00251 #ifdef pdebug
00252   G4cout<<"G4QAMNCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00253   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00254 #endif
00255   return lastCS*millibarn;
00256 }
00257 
00258 // Gives the threshold energy = the same for all nuclei (@@ can be reduced for hevy nuclei)
00259 G4double G4QANuANuNuclearCrossSection::ThresholdEnergy(G4int, G4int, G4int) {return 0.;}
00260 
00261 // The main member function giving the gamma-A cross section (E_kin in MeV, CS in mb)
00262 G4double G4QANuANuNuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I,
00263                                        G4int , G4int targZ, G4int targN, G4double Momentum)
00264 {
00265   static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2
00266   static const G4int nE=65;   // !! If change this, change it in GetFunctions() (*.hh) !!
00267   static const G4int mL=nE-1;
00268   static const G4double EMi=0.;      // Universal threshold of the reaction in GeV
00269   static const G4double EMa=300.;    // Maximum tabulated Energy of nu_mu in GeV 
00270   // *** Begin of the Associative memory for acceleration of the cross section calculations
00271   static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional)
00272   static G4bool first=true;          // Flag of initialization of the energy axis
00273   // *** End of Static Definitions (Associative Memory) ***
00274   //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
00275   //G4double TotEnergy2=Momentum;
00276   onlyCS=CS;                         // Flag to calculate only CS (not TX & QE)
00277   lastE=Momentum/GeV;                // Kinetic energy of the muon neutrino (in GeV!)
00278   if (lastE<=EMi)                    // Energy is below the minimum energy in the table
00279   {
00280     lastE=0.;
00281     lastSig=0.;
00282     return 0.;
00283   }
00284   G4int Z=targZ;                     // New Z, which can change the sign
00285   if(F<=0)                           // This isotope was not the last used isotop
00286   {
00287     if(F<0)                          // This isotope was found in DAMDB =-------=> RETRIEVE
00288     {
00289       lastTX =(*TX)[I];              // Pointer to the prepared TX function (same isotope)
00290       lastQE =(*QE)[I];              // Pointer to the prepared QE function (same isotope)
00291    }
00292    else                              // This isotope wasn't calculated previously => CREATE
00293    {
00294       if(first)
00295       {
00296         lastEN = new G4double[nE];   // This must be done only once!
00297         Z=-Z;                        // To explain GetFunctions that E-axis must be filled
00298         first=false;                 // To make it only once
00299       }
00300       lastTX = new G4double[nE];     // Allocate memory for the new TX function
00301       lastQE = new G4double[nE];     // Allocate memory for the new QE function
00302       G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
00303       if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
00304       // *** The synchronization check ***
00305       G4int sync=TX->size();
00306       if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
00307       TX->push_back(lastTX);
00308       QE->push_back(lastQE);
00309     } // End of creation of the new set of parameters
00310   } // End of parameters udate
00311   // =-----------------= NOW Calculate the Cross Section =---------------=
00312   if (lastE<=EMi)                    // Check that antiNuEnergy is higher than ThreshE
00313   {
00314     lastE=0.;
00315     lastSig=0.;
00316     return 0.;
00317   }
00318   if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
00319   {
00320     G4int chk=1;
00321     G4int ran=mL/2;
00322     G4int sep=ran;  // as a result = an index of the left edge of the interval
00323     while(ran>=2)
00324     {
00325       G4int newran=ran/2;
00326       if(lastE<=lastEN[sep]) sep-=newran;
00327       else                   sep+=newran;
00328       ran=newran;
00329       chk=chk+chk; 
00330     }
00331     if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
00332     G4double lowE=lastEN[sep];
00333     G4double highE=lastEN[sep+1];
00334     G4double lowTX=lastTX[sep];
00335     if(lastE<lowE||sep>=mL||lastE>highE)
00336       G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
00337             <<", sep="<<sep<<", mL="<<mL<<G4endl;
00338     lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
00339     if(!onlyCS)                       // Skip the differential cross-section parameters
00340     {
00341       G4double lowQE=lastQE[sep];
00342       lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
00343 #ifdef pdebug
00344       G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
00345 #endif
00346     }
00347   }
00348   else
00349   {
00350     lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
00351     lastQEL=lastQE[mL];
00352   }
00353   if(lastQEL<0.) lastQEL = 0.;
00354   if(lastSig<0.) lastSig = 0.;
00355   // The cross-sections are expected to be in mb
00356   lastSig*=mb38;
00357   if(!onlyCS) lastQEL*=mb38;
00358   return lastSig;
00359 }
00360 
00361 // Calculate the cros-section functions
00362 // ****************************************************************************************
00363 // *** This tables are the same for all lepto-nuclear reactions, only mass is different ***
00364 // ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ ***
00365 // ****************************************************************************************
00366 G4int G4QANuANuNuclearCrossSection::GetFunctions (G4int z, G4int n,
00367                                                      G4double* t, G4double* q, G4double* e)
00368 {
00369   static const G4int nE=65; // !! If change this, change it in GetCrossSection() (*.cc) !!
00370   static const G4double nuEn[nE]={0.,
00371   1.00463e-5,1.05336e-5,1.10692e-5,1.16592e-5,1.23109e-5,1.30323e-5,1.38331e-5,1.47245e-5,
00372   1.57194e-5,1.68335e-5,1.80848e-5,1.94948e-5,2.10894e-5,2.28991e-5,2.49608e-5,2.73189e-5,
00373   3.00273e-5,3.31516e-5,3.67722e-5,4.09881e-5,4.59217e-5,5.17255e-5,5.85908e-5,6.67583e-5,
00374   7.65338e-5,8.83078e-5,.000102583,.000120011,.000141441,.000167995,.000201160,.000242926,
00375   .000295985,.000364008,.000452051,.000567152,.000719210,.000922307,.001196710,.001571930,
00376   .002091530,.002820590,.003857810,.005354930,.007548840,.010815300,.015760100,.023376900,
00377   .035325600,.054430800,.085595700,.137508000,.225898000,.379892000,.654712000,1.15767000,
00378   2.10277000,3.92843000,7.55861000,14.9991000,30.7412000,65.1734000,143.155000,326.326000};
00379   static const G4double TOTX[nE]={0.,
00380   3.63538e-5,3.81165e-5,4.00539e-5,4.21884e-5,4.45454e-5,4.71548e-5,5.00511e-5,5.32747e-5,
00381   5.68730e-5,6.09016e-5,6.54261e-5,7.05246e-5,7.62894e-5,8.28315e-5,9.02838e-5,9.88066e-5,
00382   .000108594,.000119884,.000132964,.000148191,.000166007,.000186959,.000211736,.000241202,
00383   .000276453,.000318890,.000370311,.000433042,.000510116,.000605516,.000724514,.000874144,
00384   .001063870,.001306520,.001619660,.002027520,.002563780,.003275710,.004230080,.005521890,
00385   .007286920,.009719890,.013099700,.018695100,.029208400,.042095000,.059253700,.082373900,
00386   .113071000,.151041000,.191803000,.224208000,.234187000,.217774000,.187139000,.157818000,
00387   .137494000,.125872000,.120462000,.119148000,.120418000,.123027000,.126408000,.129071000};
00388   static const G4double QELX[nE]={0.,
00389   .365220e-9,.401502e-9,.443364e-9,.491885e-9,.548393e-9,.614536e-9,.692362e-9,.784441e-9,
00390   .894012e-9,1.02519e-9,1.18322e-9,1.37487e-9,1.60890e-9,1.89677e-9,2.25355e-9,2.69928e-9,
00391   3.26079e-9,3.97433e-9,4.88937e-9,6.07407e-9,7.62331e-9,9.67058e-9,1.24058e-8,1.61022e-8,
00392   2.11580e-8,2.81605e-8,3.79876e-8,5.19696e-8,7.21515e-8,1.01724e-7,1.45743e-7,2.12353e-7,
00393   3.14890e-7,4.75585e-7,7.32171e-7,1.14991e-6,1.84390e-6,3.02121e-6,5.06217e-6,8.68002e-6,
00394   1.52408e-5,2.74159e-5,5.05363e-5,.000100111,.000220489,.000455269,.000933841,.001925650,
00395   .003994300,.008221270,.016417600,.030830400,.052902400,.082519200,.115560000,.149598000,
00396   .184112000,.215102000,.238253000,.252949000,.261267000,.265626000,.267782000,.268791000};
00397   // --------------------------------
00398   G4int first=0;
00399   if(z<0.)
00400   {
00401     first=1;
00402     z=-z;
00403   }
00404   if(z<1 || z>92)             // neutron & plutonium are forbidden
00405   {
00406     G4cout<<"*G4QANuANuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
00407     return -1;
00408   }
00409   for(G4int k=0; k<nE; k++)
00410   {
00411     G4double a=n+z;
00412     G4double za=z+a;
00413     G4double dz=z+z;
00414     G4double da=a+a;
00415     G4double ta=da+a;
00416     if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
00417     t[k]=TOTX[k]*nuEn[k]*(za+za)/ta+QELX[k]*(dz+dz-da)/ta; // TotalCrossSection
00418     q[k]=QELX[k]*dz/a;                                     // QuasiElasticCrossSection
00419   }
00420   return first;
00421 }
00422 
00423 // Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
00424 G4double G4QANuANuNuclearCrossSection::GetQEL_ExchangeQ2()
00425 {
00426   static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
00427   static const G4double power=-3.5;    // direct power for the magic variable
00428   static const G4double pconv=1./power;// conversion power for the magic variable
00429   static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
00430   static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
00431   static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
00432   // Reversed table
00433   static const G4double Xl[nQ2]={5.20224e-16,
00434  .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898,
00435  .093878, .103332, .112844, .122410, .132023, .141680, .151376, .161109, .170875, .180672,
00436  .190499, .200352, .210230, .220131, .230055, .239999, .249963, .259945, .269944, .279960,
00437  .289992, .300039, .310099, .320173, .330260, .340359, .350470, .360592, .370724, .380867,
00438  .391019, .401181, .411352, .421531, .431719, .441915, .452118, .462329, .472547, .482771,
00439  .493003, .503240, .513484, .523734, .533989, .544250, .554517, .564788, .575065, .585346,
00440  .595632, .605923, .616218, .626517, .636820, .647127, .657438, .667753, .678072, .688394,
00441  .698719, .709048, .719380, .729715, .740053, .750394, .760738, .771085, .781434, .791786,
00442  .802140, .812497, .822857, .833219, .843582, .853949, .864317, .874687, .885060, .895434,
00443  .905810, .916188, .926568, .936950, .947333, .957719, .968105, .978493, .988883, .999275};
00444    // Direct table
00445   static const G4double Xmax=Xl[lQ2];
00446   static const G4double Xmin=Xl[0];
00447   static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
00448   static const G4double inl[nQ2]={0,
00449  1.52225, 2.77846, 3.96651, 5.11612, 6.23990, 7.34467, 8.43466, 9.51272, 10.5809, 11.6406,
00450  12.6932, 13.7394, 14.7801, 15.8158, 16.8471, 17.8743, 18.8979, 19.9181, 20.9353, 21.9496,
00451  22.9614, 23.9707, 24.9777, 25.9826, 26.9855, 27.9866, 28.9860, 29.9837, 30.9798, 31.9745,
00452  32.9678, 33.9598, 34.9505, 35.9400, 36.9284, 37.9158, 38.9021, 39.8874, 40.8718, 41.8553,
00453  42.8379, 43.8197, 44.8007, 45.7810, 46.7605, 47.7393, 48.7174, 49.6950, 50.6718, 51.6481,
00454  52.6238, 53.5990, 54.5736, 55.5476, 56.5212, 57.4943, 58.4670, 59.4391, 60.4109, 61.3822,
00455  62.3531, 63.3236, 64.2937, 65.2635, 66.2329, 67.2019, 68.1707, 69.1390, 70.1071, 71.0748,
00456  72.0423, 73.0095, 73.9763, 74.9429, 75.9093, 76.8754, 77.8412, 78.8068, 79.7721, 80.7373,
00457  81.7022, 82.6668, 83.6313, 84.5956, 85.5596, 86.5235, 87.4872, 88.4507, 89.4140, 90.3771,
00458  91.3401, 92.3029, 93.2656, 94.2281, 95.1904, 96.1526, 97.1147, 98.0766, 99.0384, 100.000};
00459   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00460   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00461   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00462   G4double ME=Enu*MN;                 // M*E
00463   G4double dME=ME+ME;                 // 2*M*E
00464   G4double dEMN=(dEnu+MN)*ME;
00465   G4double sqE=Enu*ME;
00466   G4double E2M=MN*Enu2;
00467   G4double ymax=(E2M+sqE)/dEMN;
00468   G4double Q2mi=0.; // Q2_min(E_nu)
00469   G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
00470   G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
00471   G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
00472   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00473   G4double rXi=(Xmi-Xmin)/dX;
00474   G4int    iXi=static_cast<int>(rXi);
00475   if(iXi<0) iXi=0;
00476   if(iXi>bQ2) iXi=bQ2;
00477   G4double dXi=rXi-iXi;
00478   G4double bnti=inl[iXi];
00479   G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
00480   //
00481   G4double rXa=(Xma-Xmin)/dX;
00482   G4int    iXa=static_cast<int>(rXa);
00483   if(iXa<0) iXa=0;
00484   if(iXa>bQ2) iXa=bQ2;
00485   G4double dXa=rXa-iXa;
00486   G4double bnta=inl[iXa];
00487   G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
00488   // *** Find X using the reversed table ***
00489   G4double intx=inti+(inta-inti)*G4UniformRand();
00490   G4int    intc=static_cast<int>(intx);
00491   if(intc<0) intc=0;
00492   if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
00493   G4double dint=intx-intc;
00494   G4double mX=Xl[intc];
00495   G4double X=mX+dint*(Xl[intc+1]-mX);
00496   G4double Q2=std::pow(X,pconv)-1.;
00497   return Q2*GeV*GeV;
00498 }
00499 
00500 // Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
00501 G4double G4QANuANuNuclearCrossSection::GetNQE_ExchangeQ2()
00502 {
00503   static const double mpi=.13957018;    // charged pi meson mass in GeV
00504   static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
00505   static const double dMN=MN+MN;        // 2*M_N in GeV
00506   static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
00507   static const G4int power=7;           // direct power for the magic variable
00508   static const G4double pconv=1./power; // conversion power for the magic variable
00509   static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
00510   static const G4int lX=nX-1;           // index of the last in the Xl table
00511   static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
00512   static const G4int nE=20;             // #Of point in the El table (in GeV^2)
00513   static const G4int bE=nE-1;           // index of the last in the El table
00514   static const G4int pE=bE-1;           // index of the before last in the El table
00515   // Reversed table
00516   static const G4double X0[nX]={5.21412e-05,
00517  .437860, .681908, .891529, 1.08434, 1.26751, 1.44494, 1.61915, 1.79198, 1.96493, 2.13937,
00518  2.31664, 2.49816, 2.68559, 2.88097, 3.08705, 3.30774, 3.54917, 3.82233, 4.15131, 4.62182};
00519   static const G4double X1[nX]={.00102591,
00520  1.00443, 1.55828, 2.03126, 2.46406, 2.87311, 3.26723, 3.65199, 4.03134, 4.40835, 4.78561,
00521  5.16549, 5.55031, 5.94252, 6.34484, 6.76049, 7.19349, 7.64917, 8.13502, 8.66246, 9.25086};
00522   static const G4double X2[nX]={.0120304,
00523  2.59903, 3.98637, 5.15131, 6.20159, 7.18024, 8.10986, 9.00426, 9.87265, 10.7217, 11.5564,
00524  12.3808, 13.1983, 14.0116, 14.8234, 15.6359, 16.4515, 17.2723, 18.1006, 18.9386, 19.7892};
00525   static const G4double X3[nX]={.060124,
00526  5.73857, 8.62595, 10.9849, 13.0644, 14.9636, 16.7340, 18.4066, 20.0019, 21.5342, 23.0142,
00527  24.4497, 25.8471, 27.2114, 28.5467, 29.8564, 31.1434, 32.4102, 33.6589, 34.8912, 36.1095};
00528   static const G4double X4[nX]={.0992363,
00529  8.23746, 12.1036, 15.1740, 17.8231, 20.1992, 22.3792, 24.4092, 26.3198, 28.1320, 29.8615,
00530  31.5200, 33.1169, 34.6594, 36.1536, 37.6044, 39.0160, 40.3920, 41.7353, 43.0485, 44.3354};
00531   static const G4double X5[nX]={.0561127,
00532  7.33661, 10.5694, 13.0778, 15.2061, 17.0893, 18.7973, 20.3717, 21.8400, 23.2211, 24.5291,
00533  25.7745, 26.9655, 28.1087, 29.2094, 30.2721, 31.3003, 32.2972, 33.2656, 34.2076, 35.1265};
00534   static const G4double X6[nX]={.0145859,
00535  4.81774, 6.83565, 8.37399, 9.66291, 10.7920, 11.8075, 12.7366, 13.5975, 14.4025, 15.1608,
00536  15.8791, 16.5628, 17.2162, 17.8427, 18.4451, 19.0259, 19.5869, 20.1300, 20.6566, 21.1706};
00537   static const G4double X7[nX]={.00241155,
00538  2.87095, 4.02492, 4.89243, 5.61207, 6.23747, 6.79613, 7.30433, 7.77270, 8.20858, 8.61732,
00539  9.00296, 9.36863, 9.71682, 10.0495, 10.3684, 10.6749, 10.9701, 11.2550, 11.5306, 11.7982};
00540   static const G4double X8[nX]={.000316863,
00541  1.76189, 2.44632, 2.95477, 3.37292, 3.73378, 4.05420, 4.34415, 4.61009, 4.85651, 5.08666,
00542  5.30299, 5.50738, 5.70134, 5.88609, 6.06262, 6.23178, 6.39425, 6.55065, 6.70149, 6.84742};
00543   static const G4double X9[nX]={3.73544e-05,
00544  1.17106, 1.61289, 1.93763, 2.20259, 2.42976, 2.63034, 2.81094, 2.97582, 3.12796, 3.26949,
00545  3.40202, 3.52680, 3.64482, 3.75687, 3.86360, 3.96557, 4.06323, 4.15697, 4.24713, 4.33413};
00546   static const G4double XA[nX]={4.19131e-06,
00547  .849573, 1.16208, 1.38955, 1.57379, 1.73079, 1.86867, 1.99221, 2.10451, 2.20770, 2.30332,
00548  2.39252, 2.47622, 2.55511, 2.62977, 2.70066, 2.76818, 2.83265, 2.89437, 2.95355, 3.01051};
00549   static const G4double XB[nX]={4.59981e-07,
00550  .666131, .905836, 1.07880, 1.21796, 1.33587, 1.43890, 1.53080, 1.61399, 1.69011, 1.76040,
00551  1.82573, 1.88682, 1.94421, 1.99834, 2.04959, 2.09824, 2.14457, 2.18878, 2.23107, 2.27162};
00552   static const G4double XC[nX]={4.99861e-08,
00553  .556280, .752730, .893387, 1.00587, 1.10070, 1.18317, 1.25643, 1.32247, 1.38269, 1.43809,
00554  1.48941, 1.53724, 1.58203, 1.62416, 1.66391, 1.70155, 1.73728, 1.77128, 1.80371, 1.83473};
00555   static const G4double XD[nX]={5.40832e-09,
00556  .488069, .657650, .778236, .874148, .954621, 1.02432, 1.08599, 1.14138, 1.19172, 1.23787,
00557  1.28049, 1.32008, 1.35705, 1.39172, 1.42434, 1.45514, 1.48429, 1.51197, 1.53829, 1.56339};
00558   static const G4double XE[nX]={5.84029e-10,
00559  .445057, .597434, .705099, .790298, .861468, .922865, .976982, 1.02542, 1.06930, 1.10939,
00560  1.14630, 1.18050, 1.21233, 1.24208, 1.27001, 1.29630, 1.32113, 1.34462, 1.36691, 1.38812};
00561   static const G4double XF[nX]={6.30137e-11,
00562  .418735, .560003, .659168, .737230, .802138, .857898, .906854, .950515, .989915, 1.02580,
00563  1.05873, 1.08913, 1.11734, 1.14364, 1.16824, 1.19133, 1.21306, 1.23358, 1.25298, 1.27139};
00564   static const G4double XG[nX]={6.79627e-12,
00565  .405286, .539651, .633227, .706417, .766929, .818642, .863824, .903931, .939963, .972639,
00566  1.00250, 1.02995, 1.05532, 1.07887, 1.10082, 1.12134, 1.14058, 1.15867, 1.17572, 1.19183};
00567   static const G4double XH[nX]={7.32882e-13,
00568  .404391, .535199, .625259, .695036, .752243, .800752, .842823, .879906, .912994, .942802,
00569  .969862, .994583, 1.01729, 1.03823, 1.05763, 1.07566, 1.09246, 1.10816, 1.12286, 1.13667};
00570   static const G4double XI[nX]={7.90251e-14,
00571  .418084, .548382, .636489, .703728, .758106, .803630, .842633, .876608, .906576, .933269,
00572  .957233, .978886, .998556, 1.01651, 1.03295, 1.04807, 1.06201, 1.07489, 1.08683, 1.09792};
00573   static const G4double XJ[nX]={8.52083e-15,
00574  .447299, .579635, .666780, .731788, .783268, .825512, .861013, .891356, .917626, .940597,
00575  .960842, .978802, .994820, 1.00917, 1.02208, 1.03373, 1.04427, 1.05383, 1.06253, 1.07046};
00576   // Direct table
00577   static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
00578                         X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
00579   static const G4double dX[nE]={
00580     (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
00581     (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
00582     (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
00583     (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
00584     (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
00585   static const G4double* Xl[nE]=
00586                              {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
00587   static const G4double I0[nX]={0,
00588  .354631, 1.08972, 2.05138, 3.16564, 4.38343, 5.66828, 6.99127, 8.32858, 9.65998, 10.9680,
00589  12.2371, 13.4536, 14.6050, 15.6802, 16.6686, 17.5609, 18.3482, 19.0221, 19.5752, 20.0000};
00590   static const G4double I1[nX]={0,
00591  .281625, .877354, 1.67084, 2.60566, 3.64420, 4.75838, 5.92589, 7.12829, 8.34989, 9.57708,
00592  10.7978, 12.0014, 13.1781, 14.3190, 15.4162, 16.4620, 17.4496, 18.3724, 19.2245, 20.0000};
00593   static const G4double I2[nX]={0,
00594  .201909, .642991, 1.24946, 1.98463, 2.82370, 3.74802, 4.74263, 5.79509, 6.89474, 8.03228,
00595  9.19947, 10.3889, 11.5938, 12.8082, 14.0262, 15.2427, 16.4527, 17.6518, 18.8356, 20.0000};
00596   static const G4double I3[nX]={0,
00597  .140937, .461189, .920216, 1.49706, 2.17728, 2.94985, 3.80580, 4.73758, 5.73867, 6.80331,
00598  7.92637, 9.10316, 10.3294, 11.6013, 12.9150, 14.2672, 15.6548, 17.0746, 18.5239, 20.0000};
00599   static const G4double I4[nX]={0,
00600  .099161, .337358, .694560, 1.16037, 1.72761, 2.39078, 3.14540, 3.98768, 4.91433, 5.92245,
00601  7.00942, 8.17287, 9.41060, 10.7206, 12.1010, 13.5500, 15.0659, 16.6472, 18.2924, 20.0000};
00602   static const G4double I5[nX]={0,
00603  .071131, .255084, .543312, .932025, 1.41892, 2.00243, 2.68144, 3.45512, 4.32283, 5.28411,
00604  6.33859, 7.48602, 8.72621, 10.0590, 11.4844, 13.0023, 14.6128, 16.3158, 18.1115, 20.0000};
00605   static const G4double I6[nX]={0,
00606  .053692, .202354, .443946, .778765, 1.20774, 1.73208, 2.35319, 3.07256, 3.89177, 4.81249,
00607  5.83641, 6.96528, 8.20092, 9.54516, 10.9999, 12.5670, 14.2486, 16.0466, 17.9630, 20.0000};
00608   static const G4double I7[nX]={0,
00609  .043065, .168099, .376879, .672273, 1.05738, 1.53543, 2.10973, 2.78364, 3.56065, 4.44429,
00610  5.43819, 6.54610, 7.77186, 9.11940, 10.5928, 12.1963, 13.9342, 15.8110, 17.8313, 20.0000};
00611   static const G4double I8[nX]={0,
00612  .036051, .143997, .327877, .592202, .941572, 1.38068, 1.91433, 2.54746, 3.28517, 4.13277,
00613  5.09574, 6.17984, 7.39106, 8.73568, 10.2203, 11.8519, 13.6377, 15.5854, 17.7033, 20.0000};
00614   static const G4double I9[nX]={0,
00615  .030977, .125727, .289605, .528146, .846967, 1.25183, 1.74871, 2.34384, 3.04376, 3.85535,
00616  4.78594, 5.84329, 7.03567, 8.37194, 9.86163, 11.5150, 13.3430, 15.3576, 17.5719, 20.0000};
00617   static const G4double IA[nX]={0,
00618  .027129, .111420, .258935, .475812, .768320, 1.14297, 1.60661, 2.16648, 2.83034, 3.60650,
00619  4.50394, 5.53238, 6.70244, 8.02569, 9.51488, 11.1841, 13.0488, 15.1264, 17.4362, 20.0000};
00620   static const G4double IB[nX]={0,
00621  .024170, .100153, .234345, .433198, .703363, 1.05184, 1.48607, 2.01409, 2.64459, 3.38708,
00622  4.25198, 5.25084, 6.39647, 7.70319, 9.18708, 10.8663, 12.7617, 14.8968, 17.2990, 20.0000};
00623   static const G4double IC[nX]={0,
00624  .021877, .091263, .214670, .398677, .650133, .976322, 1.38510, 1.88504, 2.48555, 3.19709,
00625  4.03129, 5.00127, 6.12184, 7.40989, 8.88482, 10.5690, 12.4888, 14.6748, 17.1638, 20.0000};
00626   static const G4double ID[nX]={0,
00627  .020062, .084127, .198702, .370384, .606100, .913288, 1.30006, 1.77535, 2.34912, 3.03253,
00628  3.83822, 4.78063, 5.87634, 7.14459, 8.60791, 10.2929, 12.2315, 14.4621, 17.0320, 20.0000};
00629   static const G4double IE[nX]={0,
00630  .018547, .078104, .185102, .346090, .567998, .858331, 1.22535, 1.67824, 2.22735, 2.88443,
00631  3.66294, 4.57845, 5.64911, 6.89637, 8.34578, 10.0282, 11.9812, 14.2519, 16.8993, 20.0000};
00632   static const G4double IF[nX]={0,
00633  .017143, .072466, .172271, .323007, .531545, .805393, 1.15288, 1.58338, 2.10754, 2.73758,
00634  3.48769, 4.37450, 5.41770, 6.64092, 8.07288, 9.74894, 11.7135, 14.0232, 16.7522, 20.0000};
00635   static const G4double IG[nX]={0,
00636  .015618, .066285, .158094, .297316, .490692, .745653, 1.07053, 1.47479, 1.96931, 2.56677,
00637  3.28205, 4.13289, 5.14068, 6.33158, 7.73808, 9.40133, 11.3745, 13.7279, 16.5577, 20.0000};
00638   static const G4double IH[nX]={0,
00639  .013702, .058434, .139923, .264115, .437466, .667179, .961433, 1.32965, 1.78283, 2.33399,
00640  2.99871, 3.79596, 4.74916, 5.88771, 7.24937, 8.88367, 10.8576, 13.2646, 16.2417, 20.0000};
00641   static const G4double II[nX]={0,
00642  .011264, .048311, .116235, .220381, .366634, .561656, .813132, 1.13008, 1.52322, 2.00554,
00643  2.59296, 3.30542, 4.16834, 5.21490, 6.48964, 8.05434, 9.99835, 12.4580, 15.6567, 20.0000};
00644   static const G4double IJ[nX]={0,
00645  .008628, .037206, .089928, .171242, .286114, .440251, .640343, .894382, 1.21208, 1.60544,
00646  2.08962, 2.68414, 3.41486, 4.31700, 5.44048, 6.85936, 8.69067, 11.1358, 14.5885, 20.0000};
00647   static const G4double* Il[nE]=
00648                              {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
00649   static const G4double lE[nE]={
00650 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
00651  2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
00652   static const G4double lEmi=lE[0];
00653   static const G4double lEma=lE[nE-1];
00654   static const G4double dlE=(lEma-lEmi)/bE;
00655   //***************************************************************************************
00656   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00657   G4double lEn=std::log(Enu);         // log(E) for interpolation
00658   G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
00659   G4int fE=static_cast<int>(rE);      // Left bin for interpolation
00660   if(fE<0) fE=0;
00661   if(fE>pE)fE=pE;
00662   G4int    sE=fE+1;                   // Right bin for interpolation
00663   G4double dE=rE-fE;                  // relative log shift from the left bin
00664   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00665   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00666   G4double Emu=Enu;                   // Free Energy of neutrino/anti-neutrino
00667   G4double ME=Enu*MN;                 // M*E
00668   G4double dME=ME+ME;                 // 2*M*E
00669   G4double dEMN=(dEnu+MN)*ME;
00670   G4double sqE=Enu*ME;
00671   G4double E2M=MN*Enu2;
00672   G4double ymax=(E2M+sqE)/dEMN;
00673   G4double Q2mi=0.;                   // Q2_min(E_nu)
00674   G4double Q2ma=dME*ymax;             // Q2_max(E_nu)
00675   G4double Q2nq=Emu*dMN-mcV;
00676   if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
00677   // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
00678   G4double Rmi=Q2mi/Q2ma;
00679   G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78);
00680   // --- E-interpolation must be done in a log scale ---
00681   G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
00682   G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
00683   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00684   G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
00685   G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
00686   G4double rXi=(Xmi-iXmi)/idX;
00687   G4int    iXi=static_cast<int>(rXi);
00688   if(iXi<0) iXi=0;
00689   if(iXi>bX) iXi=bX;
00690   G4double dXi=rXi-iXi;
00691   G4double bntil=Il[fE][iXi];
00692   G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
00693   G4double bntir=Il[sE][iXi];
00694   G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
00695   G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
00696   //
00697   G4double rXa=(Xma-iXmi)/idX;
00698   G4int    iXa=static_cast<int>(rXa);
00699   if(iXa<0) iXa=0;
00700   if(iXa>bX) iXa=bX;
00701   G4double dXa=rXa-iXa;
00702   G4double bntal=Il[fE][iXa];
00703   G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
00704   G4double bntar=Il[sE][iXa];
00705   G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
00706   G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
00707   //
00708   // *** Find X using the reversed table ***
00709   G4double intx=inti+(inta-inti)*G4UniformRand(); 
00710   G4int    intc=static_cast<int>(intx);
00711   if(intc<0) intc=0;
00712   if(intc>bX) intc=bX;
00713   G4double dint=intx-intc;
00714   G4double mXl=Xl[fE][intc];
00715   G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
00716   G4double mXr=Xl[sE][intc];
00717   G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
00718   G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
00719   G4double R=shift-std::pow(X,pconv);
00720   G4double Q2=R*Q2ma;
00721   return Q2*GeV*GeV;
00722 }
00723 
00724 // It returns a fraction of the direct interaction of the neutrino with quark-partons
00725 G4double G4QANuANuNuclearCrossSection::GetDirectPart(G4double Q2)
00726 {
00727   G4double f=Q2/4.62;
00728   G4double ff=f*f;
00729   G4double r=ff*ff;
00730   G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
00731   //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
00732   return 1.-s_value*(1.-s_value/2);
00733 }
00734 
00735 // #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
00736 G4double G4QANuANuNuclearCrossSection::GetNPartons(G4double Q2)
00737 {
00738   return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
00739 }
00740 
00741 // This class can provide only virtual exchange by gamma (a substitute for Z0 boson)
00742 G4int G4QANuANuNuclearCrossSection::GetExchangePDGCode() {return 22;}

Generated on Mon May 27 17:49:30 2013 for Geant4 by  doxygen 1.4.7