G4QNuMuNuclearCrossSection.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: G4QNuMuNuclearCrossSection for (nu,mu-)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: neutrino_mu -> mu 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 "G4QNuMuNuclearCrossSection.hh"
00050 #include "G4SystemOfUnits.hh"
00051 
00052 // Initialization of the
00053 G4bool    G4QNuMuNuclearCrossSection::onlyCS=true;// Flag to calculate only CS (not QE)
00054 G4double  G4QNuMuNuclearCrossSection::lastSig=0.;// Last calculated total cross section
00055 G4double  G4QNuMuNuclearCrossSection::lastQEL=0.;// Last calculated quasi-el. cross section
00056 G4int     G4QNuMuNuclearCrossSection::lastL=0;   // Last used in cross section TheLastBin
00057 G4double  G4QNuMuNuclearCrossSection::lastE=0.;  // Last used in cross section TheEnergy
00058 G4double* G4QNuMuNuclearCrossSection::lastEN=0;  // Pointer to the Energy Scale of TX & QE
00059 G4double* G4QNuMuNuclearCrossSection::lastTX=0;  // Pointer to the LastArray of TX function
00060 G4double* G4QNuMuNuclearCrossSection::lastQE=0;  // Pointer to the LastArray of QE function
00061 G4int     G4QNuMuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
00062 G4int     G4QNuMuNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
00063 G4int     G4QNuMuNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
00064 G4double  G4QNuMuNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
00065 G4double  G4QNuMuNuclearCrossSection::lastTH=0.; // Last threshold momentum
00066 G4double  G4QNuMuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
00067 G4int     G4QNuMuNuclearCrossSection::lastI=0;   // The last position in the DAMDB
00068 std::vector<G4double*>* G4QNuMuNuclearCrossSection::TX = new std::vector<G4double*>;
00069 std::vector<G4double*>* G4QNuMuNuclearCrossSection::QE = new std::vector<G4double*>;
00070 
00071 // Returns Pointer to the G4VQCrossSection class
00072 G4VQCrossSection* G4QNuMuNuclearCrossSection::GetPointer()
00073 {
00074   static G4QNuMuNuclearCrossSection theCrossSection; //**Static body of the Cross Section**
00075   return &theCrossSection;
00076 }
00077 
00078 G4QNuMuNuclearCrossSection::~G4QNuMuNuclearCrossSection()
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 G4QNuMuNuclearCrossSection::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<<"G4QNMNCS::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 pdebug
00111     G4cout<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"---G4QNMNCrossSec::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<<"G4QNMNCS::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<<"G4QNMNCrossSection::GetCrossSect:NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
00192 #endif
00193         if(pEn>lastTH)
00194         {
00195 #ifdef pdebug
00196           G4cout<<"G4QNMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00197 #endif
00198           lastTH=pEn;
00199         }
00200       }
00201 #ifdef pdebug
00202       G4cout<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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<<"G4QNMNCS::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 G4QNuMuNuclearCrossSection::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 mmu=.105658369; // Mass of a muon in GeV
00267   static const G4double mmu2=mmu*mmu;   // Squared mass of a muon in GeV^2
00268   static const G4double thresh=mmu+mmu2/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=mmu+mmu2/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 G4QNuMuNuclearCrossSection::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() (*.hh) !!
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 mmu=.105658369; // Mass of a muon in GeV
00287   static const G4double mmu2=mmu*mmu;// Squared mass of a muon in GeV^2
00288   static const G4double EMi=mmu+mmu2/dmN; // Universal threshold of the reaction in GeV
00289   static const G4double EMa=300.;    // Maximum tabulated Energy of nu_mu 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 Muon
00295   //G4double TotEnergy2=Momentum;
00296   onlyCS=CS;                         // Flag to calculate only CS (not TX & QE)
00297   lastE=Momentum/GeV;                // Kinetic energy of the muon 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*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
00324       // *** The synchronization check ***
00325       G4int sync=TX->size();
00326       if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::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*G4NuMuNuclearCS::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*G4NuMuNuclearCS::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<<"G4NuMuNuclearCS::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 G4QNuMuNuclearCrossSection::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 mmu=.105658369; // Mass of a muon in GeV
00392   static const G4double mmu2=mmu*mmu;   // Squared mass of a muon in GeV^2
00393   static const G4double thresh=mmu+mmu2/dmN; // Universal threshold in GeV
00394   static const G4int nE=65; // !! If change this, change it in GetCrossSection() (*.cc) !!
00395   static const G4double nuEn[nE]={thresh,
00396     .112039,.116079,.120416,.125076,.130090,.135494,.141324,.147626,.154445,.161838,
00397     .169864,.178594,.188105,.198485,.209836,.222272,.235923,.250941,.267497,.285789,
00398     .306045,.328530,.353552,.381466,.412689,.447710,.487101,.531538,.581820,.638893,
00399     .703886,.778147,.863293,.961275,1.07445,1.20567,1.35843,1.53701,1.74667,1.99390,
00400     2.28679,2.63542,3.05245,3.55386,4.15990,4.89644,5.79665,6.90336,8.27224,9.97606,
00401     12.1106,14.8029,18.2223,22.5968,28.2351,35.5587,45.1481,57.8086,74.6682,97.3201,
00402     128.036,170.085,228.220,309.420};
00403   static const G4double TOTX[nE]={0.,
00404     .108618,.352160,.476083,.566575,.639014,.699871,.752634,.799407,.841524,.879844,
00405     .914908,.947050,.976456,1.00321,1.02734,1.04881,1.06755,1.08349,1.09653,1.10657,
00406     1.11355,1.11739,1.11806,1.11556,1.10992,1.10124,1.08964,1.07532,1.05851,1.03950,
00407     1.01859,.996169,.972593,.948454,.923773,.899081,.874713,.850965,.828082,.806265,
00408     .785659,.766367,.748450,.731936,.716824,.703098,.690723,.679652,.669829,.661187,
00409     .653306,.646682,.640986,.636125,.631993,.628479,.625458,.622800,.620364,.616231,
00410     .614986,.612563,.609807,.606511};
00411   static const G4double QELX[nE]={0.,
00412     .012170,.040879,.057328,.070865,.083129,.094828,.106366,.118013,.129970,.142392,
00413     .155410,.169138,.183676,.199123,.215573,.233120,.251860,.271891,.293317,.316246,
00414     .340796,.367096,.395292,.425547,.458036,.491832,.524989,.556457,.585692,.612377,
00415     .636544,.657790,.676260,.692007,.705323,.716105,.724694,.731347,.736340,.740172,
00416     .742783,.744584,.745804,.746829,.747479,.747995,.748436,.749047,.749497,.749925,
00417     .750486,.750902,.751268,.751566,.752026,.752266,.752428,.752761,.752873,.753094,
00418     .753161,.753164,.753340,.753321};
00419   // --------------------------------
00420   G4int first=0;
00421   if(z<0.)
00422   {
00423     first=1;
00424     z=-z;
00425   }
00426   if(z<1 || z>92)             // neutron and plutonium are forbidden
00427   {
00428     G4cout<<"***G4QNuMuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
00429     return -1;
00430   }
00431   for(G4int k=0; k<nE; k++)
00432   {
00433     G4double a=n+z;
00434     G4double na=n+a;
00435     G4double dn=n+n;
00436     G4double da=a+a;
00437     G4double ta=da+a;
00438     if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
00439     t[k]=TOTX[k]*nuEn[k]*(na+na)/ta+QELX[k]*(dn+dn-da)/ta; // TotalCrossSection
00440     q[k]=QELX[k]*dn/a;                                     // QuasiElasticCrossSection
00441   }
00442   return first;
00443 }
00444 
00445 // Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
00446 G4double G4QNuMuNuclearCrossSection::GetQEL_ExchangeQ2()
00447 {
00448   static const G4double mmu=.105658369;// Mass of muon in GeV
00449   static const G4double mmu2=mmu*mmu;  // Squared Mass of muon in GeV^2
00450   static const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
00451   static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
00452   static const double MN2=MN*MN;       // M_N^2 in GeV^2
00453   static const G4double power=-3.5;    // direct power for the magic variable
00454   static const G4double pconv=1./power;// conversion power for the magic variable
00455   static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
00456   static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
00457   static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
00458   // Reversed table
00459   static const G4double Xl[nQ2]={1.87905e-10,
00460  .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374,
00461  .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069,
00462  .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718,
00463  .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169,
00464  .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670,
00465  .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720,
00466  .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457,
00467  .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366,
00468  .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128,
00469  .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539};
00470   // Direct table
00471   static const G4double Xmax=Xl[lQ2];
00472   static const G4double Xmin=Xl[0];
00473   static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
00474   static const G4double inl[nQ2]={0,
00475  1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811,
00476  15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039,
00477  26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876,
00478  36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950,
00479  46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377,
00480  55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162,
00481  65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852,
00482  74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767,
00483  83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106,
00484  92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000};
00485   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00486   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00487   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00488   G4double ME=Enu*MN;                 // M*E
00489   G4double dME=ME+ME;                 // 2*M*E
00490   G4double dEMN=(dEnu+MN)*ME;
00491   G4double MEm=ME-hmmu2;
00492   G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
00493   G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
00494   G4double ymax=(E2M+sqE)/dEMN;
00495   G4double ymin=(E2M-sqE)/dEMN;
00496   G4double rmin=1.-ymin;
00497   G4double rhm2E=hmmu2/Enu2;
00498   G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
00499   G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
00500   G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
00501   G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
00502   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00503   G4double rXi=(Xmi-Xmin)/dX;
00504   G4int    iXi=static_cast<int>(rXi);
00505   if(iXi<0) iXi=0;
00506   if(iXi>bQ2) iXi=bQ2;
00507   G4double dXi=rXi-iXi;
00508   G4double bnti=inl[iXi];
00509   G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
00510   //
00511   G4double rXa=(Xma-Xmin)/dX;
00512   G4int    iXa=static_cast<int>(rXa);
00513   if(iXa<0) iXa=0;
00514   if(iXa>bQ2) iXa=bQ2;
00515   G4double dXa=rXa-iXa;
00516   G4double bnta=inl[iXa];
00517   G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
00518   // *** Find X using the reversed table ***
00519   G4double intx=inti+(inta-inti)*G4UniformRand();
00520   G4int    intc=static_cast<int>(intx);
00521   if(intc<0) intc=0;
00522   if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
00523   G4double dint=intx-intc;
00524   G4double mX=Xl[intc];
00525   G4double X=mX+dint*(Xl[intc+1]-mX);
00526   G4double Q2=std::pow(X,pconv)-1.;
00527   return Q2*GeV*GeV;
00528 }
00529 
00530 // Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
00531 G4double G4QNuMuNuclearCrossSection::GetNQE_ExchangeQ2()
00532 {
00533   static const double mpi=.13957018;    // charged pi meson mass in GeV
00534   static const G4double mmu=.105658369; // Mass of muon in GeV
00535   static const G4double mmu2=mmu*mmu;   // Squared Mass of muon in GeV^2
00536   static const double hmmu2=mmu2/2;     // .5*m_mu^2 in GeV^2
00537   static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
00538   static const double MN2=MN*MN;        // M_N^2 in GeV^2
00539   static const double dMN=MN+MN;        // 2*M_N in GeV
00540   static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
00541   static const G4int power=7;           // direct power for the magic variable
00542   static const G4double pconv=1./power; // conversion power for the magic variable
00543   static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
00544   static const G4int lX=nX-1;           // index of the last in the Xl table
00545   static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
00546   static const G4int nE=20;             // #Of point in the El table (in GeV^2)
00547   static const G4int bE=nE-1;           // index of the last in the El table
00548   static const G4int pE=bE-1;           // index of the before last in the El table
00549   // Reversed table
00550   static const G4double X0[nX]={6.14081e-05,
00551  .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244,
00552  2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419};
00553   static const G4double X1[nX]={.00125268,
00554  .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049,
00555  4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330};
00556   static const G4double X2[nX]={.015694,
00557  1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728,
00558  10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276};
00559   static const G4double X3[nX]={.0866877,
00560  4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302,
00561  21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115};
00562   static const G4double X4[nX]={.160483,
00563  5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127,
00564  29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265};
00565   static const G4double X5[nX]={.0999307,
00566  5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398,
00567  26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649};
00568   static const G4double X6[nX]={.0276367,
00569  3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519,
00570  16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148};
00571   static const G4double X7[nX]={.00472383,
00572  2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473,
00573  9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057};
00574   static const G4double X8[nX]={.000630783,
00575  1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480,
00576  5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037};
00577   static const G4double X9[nX]={7.49179e-05,
00578  .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823,
00579  3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028};
00580   static const G4double XA[nX]={8.43437e-06,
00581  .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704,
00582  2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823};
00583   static const G4double XB[nX]={9.27028e-07,
00584  .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059,
00585  1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341};
00586   static const G4double XC[nX]={1.00807e-07,
00587  .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668,
00588  1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014};
00589   static const G4double XD[nX]={1.09102e-08,
00590  .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878,
00591  1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617};
00592   static const G4double XE[nX]={1.17831e-09,
00593  .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374,
00594  .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571};
00595   static const G4double XF[nX]={1.27141e-10,
00596  .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046,
00597  .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338};
00598   static const G4double XG[nX]={1.3713e-11,
00599  .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870,
00600  .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379};
00601   static const G4double XH[nX]={1.47877e-12,
00602  .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421,
00603  .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197};
00604   static const G4double XI[nX]={1.59454e-13,
00605  .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806,
00606  .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872};
00607   static const G4double XJ[nX]={1.71931e-14,
00608  .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147,
00609  .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813};
00610   // Direct table
00611   static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
00612                         X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
00613   static const G4double dX[nE]={
00614     (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
00615     (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
00616     (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
00617     (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
00618     (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
00619   static const G4double* Xl[nE]=
00620                              {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
00621   static const G4double I0[nX]={0,
00622  .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555,
00623  13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000};
00624   static const G4double I1[nX]={0,
00625  .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365,
00626  13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000};
00627   static const G4double I2[nX]={0,
00628  .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586,
00629  12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000};
00630   static const G4double I3[nX]={0,
00631  .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461,
00632  12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000};
00633   static const G4double I4[nX]={0,
00634  .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515,
00635  11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000};
00636   static const G4double I5[nX]={0,
00637  .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040,
00638  10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000};
00639   static const G4double I6[nX]={0,
00640  .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243,
00641  9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000};
00642   static const G4double I7[nX]={0,
00643  .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717,
00644  8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000};
00645   static const G4double I8[nX]={0,
00646  .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415,
00647  8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000};
00648   static const G4double I9[nX]={0,
00649  .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672,
00650  8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000};
00651   static const G4double IA[nX]={0,
00652  .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265,
00653  7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000};
00654   static const G4double IB[nX]={0,
00655  .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298,
00656  7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000};
00657   static const G4double IC[nX]={0,
00658  .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696,
00659  7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000};
00660   static const G4double ID[nX]={0,
00661  .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212,
00662  7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000};
00663   static const G4double IE[nX]={0,
00664  .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949,
00665  7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000};
00666   static const G4double IF[nX]={0,
00667  .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164,
00668  7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000};
00669   static const G4double IG[nX]={0,
00670  .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249,
00671  7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000};
00672   static const G4double IH[nX]={0,
00673  .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804,
00674  7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000};
00675   static const G4double II[nX]={0,
00676  .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722,
00677  7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000};
00678   static const G4double IJ[nX]={0,
00679  .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196,
00680  6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000};
00681   static const G4double* Il[nE]=
00682                              {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
00683   static const G4double lE[nE]={
00684 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
00685  2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
00686   static const G4double lEmi=lE[0];
00687   static const G4double lEma=lE[nE-1];
00688   static const G4double dlE=(lEma-lEmi)/bE;
00689   //***************************************************************************************
00690   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00691   G4double lEn=std::log(Enu);         // log(E) for interpolation
00692   G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
00693   G4int fE=static_cast<int>(rE);      // Left bin for interpolation
00694   if(fE<0) fE=0;
00695   if(fE>pE)fE=pE;
00696   G4int    sE=fE+1;                   // Right bin for interpolation
00697   G4double dE=rE-fE;                  // relative log shift from the left bin
00698   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00699   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00700   G4double Emu=Enu-mmu;               // Free Energy of neutrino/anti-neutrino
00701   G4double ME=Enu*MN;                 // M*E
00702   G4double dME=ME+ME;                 // 2*M*E
00703   G4double dEMN=(dEnu+MN)*ME;
00704   G4double MEm=ME-hmmu2;
00705   G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
00706   G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
00707   G4double ymax=(E2M+sqE)/dEMN;
00708   G4double ymin=(E2M-sqE)/dEMN;
00709   G4double rmin=1.-ymin;
00710   G4double rhm2E=hmmu2/Enu2;
00711   G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
00712   G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
00713   G4double Q2nq=Emu*dMN-mcV;
00714   if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
00715   // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
00716   G4double Rmi=Q2mi/Q2ma;
00717   G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu
00718   // --- E-interpolation must be done in a log scale ---
00719   G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
00720   G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
00721   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00722   G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
00723   G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
00724   G4double rXi=(Xmi-iXmi)/idX;
00725   G4int    iXi=static_cast<int>(rXi);
00726   if(iXi<0) iXi=0;
00727   if(iXi>bX) iXi=bX;
00728   G4double dXi=rXi-iXi;
00729   G4double bntil=Il[fE][iXi];
00730   G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
00731   G4double bntir=Il[sE][iXi];
00732   G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
00733   G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
00734   //
00735   G4double rXa=(Xma-iXmi)/idX;
00736   G4int    iXa=static_cast<int>(rXa);
00737   if(iXa<0) iXa=0;
00738   if(iXa>bX) iXa=bX;
00739   G4double dXa=rXa-iXa;
00740   G4double bntal=Il[fE][iXa];
00741   G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
00742   G4double bntar=Il[sE][iXa];
00743   G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
00744   G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
00745   //
00746   // *** Find X using the reversed table ***
00747   G4double intx=inti+(inta-inti)*G4UniformRand(); 
00748   G4int    intc=static_cast<int>(intx);
00749   if(intc<0) intc=0;
00750   if(intc>bX) intc=bX;
00751   G4double dint=intx-intc;
00752   G4double mXl=Xl[fE][intc];
00753   G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
00754   G4double mXr=Xl[sE][intc];
00755   G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
00756   G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
00757   G4double R=shift-std::pow(X,pconv);
00758   G4double Q2=R*Q2ma;
00759   return Q2*GeV*GeV;
00760 }
00761 
00762 // It returns a fraction of the direct interaction of the neutrino with quark-partons
00763 G4double G4QNuMuNuclearCrossSection::GetDirectPart(G4double Q2)
00764 {
00765   G4double f=Q2/4.62;
00766   G4double ff=f*f;
00767   G4double r=ff*ff;
00768   G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
00769   //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
00770   return 1.-s_value*(1.-s_value/2);
00771 }
00772 
00773 // #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
00774 G4double G4QNuMuNuclearCrossSection::GetNPartons(G4double Q2)
00775 {
00776   return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
00777 }
00778 
00779 // This class can provide only virtual exchange pi+ (a substitute for W+ boson)
00780 G4int G4QNuMuNuclearCrossSection::GetExchangePDGCode() {return 211;}

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