G4QANuMuNuclearCrossSection.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: G4QANuMuNuclearCrossSection for (anu_mu,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 desctription: antinu_mu -> mu cross-section
00040 // ----------------------------------------------------------------------------------------
00041 
00042 //#define debug
00043 //#define edebug
00044 //#define pdebug
00045 //#define ppdebug
00046 //#define tdebug
00047 //#define sdebug
00048 
00049 #include "G4QANuMuNuclearCrossSection.hh"
00050 #include "G4SystemOfUnits.hh"
00051 
00052 // Initialization of the
00053 G4bool    G4QANuMuNuclearCrossSection::onlyCS=true;//Flag to calculate only CS (not QE)
00054 G4double  G4QANuMuNuclearCrossSection::lastSig=0.;//Last calculated total cross section
00055 G4double  G4QANuMuNuclearCrossSection::lastQEL=0.;//Last calculated quasi-el. cross section
00056 G4int     G4QANuMuNuclearCrossSection::lastL=0;   //Last used in cross section TheLastBin
00057 G4double  G4QANuMuNuclearCrossSection::lastE=0.;  //Last used in cross section TheEnergy
00058 G4double* G4QANuMuNuclearCrossSection::lastEN=0;  //Pointer to the Energy Scale of TX & QE
00059 G4double* G4QANuMuNuclearCrossSection::lastTX=0;  //Pointer to the LastArray of TX function
00060 G4double* G4QANuMuNuclearCrossSection::lastQE=0;  //Pointer to the LastArray of QE function
00061 G4int     G4QANuMuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
00062 G4int     G4QANuMuNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
00063 G4int     G4QANuMuNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
00064 G4double  G4QANuMuNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
00065 G4double  G4QANuMuNuclearCrossSection::lastTH=0.; // Last threshold momentum
00066 G4double  G4QANuMuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
00067 G4int     G4QANuMuNuclearCrossSection::lastI=0;   // The last position in the DAMDB
00068 std::vector<G4double*>* G4QANuMuNuclearCrossSection::TX = new std::vector<G4double*>;
00069 std::vector<G4double*>* G4QANuMuNuclearCrossSection::QE = new std::vector<G4double*>;
00070 
00071 // Returns Pointer to the G4VQCrossSection class
00072 G4VQCrossSection* G4QANuMuNuclearCrossSection::GetPointer()
00073 {
00074   static G4QANuMuNuclearCrossSection theCrossSection;//**Static body of the Cross Section**
00075   return &theCrossSection;
00076 }
00077 
00078 G4QANuMuNuclearCrossSection::~G4QANuMuNuclearCrossSection()
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 G4QANuMuNuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
00091                                                       G4int tgZ, G4int tgN, G4int pPDG)
00092 {
00093   static G4int j;                      // A#0f records found in DB for this projectile
00094   static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
00095   static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
00096   static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
00097   static std::vector <G4double> colP;  // Vector of last momenta for the reaction
00098   static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
00099   static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
00100   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
00101   G4double pEn=pMom;
00102 #ifdef debug
00103   G4cout<<"G4QAMNCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
00104         <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
00105         <<colN.size()<<G4endl;
00106   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00107 #endif
00108   if(pPDG!=-14)
00109   {
00110 #ifdef debug
00111     G4cout<<"G4QAMNCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00112     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00113 #endif
00114     return 0.;                         // projectile PDG=0 is a mistake (?!) @@
00115   }
00116   G4bool in=false;                     // By default the isotope must be found in the AMDB
00117   if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
00118   {
00119     in = false;                        // By default the isotope haven't be found in AMDB  
00120     lastP   = 0.;                      // New momentum history (nothing to compare with)
00121     lastPDG = pPDG;                    // The last PDG of the projectile
00122     lastN   = tgN;                     // The last N of the calculated nucleus
00123     lastZ   = tgZ;                     // The last Z of the calculated nucleus
00124     lastI   = colN.size();             // Size of the Associative Memory DB in the heap
00125     j  = 0;                            // A#0f records found in DB for this projectile
00126     if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found
00127     {                                  // The nucleus with projPDG is found in AMDB
00128       if(colN[i]==tgN && colZ[i]==tgZ)
00129       {
00130         lastI=i;
00131         lastTH =colTH[i];                // Last THreshold (A-dependent)
00132 #ifdef pdebug
00133         G4cout<<"G4QAMNCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00134         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00135 #endif
00136         if(pEn<=lastTH)
00137         {
00138 #ifdef pdebug
00139           G4cout<<"G4QAMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl;
00140           //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00141 #endif
00142           return 0.;                     // Energy is below the Threshold value
00143         }
00144         lastP  =colP [i];                // Last Momentum  (A-dependent)
00145         lastCS =colCS[i];                // Last CrossSect (A-dependent)
00146         if(std::fabs(lastP/pMom-1.)<tolerance)
00147         {
00148 #ifdef pdebug
00149           G4cout<<"G4QAMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00150           //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00151 #endif
00152           return lastCS*millibarn;     // Use theLastCS
00153         }
00154         in = true;                       // This is the case when the isotop is found in DB
00155         // Momentum pMom is in IU ! @@ Units
00156 #ifdef pdebug
00157         G4cout<<"G4QAMNCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00158 #endif
00159         lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
00160 #ifdef pdebug
00161         G4cout<<"G4QAMNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00162         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00163 #endif
00164         if(lastCS<=0. && pEn>lastTH)    // Correct the threshold
00165         {
00166 #ifdef pdebug
00167           G4cout<<"G4QAMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00168 #endif
00169           lastTH=pEn;
00170         }
00171         break;                           // Go out of the LOOP
00172       }
00173 #ifdef pdebug
00174       G4cout<<"---G4QAMNCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
00175             <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl;
00176       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00177 #endif
00178       j++;                             // Increment a#0f records found in DB for this pPDG
00179     }
00180     if(!in)                            // This nucleus has not been calculated previously
00181     {
00182 #ifdef pdebug
00183       G4cout<<"G4QAMNCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<lastI<<G4endl;
00184 #endif
00186       lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
00187       if(lastCS<=0.)
00188       {
00189         lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
00190 #ifdef pdebug
00191         G4cout<<"G4QAMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
00192 #endif
00193         if(pEn>lastTH)
00194         {
00195 #ifdef pdebug
00196           G4cout<<"G4QAMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00197 #endif
00198           lastTH=pEn;
00199         }
00200       }
00201 #ifdef pdebug
00202       G4cout<<"G4QAMNCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00203       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00204 #endif
00205       colN.push_back(tgN);
00206       colZ.push_back(tgZ);
00207       colPDG.push_back(pPDG);
00208       colP.push_back(pMom);
00209       colTH.push_back(lastTH);
00210       colCS.push_back(lastCS);
00211 #ifdef pdebug
00212       G4cout<<"G4QAMNCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl;
00213       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00214 #endif
00215       return lastCS*millibarn;
00216     } // End of creation of the new set of parameters
00217     else
00218     {
00219 #ifdef pdebug
00220       G4cout<<"G4QAMNCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00221 #endif
00222       colP[lastI]=pMom;
00223       colPDG[lastI]=pPDG;
00224       colCS[lastI]=lastCS;
00225     }
00226   } // End of parameters udate
00227   else if(pEn<=lastTH)
00228   {
00229 #ifdef pdebug
00230     G4cout<<"G4QAMNCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00231     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00232 #endif
00233     return 0.;                         // Momentum is below the Threshold Value -> CS=0
00234   }
00235   else if(std::fabs(lastP/pMom-1.)<tolerance)
00236   {
00237 #ifdef pdebug
00238     G4cout<<"G4QAMNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00239     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00240 #endif
00241     return lastCS*millibarn;     // Use theLastCS
00242   }
00243   else
00244   {
00245 #ifdef pdebug
00246     G4cout<<"G4QAMNCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00247 #endif
00248     lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00249     lastP=pMom;
00250   }
00251 #ifdef pdebug
00252   G4cout<<"G4QAMNCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00253   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00254 #endif
00255   return lastCS*millibarn;
00256 }
00257 
00258 // Gives the threshold energy = the same for all nuclei (@@ can be reduced for hevy nuclei)
00259 G4double G4QANuMuNuclearCrossSection::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 G4QANuMuNuclearCrossSection::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 antiNuEnergy 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 G4QANuMuNuclearCrossSection::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     .077498,.247583,.329691,.386384,.429087,.462699,.489899,.512316,.530996,.546614,
00405     .559616,.570292,.578840,.585395,.590053,.593083,.594197,.593614,.591396,.587611,
00406     .582335,.575653,.567667,.558490,.548417,.537270,.525352,.512825,.499857,.486620,
00407     .473283,.460014,.446970,.434294,.422116,.410656,.399782,.389665,.380349,.371860,
00408     .364207,.357387,.351388,.346192,.341778,.338122,.335198,.332980,.331439,.330544,
00409     .330263,.330558,.331391,.332718,.334494,.336667,.339182,.341697,.344470,.348125,
00410     .351322,.354481,.357507,.359239};
00411   static const G4double QELX[nE]={0.,
00412     .008683,.028739,.039700,.048327,.055820,.062693,.069235,.075631,.082010,.088463,
00413     .095059,.101851,.108883,.116192,.123814,.131826,.140185,.148962,.158197,.167933,
00414     .178221,.189119,.200700,.213045,.226326,.240454,.255277,.270612,.286388,.302608,
00415     .319318,.336582,.354468,.373031,.392427,.412445,.433146,.454448,.476222,.498289,
00416     .520430,.542558,.564130,.585003,.604928,.623680,.641266,.657255,.671704,.684586,
00417     .696111,.706028,.714553,.721951,.728085,.733182,.737348,.740958,.743716,.746059,
00418     .747806,.749129,.750331,.751100};
00419 
00420   // --------------------------------
00421   G4int first=0;
00422   if(z<0.)
00423   {
00424     first=1;
00425     z=-z;
00426   }
00427   if(z<1 || z>92)             // neutron & plutonium are forbidden
00428   {
00429     G4cout<<"**G4QANuMuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
00430     return -1;
00431   }
00432   for(G4int k=0; k<nE; k++)
00433   {
00434     G4double a=n+z;
00435     G4double za=z+a;
00436     G4double dz=z+z;
00437     G4double da=a+a;
00438     G4double ta=da+a;
00439     if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
00440     t[k]=TOTX[k]*nuEn[k]*(za+za)/ta+QELX[k]*(dz+dz-da)/ta; // TotalCrossSection
00441     q[k]=QELX[k]*dz/a;                                     // QuasiElasticCrossSection
00442   }
00443   return first;
00444 }
00445 
00446 // Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
00447 G4double G4QANuMuNuclearCrossSection::GetQEL_ExchangeQ2()
00448 {
00449   static const G4double mmu=.105658369;// Mass of muon in GeV
00450   static const G4double mmu2=mmu*mmu;  // Squared Mass of muon in GeV^2
00451   static const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
00452   static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
00453   static const double MN2=MN*MN;       // M_N^2 in GeV^2
00454   static const G4double power=-3.5;    // direct power for the magic variable
00455   static const G4double pconv=1./power;// conversion power for the magic variable
00456   static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
00457   static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
00458   static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
00459   // Reversed table
00460   static const G4double Xl[nQ2]={5.20224e-16,
00461  .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898,
00462  .093878, .103332, .112844, .122410, .132023, .141680, .151376, .161109, .170875, .180672,
00463  .190499, .200352, .210230, .220131, .230055, .239999, .249963, .259945, .269944, .279960,
00464  .289992, .300039, .310099, .320173, .330260, .340359, .350470, .360592, .370724, .380867,
00465  .391019, .401181, .411352, .421531, .431719, .441915, .452118, .462329, .472547, .482771,
00466  .493003, .503240, .513484, .523734, .533989, .544250, .554517, .564788, .575065, .585346,
00467  .595632, .605923, .616218, .626517, .636820, .647127, .657438, .667753, .678072, .688394,
00468  .698719, .709048, .719380, .729715, .740053, .750394, .760738, .771085, .781434, .791786,
00469  .802140, .812497, .822857, .833219, .843582, .853949, .864317, .874687, .885060, .895434,
00470  .905810, .916188, .926568, .936950, .947333, .957719, .968105, .978493, .988883, .999275};
00471    // Direct table
00472   static const G4double Xmax=Xl[lQ2];
00473   static const G4double Xmin=Xl[0];
00474   static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
00475   static const G4double inl[nQ2]={0,
00476  1.52225, 2.77846, 3.96651, 5.11612, 6.23990, 7.34467, 8.43466, 9.51272, 10.5809, 11.6406,
00477  12.6932, 13.7394, 14.7801, 15.8158, 16.8471, 17.8743, 18.8979, 19.9181, 20.9353, 21.9496,
00478  22.9614, 23.9707, 24.9777, 25.9826, 26.9855, 27.9866, 28.9860, 29.9837, 30.9798, 31.9745,
00479  32.9678, 33.9598, 34.9505, 35.9400, 36.9284, 37.9158, 38.9021, 39.8874, 40.8718, 41.8553,
00480  42.8379, 43.8197, 44.8007, 45.7810, 46.7605, 47.7393, 48.7174, 49.6950, 50.6718, 51.6481,
00481  52.6238, 53.5990, 54.5736, 55.5476, 56.5212, 57.4943, 58.4670, 59.4391, 60.4109, 61.3822,
00482  62.3531, 63.3236, 64.2937, 65.2635, 66.2329, 67.2019, 68.1707, 69.1390, 70.1071, 71.0748,
00483  72.0423, 73.0095, 73.9763, 74.9429, 75.9093, 76.8754, 77.8412, 78.8068, 79.7721, 80.7373,
00484  81.7022, 82.6668, 83.6313, 84.5956, 85.5596, 86.5235, 87.4872, 88.4507, 89.4140, 90.3771,
00485  91.3401, 92.3029, 93.2656, 94.2281, 95.1904, 96.1526, 97.1147, 98.0766, 99.0384, 100.000};
00486   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00487   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00488   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00489   G4double ME=Enu*MN;                 // M*E
00490   G4double dME=ME+ME;                 // 2*M*E
00491   G4double dEMN=(dEnu+MN)*ME;
00492   G4double MEm=ME-hmmu2;
00493   G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
00494   G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
00495   G4double ymax=(E2M+sqE)/dEMN;
00496   G4double ymin=(E2M-sqE)/dEMN;
00497   G4double rmin=1.-ymin;
00498   G4double rhm2E=hmmu2/Enu2;
00499   G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
00500   G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
00501   G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
00502   G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
00503   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00504   G4double rXi=(Xmi-Xmin)/dX;
00505   G4int    iXi=static_cast<int>(rXi);
00506   if(iXi<0) iXi=0;
00507   if(iXi>bQ2) iXi=bQ2;
00508   G4double dXi=rXi-iXi;
00509   G4double bnti=inl[iXi];
00510   G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
00511   //
00512   G4double rXa=(Xma-Xmin)/dX;
00513   G4int    iXa=static_cast<int>(rXa);
00514   if(iXa<0) iXa=0;
00515   if(iXa>bQ2) iXa=bQ2;
00516   G4double dXa=rXa-iXa;
00517   G4double bnta=inl[iXa];
00518   G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
00519   // *** Find X using the reversed table ***
00520   G4double intx=inti+(inta-inti)*G4UniformRand();
00521   G4int    intc=static_cast<int>(intx);
00522   if(intc<0) intc=0;
00523   if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
00524   G4double dint=intx-intc;
00525   G4double mX=Xl[intc];
00526   G4double X=mX+dint*(Xl[intc+1]-mX);
00527   G4double Q2=std::pow(X,pconv)-1.;
00528   return Q2*GeV*GeV;
00529 }
00530 
00531 // Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
00532 G4double G4QANuMuNuclearCrossSection::GetNQE_ExchangeQ2()
00533 {
00534   static const double mpi=.13957018;    // charged pi meson mass in GeV
00535   static const G4double mmu=.105658369; // Mass of muon in GeV
00536   static const G4double mmu2=mmu*mmu;   // Squared Mass of muon in GeV^2
00537   static const double hmmu2=mmu2/2;     // .5*m_mu^2 in GeV^2
00538   static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
00539   static const double MN2=MN*MN;        // M_N^2 in GeV^2
00540   static const double dMN=MN+MN;        // 2*M_N in GeV
00541   static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
00542   static const G4int power=7;           // direct power for the magic variable
00543   static const G4double pconv=1./power; // conversion power for the magic variable
00544   static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
00545   static const G4int lX=nX-1;           // index of the last in the Xl table
00546   static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
00547   static const G4int nE=20;             // #Of point in the El table (in GeV^2)
00548   static const G4int bE=nE-1;           // index of the last in the El table
00549   static const G4int pE=bE-1;           // index of the before last in the El table
00550   // Reversed table
00551   static const G4double X0[nX]={5.21412e-05,
00552  .437860, .681908, .891529, 1.08434, 1.26751, 1.44494, 1.61915, 1.79198, 1.96493, 2.13937,
00553  2.31664, 2.49816, 2.68559, 2.88097, 3.08705, 3.30774, 3.54917, 3.82233, 4.15131, 4.62182};
00554   static const G4double X1[nX]={.00102591,
00555  1.00443, 1.55828, 2.03126, 2.46406, 2.87311, 3.26723, 3.65199, 4.03134, 4.40835, 4.78561,
00556  5.16549, 5.55031, 5.94252, 6.34484, 6.76049, 7.19349, 7.64917, 8.13502, 8.66246, 9.25086};
00557   static const G4double X2[nX]={.0120304,
00558  2.59903, 3.98637, 5.15131, 6.20159, 7.18024, 8.10986, 9.00426, 9.87265, 10.7217, 11.5564,
00559  12.3808, 13.1983, 14.0116, 14.8234, 15.6359, 16.4515, 17.2723, 18.1006, 18.9386, 19.7892};
00560   static const G4double X3[nX]={.060124,
00561  5.73857, 8.62595, 10.9849, 13.0644, 14.9636, 16.7340, 18.4066, 20.0019, 21.5342, 23.0142,
00562  24.4497, 25.8471, 27.2114, 28.5467, 29.8564, 31.1434, 32.4102, 33.6589, 34.8912, 36.1095};
00563   static const G4double X4[nX]={.0992363,
00564  8.23746, 12.1036, 15.1740, 17.8231, 20.1992, 22.3792, 24.4092, 26.3198, 28.1320, 29.8615,
00565  31.5200, 33.1169, 34.6594, 36.1536, 37.6044, 39.0160, 40.3920, 41.7353, 43.0485, 44.3354};
00566   static const G4double X5[nX]={.0561127,
00567  7.33661, 10.5694, 13.0778, 15.2061, 17.0893, 18.7973, 20.3717, 21.8400, 23.2211, 24.5291,
00568  25.7745, 26.9655, 28.1087, 29.2094, 30.2721, 31.3003, 32.2972, 33.2656, 34.2076, 35.1265};
00569   static const G4double X6[nX]={.0145859,
00570  4.81774, 6.83565, 8.37399, 9.66291, 10.7920, 11.8075, 12.7366, 13.5975, 14.4025, 15.1608,
00571  15.8791, 16.5628, 17.2162, 17.8427, 18.4451, 19.0259, 19.5869, 20.1300, 20.6566, 21.1706};
00572   static const G4double X7[nX]={.00241155,
00573  2.87095, 4.02492, 4.89243, 5.61207, 6.23747, 6.79613, 7.30433, 7.77270, 8.20858, 8.61732,
00574  9.00296, 9.36863, 9.71682, 10.0495, 10.3684, 10.6749, 10.9701, 11.2550, 11.5306, 11.7982};
00575   static const G4double X8[nX]={.000316863,
00576  1.76189, 2.44632, 2.95477, 3.37292, 3.73378, 4.05420, 4.34415, 4.61009, 4.85651, 5.08666,
00577  5.30299, 5.50738, 5.70134, 5.88609, 6.06262, 6.23178, 6.39425, 6.55065, 6.70149, 6.84742};
00578   static const G4double X9[nX]={3.73544e-05,
00579  1.17106, 1.61289, 1.93763, 2.20259, 2.42976, 2.63034, 2.81094, 2.97582, 3.12796, 3.26949,
00580  3.40202, 3.52680, 3.64482, 3.75687, 3.86360, 3.96557, 4.06323, 4.15697, 4.24713, 4.33413};
00581   static const G4double XA[nX]={4.19131e-06,
00582  .849573, 1.16208, 1.38955, 1.57379, 1.73079, 1.86867, 1.99221, 2.10451, 2.20770, 2.30332,
00583  2.39252, 2.47622, 2.55511, 2.62977, 2.70066, 2.76818, 2.83265, 2.89437, 2.95355, 3.01051};
00584   static const G4double XB[nX]={4.59981e-07,
00585  .666131, .905836, 1.07880, 1.21796, 1.33587, 1.43890, 1.53080, 1.61399, 1.69011, 1.76040,
00586  1.82573, 1.88682, 1.94421, 1.99834, 2.04959, 2.09824, 2.14457, 2.18878, 2.23107, 2.27162};
00587   static const G4double XC[nX]={4.99861e-08,
00588  .556280, .752730, .893387, 1.00587, 1.10070, 1.18317, 1.25643, 1.32247, 1.38269, 1.43809,
00589  1.48941, 1.53724, 1.58203, 1.62416, 1.66391, 1.70155, 1.73728, 1.77128, 1.80371, 1.83473};
00590   static const G4double XD[nX]={5.40832e-09,
00591  .488069, .657650, .778236, .874148, .954621, 1.02432, 1.08599, 1.14138, 1.19172, 1.23787,
00592  1.28049, 1.32008, 1.35705, 1.39172, 1.42434, 1.45514, 1.48429, 1.51197, 1.53829, 1.56339};
00593   static const G4double XE[nX]={5.84029e-10,
00594  .445057, .597434, .705099, .790298, .861468, .922865, .976982, 1.02542, 1.06930, 1.10939,
00595  1.14630, 1.18050, 1.21233, 1.24208, 1.27001, 1.29630, 1.32113, 1.34462, 1.36691, 1.38812};
00596   static const G4double XF[nX]={6.30137e-11,
00597  .418735, .560003, .659168, .737230, .802138, .857898, .906854, .950515, .989915, 1.02580,
00598  1.05873, 1.08913, 1.11734, 1.14364, 1.16824, 1.19133, 1.21306, 1.23358, 1.25298, 1.27139};
00599   static const G4double XG[nX]={6.79627e-12,
00600  .405286, .539651, .633227, .706417, .766929, .818642, .863824, .903931, .939963, .972639,
00601  1.00250, 1.02995, 1.05532, 1.07887, 1.10082, 1.12134, 1.14058, 1.15867, 1.17572, 1.19183};
00602   static const G4double XH[nX]={7.32882e-13,
00603  .404391, .535199, .625259, .695036, .752243, .800752, .842823, .879906, .912994, .942802,
00604  .969862, .994583, 1.01729, 1.03823, 1.05763, 1.07566, 1.09246, 1.10816, 1.12286, 1.13667};
00605   static const G4double XI[nX]={7.90251e-14,
00606  .418084, .548382, .636489, .703728, .758106, .803630, .842633, .876608, .906576, .933269,
00607  .957233, .978886, .998556, 1.01651, 1.03295, 1.04807, 1.06201, 1.07489, 1.08683, 1.09792};
00608   static const G4double XJ[nX]={8.52083e-15,
00609  .447299, .579635, .666780, .731788, .783268, .825512, .861013, .891356, .917626, .940597,
00610  .960842, .978802, .994820, 1.00917, 1.02208, 1.03373, 1.04427, 1.05383, 1.06253, 1.07046};
00611   // Direct table
00612   static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
00613                         X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
00614   static const G4double dX[nE]={
00615     (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
00616     (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
00617     (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
00618     (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
00619     (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
00620   static const G4double* Xl[nE]=
00621                              {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
00622   static const G4double I0[nX]={0,
00623  .354631, 1.08972, 2.05138, 3.16564, 4.38343, 5.66828, 6.99127, 8.32858, 9.65998, 10.9680,
00624  12.2371, 13.4536, 14.6050, 15.6802, 16.6686, 17.5609, 18.3482, 19.0221, 19.5752, 20.0000};
00625   static const G4double I1[nX]={0,
00626  .281625, .877354, 1.67084, 2.60566, 3.64420, 4.75838, 5.92589, 7.12829, 8.34989, 9.57708,
00627  10.7978, 12.0014, 13.1781, 14.3190, 15.4162, 16.4620, 17.4496, 18.3724, 19.2245, 20.0000};
00628   static const G4double I2[nX]={0,
00629  .201909, .642991, 1.24946, 1.98463, 2.82370, 3.74802, 4.74263, 5.79509, 6.89474, 8.03228,
00630  9.19947, 10.3889, 11.5938, 12.8082, 14.0262, 15.2427, 16.4527, 17.6518, 18.8356, 20.0000};
00631   static const G4double I3[nX]={0,
00632  .140937, .461189, .920216, 1.49706, 2.17728, 2.94985, 3.80580, 4.73758, 5.73867, 6.80331,
00633  7.92637, 9.10316, 10.3294, 11.6013, 12.9150, 14.2672, 15.6548, 17.0746, 18.5239, 20.0000};
00634   static const G4double I4[nX]={0,
00635  .099161, .337358, .694560, 1.16037, 1.72761, 2.39078, 3.14540, 3.98768, 4.91433, 5.92245,
00636  7.00942, 8.17287, 9.41060, 10.7206, 12.1010, 13.5500, 15.0659, 16.6472, 18.2924, 20.0000};
00637   static const G4double I5[nX]={0,
00638  .071131, .255084, .543312, .932025, 1.41892, 2.00243, 2.68144, 3.45512, 4.32283, 5.28411,
00639  6.33859, 7.48602, 8.72621, 10.0590, 11.4844, 13.0023, 14.6128, 16.3158, 18.1115, 20.0000};
00640   static const G4double I6[nX]={0,
00641  .053692, .202354, .443946, .778765, 1.20774, 1.73208, 2.35319, 3.07256, 3.89177, 4.81249,
00642  5.83641, 6.96528, 8.20092, 9.54516, 10.9999, 12.5670, 14.2486, 16.0466, 17.9630, 20.0000};
00643   static const G4double I7[nX]={0,
00644  .043065, .168099, .376879, .672273, 1.05738, 1.53543, 2.10973, 2.78364, 3.56065, 4.44429,
00645  5.43819, 6.54610, 7.77186, 9.11940, 10.5928, 12.1963, 13.9342, 15.8110, 17.8313, 20.0000};
00646   static const G4double I8[nX]={0,
00647  .036051, .143997, .327877, .592202, .941572, 1.38068, 1.91433, 2.54746, 3.28517, 4.13277,
00648  5.09574, 6.17984, 7.39106, 8.73568, 10.2203, 11.8519, 13.6377, 15.5854, 17.7033, 20.0000};
00649   static const G4double I9[nX]={0,
00650  .030977, .125727, .289605, .528146, .846967, 1.25183, 1.74871, 2.34384, 3.04376, 3.85535,
00651  4.78594, 5.84329, 7.03567, 8.37194, 9.86163, 11.5150, 13.3430, 15.3576, 17.5719, 20.0000};
00652   static const G4double IA[nX]={0,
00653  .027129, .111420, .258935, .475812, .768320, 1.14297, 1.60661, 2.16648, 2.83034, 3.60650,
00654  4.50394, 5.53238, 6.70244, 8.02569, 9.51488, 11.1841, 13.0488, 15.1264, 17.4362, 20.0000};
00655   static const G4double IB[nX]={0,
00656  .024170, .100153, .234345, .433198, .703363, 1.05184, 1.48607, 2.01409, 2.64459, 3.38708,
00657  4.25198, 5.25084, 6.39647, 7.70319, 9.18708, 10.8663, 12.7617, 14.8968, 17.2990, 20.0000};
00658   static const G4double IC[nX]={0,
00659  .021877, .091263, .214670, .398677, .650133, .976322, 1.38510, 1.88504, 2.48555, 3.19709,
00660  4.03129, 5.00127, 6.12184, 7.40989, 8.88482, 10.5690, 12.4888, 14.6748, 17.1638, 20.0000};
00661   static const G4double ID[nX]={0,
00662  .020062, .084127, .198702, .370384, .606100, .913288, 1.30006, 1.77535, 2.34912, 3.03253,
00663  3.83822, 4.78063, 5.87634, 7.14459, 8.60791, 10.2929, 12.2315, 14.4621, 17.0320, 20.0000};
00664   static const G4double IE[nX]={0,
00665  .018547, .078104, .185102, .346090, .567998, .858331, 1.22535, 1.67824, 2.22735, 2.88443,
00666  3.66294, 4.57845, 5.64911, 6.89637, 8.34578, 10.0282, 11.9812, 14.2519, 16.8993, 20.0000};
00667   static const G4double IF[nX]={0,
00668  .017143, .072466, .172271, .323007, .531545, .805393, 1.15288, 1.58338, 2.10754, 2.73758,
00669  3.48769, 4.37450, 5.41770, 6.64092, 8.07288, 9.74894, 11.7135, 14.0232, 16.7522, 20.0000};
00670   static const G4double IG[nX]={0,
00671  .015618, .066285, .158094, .297316, .490692, .745653, 1.07053, 1.47479, 1.96931, 2.56677,
00672  3.28205, 4.13289, 5.14068, 6.33158, 7.73808, 9.40133, 11.3745, 13.7279, 16.5577, 20.0000};
00673   static const G4double IH[nX]={0,
00674  .013702, .058434, .139923, .264115, .437466, .667179, .961433, 1.32965, 1.78283, 2.33399,
00675  2.99871, 3.79596, 4.74916, 5.88771, 7.24937, 8.88367, 10.8576, 13.2646, 16.2417, 20.0000};
00676   static const G4double II[nX]={0,
00677  .011264, .048311, .116235, .220381, .366634, .561656, .813132, 1.13008, 1.52322, 2.00554,
00678  2.59296, 3.30542, 4.16834, 5.21490, 6.48964, 8.05434, 9.99835, 12.4580, 15.6567, 20.0000};
00679   static const G4double IJ[nX]={0,
00680  .008628, .037206, .089928, .171242, .286114, .440251, .640343, .894382, 1.21208, 1.60544,
00681  2.08962, 2.68414, 3.41486, 4.31700, 5.44048, 6.85936, 8.69067, 11.1358, 14.5885, 20.0000};
00682   static const G4double* Il[nE]=
00683                              {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
00684   static const G4double lE[nE]={
00685 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
00686  2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
00687   static const G4double lEmi=lE[0];
00688   static const G4double lEma=lE[nE-1];
00689   static const G4double dlE=(lEma-lEmi)/bE;
00690   //***************************************************************************************
00691   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00692   G4double lEn=std::log(Enu);         // log(E) for interpolation
00693   G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
00694   G4int fE=static_cast<int>(rE);      // Left bin for interpolation
00695   if(fE<0) fE=0;
00696   if(fE>pE)fE=pE;
00697   G4int    sE=fE+1;                   // Right bin for interpolation
00698   G4double dE=rE-fE;                  // relative log shift from the left bin
00699   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00700   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00701   G4double Emu=Enu-mmu;               // Free Energy of neutrino/anti-neutrino
00702   G4double ME=Enu*MN;                 // M*E
00703   G4double dME=ME+ME;                 // 2*M*E
00704   G4double dEMN=(dEnu+MN)*ME;
00705   G4double MEm=ME-hmmu2;
00706   G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
00707   G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
00708   G4double ymax=(E2M+sqE)/dEMN;
00709   G4double ymin=(E2M-sqE)/dEMN;
00710   G4double rmin=1.-ymin;
00711   G4double rhm2E=hmmu2/Enu2;
00712   G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
00713   G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
00714   G4double Q2nq=Emu*dMN-mcV;
00715   if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
00716   // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
00717   G4double Rmi=Q2mi/Q2ma;
00718   G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78);
00719   // --- E-interpolation must be done in a log scale ---
00720   G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
00721   G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
00722   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00723   G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
00724   G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
00725   G4double rXi=(Xmi-iXmi)/idX;
00726   G4int    iXi=static_cast<int>(rXi);
00727   if(iXi<0) iXi=0;
00728   if(iXi>bX) iXi=bX;
00729   G4double dXi=rXi-iXi;
00730   G4double bntil=Il[fE][iXi];
00731   G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
00732   G4double bntir=Il[sE][iXi];
00733   G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
00734   G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
00735   //
00736   G4double rXa=(Xma-iXmi)/idX;
00737   G4int    iXa=static_cast<int>(rXa);
00738   if(iXa<0) iXa=0;
00739   if(iXa>bX) iXa=bX;
00740   G4double dXa=rXa-iXa;
00741   G4double bntal=Il[fE][iXa];
00742   G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
00743   G4double bntar=Il[sE][iXa];
00744   G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
00745   G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
00746   //
00747   // *** Find X using the reversed table ***
00748   G4double intx=inti+(inta-inti)*G4UniformRand(); 
00749   G4int    intc=static_cast<int>(intx);
00750   if(intc<0) intc=0;
00751   if(intc>bX) intc=bX;
00752   G4double dint=intx-intc;
00753   G4double mXl=Xl[fE][intc];
00754   G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
00755   G4double mXr=Xl[sE][intc];
00756   G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
00757   G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
00758   G4double R=shift-std::pow(X,pconv);
00759   G4double Q2=R*Q2ma;
00760   return Q2*GeV*GeV;
00761 }
00762 
00763 // It returns a fraction of the direct interaction of the neutrino with quark-partons
00764 G4double G4QANuMuNuclearCrossSection::GetDirectPart(G4double Q2)
00765 {
00766   G4double f=Q2/4.62;
00767   G4double ff=f*f;
00768   G4double r=ff*ff;
00769   G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
00770   //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
00771   return 1.-s_value*(1.-s_value/2);
00772 }
00773 
00774 // #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
00775 G4double G4QANuMuNuclearCrossSection::GetNPartons(G4double Q2)
00776 {
00777   return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
00778 }
00779 
00780 // This class can provide only virtual exchange pi- (a substitute for W- boson)
00781 G4int G4QANuMuNuclearCrossSection::GetExchangePDGCode() {return -211;}

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