G4QNuENuclearCrossSection.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: G4QNuENuclearCrossSection for A(nu_e,e-) cross sections
00031 // Created: M.V. Kossov, CERN/ITEP(Moscow), 24-SEP-2007
00032 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 24-SEP-2007
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: nu_e -> e nuclear XS
00040 // --------------------------------------------------------------------------------------
00041 
00042 //#define debug
00043 //#define edebug
00044 //#define pdebug
00045 //#define ppdebug
00046 //#define tdebug
00047 //#define sdebug
00048 
00049 #include "G4QNuENuclearCrossSection.hh"
00050 #include "G4SystemOfUnits.hh"
00051 
00052 // Initialization of the
00053 G4bool    G4QNuENuclearCrossSection::onlyCS=true;// Flag to calculate only CS (not QE)
00054 G4double  G4QNuENuclearCrossSection::lastSig=0.;// Last calculated total cross section
00055 G4double  G4QNuENuclearCrossSection::lastQEL=0.;// Last calculated quasi-el. cross section
00056 G4int     G4QNuENuclearCrossSection::lastL=0;   // Last used in cross section TheLastBin
00057 G4double  G4QNuENuclearCrossSection::lastE=0.;  // Last used in cross section TheEnergy
00058 G4double* G4QNuENuclearCrossSection::lastEN=0;  // Pointer to the Energy Scale of TX & QE
00059 G4double* G4QNuENuclearCrossSection::lastTX=0;  // Pointer to the LastArray of TX function
00060 G4double* G4QNuENuclearCrossSection::lastQE=0;  // Pointer to the LastArray of QE function
00061 G4int     G4QNuENuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
00062 G4int     G4QNuENuclearCrossSection::lastN=0;   // The last N of calculated nucleus
00063 G4int     G4QNuENuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
00064 G4double  G4QNuENuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
00065 G4double  G4QNuENuclearCrossSection::lastTH=0.; // Last threshold momentum
00066 G4double  G4QNuENuclearCrossSection::lastCS=0.; // Last value of the Cross Section
00067 G4int     G4QNuENuclearCrossSection::lastI=0;   // The last position in the DAMDB
00068 std::vector<G4double*>* G4QNuENuclearCrossSection::TX = new std::vector<G4double*>;
00069 std::vector<G4double*>* G4QNuENuclearCrossSection::QE = new std::vector<G4double*>;
00070 
00071 // Returns Pointer to the G4VQCrossSection class
00072 G4VQCrossSection* G4QNuENuclearCrossSection::GetPointer()
00073 {
00074   static G4QNuENuclearCrossSection theCrossSection; //**Static body of the Cross Section**
00075   return &theCrossSection;
00076 }
00077 
00078 G4QNuENuclearCrossSection::~G4QNuENuclearCrossSection()
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 G4QNuENuclearCrossSection::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<<"G4QNENCS::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!=12)
00109   {
00110 #ifdef debug
00111     G4cout<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"---G4QNENCrossSec::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<<"G4QNENCS::GetCrSec: CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<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<<"G4QNENCrossSection::GetCrossSect:NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
00192 #endif
00193         if(pEn>lastTH)
00194         {
00195 #ifdef pdebug
00196           G4cout<<"G4QNENCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00197 #endif
00198           lastTH=pEn;
00199         }
00200       }
00201 #ifdef pdebug
00202       G4cout<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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 G4QNuENuclearCrossSection::ThresholdEnergy(G4int Z, G4int N, G4int)
00260 {
00261   //static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/GeV;
00262   //static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/GeV;
00263   //static const G4double mDeut = G4NucleiProperties::GetNuclearMass(2,1)/GeV/2.;
00264   static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
00265   static const G4double dmN=mN+mN;    // Doubled nucleon mass (2*AtomicMassUnit, GeV)
00266   static const G4double me=.00051099892; // electron mass in GeV
00267   static const G4double me2=me*me;    // Squared mass of an electron in GeV^2
00268   static const G4double thresh=me+me2/dmN; // Universal threshold in GeV
00269   // ---------
00270   //static const G4double infEn = 9.e27;
00271   G4double dN=0.;
00272   if(Z>0||N>0) dN=thresh*GeV; // @@ if upgraded, change it in a total cross section
00273   //@@ "dN=me+me2/G4NucleiProperties::GetNuclearMass(<G4double>(Z+N),<G4double>(Z)/GeV"
00274   return dN;
00275 }
00276 
00277 // The main member function giving the gamma-A cross section (E_kin in MeV, CS in mb)
00278 G4double G4QNuENuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I,
00279                                         G4int, G4int targZ, G4int targN, G4double Momentum)
00280 {
00281   static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2
00282   static const G4int nE=65;        // !! If change this, change it in GetFunctions()33/65!!
00283   static const G4int mL=nE-1;
00284   static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
00285   static const G4double dmN=mN+mN;   // Doubled nucleon mass (2*AtomicMassUnit, GeV)
00286   static const G4double me=.00051099892;// electron mass in GeV
00287   static const G4double me2=me*me;   // Squared mass of an electron in GeV^2
00288   static const G4double EMi=me+me2/dmN; // Universal threshold of the reaction in GeV
00289   static const G4double EMa=300.;    // Maximum tabulated Energy of nu_e in GeV 
00290   // *** Begin of the Associative memory for acceleration of the cross section calculations
00291   static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional)
00292   static G4bool first=true;          // Flag of initialization of the energy axis
00293   // *** End of Static Definitions (Associative Memory) ***
00294   //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Electron
00295   //G4double TotEnergy2=Momentum;
00296   onlyCS=CS;                         // Flag to calculate only CS (not TX & QE)
00297   lastE=Momentum/GeV;                // Kinetic energy of the electron neutrino (in GeV)
00298   if (lastE<=EMi)                    // Energy is below the minimum energy in the table
00299   {
00300     lastE=0.;
00301     lastSig=0.;
00302     return 0.;
00303   }
00304   G4int Z=targZ;                     // New Z, which can change the sign
00305   if(F<=0)                           // This isotope was not the last used isotop
00306   {
00307     if(F<0)                          // This isotope was found in DAMDB =-------=> RETRIEVE
00308     {
00309       lastTX =(*TX)[I];              // Pointer to the prepared TX function (same isotope)
00310       lastQE =(*QE)[I];              // Pointer to the prepared QE function (same isotope)
00311    }
00312    else                              // This isotope wasn't calculated previously => CREATE
00313    {
00314       if(first)
00315       {
00316         lastEN = new G4double[nE];   // This must be done only once!
00317         Z=-Z;                        // To explain GetFunctions that E-axis must be filled
00318         first=false;                 // To make it only once
00319       }
00320       lastTX = new G4double[nE];     // Allocate memory for the new TX function
00321       lastQE = new G4double[nE];     // Allocate memory for the new QE function
00322       G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
00323       if(res<0) G4cerr<<"*W*G4NuENuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
00324       // *** The synchronization check ***
00325       G4int sync=TX->size();
00326       if(sync!=I) G4cerr<<"***G4NuENuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
00327       TX->push_back(lastTX);
00328       QE->push_back(lastQE);
00329     } // End of creation of the new set of parameters
00330   } // End of parameters udate
00331   // =----------------= NOW Calculate the Cross Section =-----------------=
00332   if (lastE<=EMi)                   // Check that the neutrinoEnergy is higher than ThreshE
00333   {
00334     lastE=0.;
00335     lastSig=0.;
00336     return 0.;
00337   }
00338   if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
00339   {
00340     G4int chk=1;
00341     G4int ran=mL/2;
00342     G4int sep=ran;  // as a result = an index of the left edge of the interval
00343     while(ran>=2)
00344     {
00345       G4int newran=ran/2;
00346       if(lastE<=lastEN[sep]) sep-=newran;
00347       else                   sep+=newran;
00348       ran=newran;
00349       chk=chk+chk; 
00350     }
00351     if(chk+chk!=mL) G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
00352     G4double lowE=lastEN[sep];
00353     G4double highE=lastEN[sep+1];
00354     G4double lowTX=lastTX[sep];
00355     if(lastE<lowE||sep>=mL||lastE>highE)
00356       G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
00357             <<", sep="<<sep<<", mL="<<mL<<G4endl;
00358     lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
00359     if(!onlyCS)                       // Skip the differential cross-section parameters
00360     {
00361       G4double lowQE=lastQE[sep];
00362       lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
00363 #ifdef pdebug
00364       G4cout<<"G4NuENuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
00365 #endif
00366     }
00367   }
00368   else
00369   {
00370     lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
00371     lastQEL=lastQE[mL];
00372   }
00373   if(lastQEL<0.) lastQEL = 0.;
00374   if(lastSig<0.) lastSig = 0.;
00375   // The cross-sections are expected to be in mb
00376   lastSig*=mb38;
00377   if(!onlyCS) lastQEL*=mb38;
00378   return lastSig;
00379 }
00380 
00381 // Calculate the cros-section functions
00382 // ****************************************************************************************
00383 // *** This tables are the same for all lepto-nuclear reactions, only mass is different ***
00384 // ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ ***
00385 // ****************************************************************************************
00386 G4int G4QNuENuclearCrossSection::GetFunctions(G4int z, G4int n,
00387                                                      G4double* t, G4double* q, G4double* e)
00388 {
00389   static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
00390   static const G4double dmN=mN+mN;    // Doubled nucleon mass (2*AtomicMassUnit, GeV)
00391   static const G4double me=.00051099892; // electron mass in GeV
00392   static const G4double me2=me*me;    // Squared mass of an electron in GeV^2
00393   static const G4double thresh=me+me2/dmN; // Universal threshold in GeV
00394   static const G4int nE=65; // !! If change this, change it in CalculateCrossSection() !!
00395   static const G4double nuEn[nE]={thresh,
00396     .00051331,.00053602,.00056078,.00058783,.00061743,.00064990,.00068559,.00072492,
00397     .00076834,.00081641,.00086975,.00092912,.00099536,.00106950,.00115273,.00124646,
00398     .00135235,.00147241,.00160901,.00176503,.00194392,.00214986,.00238797,.00266448,
00399     .00298709,.00336531,.00381094,.00433879,.00496745,.00572047,.00662785,.00772806,
00400     .00907075,.01072050,.01276190,.01530660,.01850330,.02255110,.02771990,.03437780,
00401     .04303240,.05438970,.06944210,.08959920,.11688400,.15423600,.20597200,.27851200,
00402     .38153100,.52979600,.74616300,1.0665200,1.5480900,2.2834800,3.4251100,5.2281000,
00403     8.1270200,12.875900,20.808500,34.331200,57.877800,99.796200,176.16300,318.68200};
00404   static const G4double TOTX[nE]={0.,
00405     .00047551,.00162896,.00232785,.00292938,.00349456,.00404939,.00460908,.00518455,
00406     .00578488,.00641848,.00709376,.00781964,.00860585,.00946334,.01040460,.01144420,
00407     .01259910,.01388920,.01533830,.01697480,.01883260,.02095280,.02338500,.02618960,
00408     .02944060,.03322870,.03766580,.04289050,.04907540,.06123530,.07521120,.09034730,
00409     .10803800,.12856800,.15277600,.18174900,.21764300,.26083100,.31424900,.37935600,
00410     .45871900,.55375100,.66506600,.79039400,.92276600,1.0489000,1.1500300,1.2071700,
00411     1.2096800,1.1612200,1.0782900,.98251100,.89137400,.81472700,.75500100,.71061200,
00412     .67873900,.65619000,.64098400,.63085100,.62389900,.61664900,.61261000,.60635700};
00413   static const G4double QELX[nE]={0.,
00414   2.44084e-7,8.73147e-7,1.30540e-6,1.72196e-6,2.15765e-6,2.63171e-6,3.15996e-6,3.75836e-6,
00415   4.44474e-6,5.24008e-6,6.16982e-6,7.26537e-6,8.56595e-6,1.01211e-5,1.19937e-5,1.42647e-5,
00416   1.70384e-5,2.04506e-5,2.46796e-5,2.99611e-5,3.66090e-5,4.50456e-5,5.58425e-5,6.97817e-5,
00417   8.79417e-5,.000111825,.000143542,.000186093,.000243780,.000350295,.000498489,.000698209,
00418   .000979989,.001378320,.001949710,.002781960,.004027110,.005882030,.008710940,.013041400,
00419   .019739800,.030118300,.046183600,.070818700,.107857000,.161778000,.236873000,.336212000,
00420   .455841000,.565128000,.647837000,.701208000,.729735000,.742062000,.746495000,.748182000,
00421   .749481000,.750637000,.751471000,.752237000,.752763000,.753103000,.753159000,.753315000};
00422 
00423   // --------------------------------
00424   G4int first=0;
00425   if(z<0.)
00426   {
00427     first=1;
00428     z=-z;
00429   }
00430   if(z<1 || z>92)             // neutron and plutonium are forbidden
00431   {
00432     G4cout<<"***G4QNuENuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
00433     return -1;
00434   }
00435   for(G4int k=0; k<nE; k++)
00436   {
00437     G4double a=n+z;
00438     G4double na=n+a;
00439     G4double dn=n+n;
00440     G4double da=a+a;
00441     G4double ta=da+a;
00442     if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
00443     t[k]=TOTX[k]*nuEn[k]*(na+na)/ta+QELX[k]*(dn+dn-da)/ta; // TotalCrossSection
00444     q[k]=QELX[k]*dn/a;                                     // QuasiElasticCrossSection
00445   }
00446   return first;
00447 }
00448 
00449 // Randomize Q2 from neutrino to the scattered electron when scattering is quasi-elastic
00450 G4double G4QNuENuclearCrossSection::GetQEL_ExchangeQ2()
00451 {
00452   static const G4double me=.00051099892; // electron mass in GeV
00453   static const G4double me2=me*me;     // Squared mass of an electron in GeV^2
00454   static const G4double hme2=me2/2;    // .5*m_e^2 in GeV^2
00455   static const G4double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
00456   static const double MN2=MN*MN;       // M_N^2 in GeV^2
00457   static const G4double power=-3.5;    // direct power for the magic variable
00458   static const G4double pconv=1./power;// conversion power for the magic variable
00459   static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
00460   static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
00461   static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
00462   // Reversed table
00463   static const G4double Xl[nQ2]={1.87905e-10,
00464  .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374,
00465  .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069,
00466  .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718,
00467  .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169,
00468  .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670,
00469  .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720,
00470  .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457,
00471  .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366,
00472  .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128,
00473  .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539};
00474   // Direct table
00475   static const G4double Xmax=Xl[lQ2];
00476   static const G4double Xmin=Xl[0];
00477   static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
00478   static const G4double inl[nQ2]={0,
00479  1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811,
00480  15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039,
00481  26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876,
00482  36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950,
00483  46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377,
00484  55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162,
00485  65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852,
00486  74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767,
00487  83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106,
00488  92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000};
00489   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00490   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00491   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00492   G4double ME=Enu*MN;                 // M*E
00493   G4double dME=ME+ME;                 // 2*M*E
00494   G4double dEMN=(dEnu+MN)*ME;
00495   G4double MEm=ME-hme2;
00496   G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
00497   G4double E2M=MN*Enu2-(Enu+MN)*hme2;
00498   G4double ymax=(E2M+sqE)/dEMN;
00499   G4double ymin=(E2M-sqE)/dEMN;
00500   G4double rmin=1.-ymin;
00501   G4double rhm2E=hme2/Enu2;
00502   G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
00503   G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
00504   G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
00505   G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
00506   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00507   G4double rXi=(Xmi-Xmin)/dX;
00508   G4int    iXi=static_cast<int>(rXi);
00509   if(iXi<0) iXi=0;
00510   if(iXi>bQ2) iXi=bQ2;
00511   G4double dXi=rXi-iXi;
00512   G4double bnti=inl[iXi];
00513   G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
00514   //
00515   G4double rXa=(Xma-Xmin)/dX;
00516   G4int    iXa=static_cast<int>(rXa);
00517   if(iXa<0) iXa=0;
00518   if(iXa>bQ2) iXa=bQ2;
00519   G4double dXa=rXa-iXa;
00520   G4double bnta=inl[iXa];
00521   G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
00522   // *** Find X using the reversed table ***
00523   G4double intx=inti+(inta-inti)*G4UniformRand();
00524   G4int    intc=static_cast<int>(intx);
00525   if(intc<0) intc=0;
00526   if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
00527   G4double dint=intx-intc;
00528   G4double mX=Xl[intc];
00529   G4double X=mX+dint*(Xl[intc+1]-mX);
00530   G4double Q2=std::pow(X,pconv)-1.;
00531   return Q2*GeV*GeV;
00532 }
00533 
00534 // Randomize Q2 from neutrino to the scattered electron when scattering is not quasiElastic
00535 G4double G4QNuENuclearCrossSection::GetNQE_ExchangeQ2()
00536 {
00537   static const double mpi=.13957018;    // charged pi meson mass in GeV
00538   static const G4double me=.00051099892;// electron mass in GeV
00539   static const G4double me2=me*me;      // Squared mass of an electron in GeV^2
00540   static const G4double hme2=me2/2;     // .5*m_e^2 in GeV^2
00541   static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
00542   static const double MN2=MN*MN;        // M_N^2 in GeV^2
00543   static const double dMN=MN+MN;        // 2*M_N in GeV
00544   static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
00545   static const G4int power=7;           // direct power for the magic variable
00546   static const G4double pconv=1./power; // conversion power for the magic variable
00547   static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
00548   static const G4int lX=nX-1;           // index of the last in the Xl table
00549   static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
00550   static const G4int nE=20;             // #Of point in the El table (in GeV^2)
00551   static const G4int bE=nE-1;           // index of the last in the El table
00552   static const G4int pE=bE-1;           // index of the before last in the El table
00553   // Reversed table
00554   static const G4double X0[nX]={6.14081e-05,
00555  .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244,
00556  2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419};
00557   static const G4double X1[nX]={.00125268,
00558  .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049,
00559  4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330};
00560   static const G4double X2[nX]={.015694,
00561  1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728,
00562  10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276};
00563   static const G4double X3[nX]={.0866877,
00564  4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302,
00565  21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115};
00566   static const G4double X4[nX]={.160483,
00567  5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127,
00568  29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265};
00569   static const G4double X5[nX]={.0999307,
00570  5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398,
00571  26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649};
00572   static const G4double X6[nX]={.0276367,
00573  3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519,
00574  16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148};
00575   static const G4double X7[nX]={.00472383,
00576  2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473,
00577  9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057};
00578   static const G4double X8[nX]={.000630783,
00579  1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480,
00580  5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037};
00581   static const G4double X9[nX]={7.49179e-05,
00582  .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823,
00583  3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028};
00584   static const G4double XA[nX]={8.43437e-06,
00585  .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704,
00586  2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823};
00587   static const G4double XB[nX]={9.27028e-07,
00588  .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059,
00589  1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341};
00590   static const G4double XC[nX]={1.00807e-07,
00591  .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668,
00592  1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014};
00593   static const G4double XD[nX]={1.09102e-08,
00594  .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878,
00595  1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617};
00596   static const G4double XE[nX]={1.17831e-09,
00597  .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374,
00598  .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571};
00599   static const G4double XF[nX]={1.27141e-10,
00600  .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046,
00601  .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338};
00602   static const G4double XG[nX]={1.3713e-11,
00603  .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870,
00604  .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379};
00605   static const G4double XH[nX]={1.47877e-12,
00606  .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421,
00607  .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197};
00608   static const G4double XI[nX]={1.59454e-13,
00609  .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806,
00610  .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872};
00611   static const G4double XJ[nX]={1.71931e-14,
00612  .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147,
00613  .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813};
00614   // Direct table
00615   static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
00616                         X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
00617   static const G4double dX[nE]={
00618     (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
00619     (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
00620     (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
00621     (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
00622     (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
00623   static const G4double* Xl[nE]=
00624                              {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
00625   static const G4double I0[nX]={0,
00626  .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555,
00627  13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000};
00628   static const G4double I1[nX]={0,
00629  .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365,
00630  13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000};
00631   static const G4double I2[nX]={0,
00632  .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586,
00633  12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000};
00634   static const G4double I3[nX]={0,
00635  .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461,
00636  12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000};
00637   static const G4double I4[nX]={0,
00638  .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515,
00639  11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000};
00640   static const G4double I5[nX]={0,
00641  .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040,
00642  10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000};
00643   static const G4double I6[nX]={0,
00644  .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243,
00645  9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000};
00646   static const G4double I7[nX]={0,
00647  .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717,
00648  8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000};
00649   static const G4double I8[nX]={0,
00650  .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415,
00651  8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000};
00652   static const G4double I9[nX]={0,
00653  .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672,
00654  8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000};
00655   static const G4double IA[nX]={0,
00656  .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265,
00657  7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000};
00658   static const G4double IB[nX]={0,
00659  .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298,
00660  7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000};
00661   static const G4double IC[nX]={0,
00662  .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696,
00663  7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000};
00664   static const G4double ID[nX]={0,
00665  .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212,
00666  7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000};
00667   static const G4double IE[nX]={0,
00668  .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949,
00669  7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000};
00670   static const G4double IF[nX]={0,
00671  .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164,
00672  7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000};
00673   static const G4double IG[nX]={0,
00674  .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249,
00675  7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000};
00676   static const G4double IH[nX]={0,
00677  .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804,
00678  7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000};
00679   static const G4double II[nX]={0,
00680  .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722,
00681  7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000};
00682   static const G4double IJ[nX]={0,
00683  .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196,
00684  6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000};
00685   static const G4double* Il[nE]=
00686                              {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
00687   static const G4double lE[nE]={
00688 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
00689  2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
00690   static const G4double lEmi=lE[0];
00691   static const G4double lEma=lE[nE-1];
00692   static const G4double dlE=(lEma-lEmi)/bE;
00693   //***************************************************************************************
00694   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00695   G4double lEn=std::log(Enu);         // log(E) for interpolation
00696   G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
00697   G4int fE=static_cast<int>(rE);      // Left bin for interpolation
00698   if(fE<0) fE=0;
00699   if(fE>pE)fE=pE;
00700   G4int    sE=fE+1;                   // Right bin for interpolation
00701   G4double dE=rE-fE;                  // relative log shift from the left bin
00702   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00703   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00704   G4double Ee=Enu-me;               // Free Energy of neutrino/anti-neutrino
00705   G4double ME=Enu*MN;                 // M*E
00706   G4double dME=ME+ME;                 // 2*M*E
00707   G4double dEMN=(dEnu+MN)*ME;
00708   G4double MEm=ME-hme2;
00709   G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
00710   G4double E2M=MN*Enu2-(Enu+MN)*hme2;
00711   G4double ymax=(E2M+sqE)/dEMN;
00712   G4double ymin=(E2M-sqE)/dEMN;
00713   G4double rmin=1.-ymin;
00714   G4double rhm2E=hme2/Enu2;
00715   G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
00716   G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
00717   G4double Q2nq=Ee*dMN-mcV;
00718   if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
00719   // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
00720   G4double Rmi=Q2mi/Q2ma;
00721   G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu
00722   // --- E-interpolation must be done in a log scale ---
00723   G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
00724   G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
00725   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00726   G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
00727   G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
00728   G4double rXi=(Xmi-iXmi)/idX;
00729   G4int    iXi=static_cast<int>(rXi);
00730   if(iXi<0) iXi=0;
00731   if(iXi>bX) iXi=bX;
00732   G4double dXi=rXi-iXi;
00733   G4double bntil=Il[fE][iXi];
00734   G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
00735   G4double bntir=Il[sE][iXi];
00736   G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
00737   G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
00738   //
00739   G4double rXa=(Xma-iXmi)/idX;
00740   G4int    iXa=static_cast<int>(rXa);
00741   if(iXa<0) iXa=0;
00742   if(iXa>bX) iXa=bX;
00743   G4double dXa=rXa-iXa;
00744   G4double bntal=Il[fE][iXa];
00745   G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
00746   G4double bntar=Il[sE][iXa];
00747   G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
00748   G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
00749   //
00750   // *** Find X using the reversed table ***
00751   G4double intx=inti+(inta-inti)*G4UniformRand(); 
00752   G4int    intc=static_cast<int>(intx);
00753   if(intc<0) intc=0;
00754   if(intc>bX) intc=bX;
00755   G4double dint=intx-intc;
00756   G4double mXl=Xl[fE][intc];
00757   G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
00758   G4double mXr=Xl[sE][intc];
00759   G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
00760   G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
00761   G4double R=shift-std::pow(X,pconv);
00762   G4double Q2=R*Q2ma;
00763   return Q2*GeV*GeV;
00764 }
00765 
00766 // It returns a fraction of the direct interaction of the neutrino with quark-partons
00767 G4double G4QNuENuclearCrossSection::GetDirectPart(G4double Q2)
00768 {
00769   G4double f=Q2/4.62;
00770   G4double ff=f*f;
00771   G4double r=ff*ff;
00772   G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
00773   //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
00774   return 1.-s_value*(1.-s_value/2);
00775 }
00776 
00777 // #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
00778 G4double G4QNuENuclearCrossSection::GetNPartons(G4double Q2)
00779 {
00780   return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
00781 }
00782 
00783 // This class can provide only virtual exchange pi+ (a substitute for W+ boson)
00784 G4int G4QNuENuclearCrossSection::GetExchangePDGCode() {return 211;}

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