//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: G4QANuMuNuclearCrossSection.cc,v 1.11 2007/11/01 16:09:38 mkossov Exp $
// GEANT4 tag $Name: geant4-09-02-ref-02 $
//
//
// G4 Physics class: G4QANuMuNuclearCrossSection for gamma+A cross sections
// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
// The last update: M.V. Kossov, CERN/ITEP (Moscow) 17-Oct-03
// 
// ****************************************************************************************
// ********** This CLASS is temporary moved from the photolepton_hadron directory *********
// ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ******
// ****************************************************************************************
//===============================================================================================

//#define debug
//#define edebug
//#define pdebug
//#define ppdebug
//#define tdebug
//#define sdebug

#include "G4QANuMuNuclearCrossSection.hh"

// Initialization of the
G4bool    G4QANuMuNuclearCrossSection::onlyCS=true;//Flag to calculate only CS (not QE)
G4double  G4QANuMuNuclearCrossSection::lastSig=0.;//Last calculated total cross section
G4double  G4QANuMuNuclearCrossSection::lastQEL=0.;//Last calculated quasi-el. cross section
G4int     G4QANuMuNuclearCrossSection::lastL=0;   //Last used in cross section TheLastBin
G4double  G4QANuMuNuclearCrossSection::lastE=0.;  //Last used in cross section TheEnergy
G4double* G4QANuMuNuclearCrossSection::lastEN=0;  //Pointer to the Energy Scale of TX & QE
G4double* G4QANuMuNuclearCrossSection::lastTX=0;  //Pointer to the LastArray of TX function
G4double* G4QANuMuNuclearCrossSection::lastQE=0;  //Pointer to the LastArray of QE function
G4int     G4QANuMuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
G4int     G4QANuMuNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
G4int     G4QANuMuNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
G4double  G4QANuMuNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
G4double  G4QANuMuNuclearCrossSection::lastTH=0.; // Last threshold momentum
G4double  G4QANuMuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
G4int     G4QANuMuNuclearCrossSection::lastI=0;   // The last position in the DAMDB

// Returns Pointer to the G4VQCrossSection class
G4VQCrossSection* G4QANuMuNuclearCrossSection::GetPointer()
{
  static G4QANuMuNuclearCrossSection theCrossSection;//**Static body of the Cross Section**
  return &theCrossSection;
}

// The main member function giving the collision cross section (P is in IU, CS is in mb)
// Make pMom in independent units ! (Now it is MeV)
G4double G4QANuMuNuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
                                                      G4int tgZ, G4int tgN, G4int pPDG)
{
  static G4int j;                      // A#0f records found in DB for this projectile
  static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
  static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
  static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
  static std::vector <G4double> colP;  // Vector of last momenta for the reaction
  static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
  static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
  // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
  G4double pEn=pMom;
#ifdef debug
  G4cout<<"G4QAMNCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
        <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
        <<colN.size()<<G4endl;
		//CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
  if(pPDG!=-14)
  {
#ifdef debug
    G4cout<<"G4QAMNCS::GetCS: *** Found pPDG="<<pPDG<<" ====> CS=0"<<G4endl;
    //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
    return 0.;                         // projectile PDG=0 is a mistake (?!) @@
  }
  G4bool in=false;                     // By default the isotope must be found in the AMDB
  if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
  {
    in = false;                        // By default the isotope haven't be found in AMDB  
    lastP   = 0.;                      // New momentum history (nothing to compare with)
    lastPDG = pPDG;                    // The last PDG of the projectile
    lastN   = tgN;                     // The last N of the calculated nucleus
    lastZ   = tgZ;                     // The last Z of the calculated nucleus
    lastI   = colN.size();             // Size of the Associative Memory DB in the heap
    j  = 0;                            // A#0f records found in DB for this projectile
    if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found
	   {                                  // The nucleus with projPDG is found in AMDB
      if(colN[i]==tgN && colZ[i]==tgZ)
						{
        lastI=i;
        lastTH =colTH[i];                // Last THreshold (A-dependent)
#ifdef pdebug
        G4cout<<"G4QAMNCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
        //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
        if(pEn<=lastTH)
        {
#ifdef pdebug
          G4cout<<"G4QAMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl;
          //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
          return 0.;                     // Energy is below the Threshold value
        }
        lastP  =colP [i];                // Last Momentum  (A-dependent)
        lastCS =colCS[i];                // Last CrossSect (A-dependent)
        if(std::fabs(lastP/pMom-1.)<tolerance)
        {
#ifdef pdebug
          G4cout<<"G4QAMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
          //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
          return lastCS*millibarn;     // Use theLastCS
        }
        in = true;                       // This is the case when the isotop is found in DB
        // Momentum pMom is in IU ! @@ Units
#ifdef pdebug
        G4cout<<"G4QAMNCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
#endif
        lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
#ifdef pdebug
        G4cout<<"G4QAMNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
        //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
        if(lastCS<=0. && pEn>lastTH)    // Correct the threshold
        {
#ifdef pdebug
          G4cout<<"G4QAMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pEn;
        }
        break;                           // Go out of the LOOP
      }
#ifdef pdebug
      G4cout<<"---G4QAMNCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
            <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl;
      //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
      j++;                             // Increment a#0f records found in DB for this pPDG
	   }
	   if(!in)                            // This nucleus has not been calculated previously
	   {
#ifdef pdebug
      G4cout<<"G4QAMNCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<lastI<<G4endl;
#endif
      //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
      lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
      if(lastCS<=0.)
						{
        lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
#ifdef pdebug
        G4cout<<"G4QAMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
#endif
        if(pEn>lastTH)
        {
#ifdef pdebug
          G4cout<<"G4QAMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pEn;
        }
						}
#ifdef pdebug
      G4cout<<"G4QAMNCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
      //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
      colN.push_back(tgN);
      colZ.push_back(tgZ);
      colPDG.push_back(pPDG);
      colP.push_back(pMom);
      colTH.push_back(lastTH);
      colCS.push_back(lastCS);
#ifdef pdebug
      G4cout<<"G4QAMNCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl;
      //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
      return lastCS*millibarn;
	   } // End of creation of the new set of parameters
    else
				{
#ifdef pdebug
      G4cout<<"G4QAMNCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
#endif
      colP[lastI]=pMom;
      colPDG[lastI]=pPDG;
      colCS[lastI]=lastCS;
    }
  } // End of parameters udate
  else if(pEn<=lastTH)
  {
#ifdef pdebug
    G4cout<<"G4QAMNCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
    //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
    return 0.;                         // Momentum is below the Threshold Value -> CS=0
  }
  else if(std::fabs(lastP/pMom-1.)<tolerance)
  {
#ifdef pdebug
    G4cout<<"G4QAMNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
    //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
    return lastCS*millibarn;     // Use theLastCS
  }
  else
  {
#ifdef pdebug
    G4cout<<"G4QAMNCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
#endif
    lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
    lastP=pMom;
  }
#ifdef pdebug
  G4cout<<"G4QAMNCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
  //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
  return lastCS*millibarn;
}

// Gives the threshold energy = the same for all nuclei (@@ can be reduced for hevy nuclei)
G4double G4QANuMuNuclearCrossSection::ThresholdEnergy(G4int Z, G4int N, G4int)
{
  //static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1.,0.)/GeV;
  //static const G4double mProt = G4NucleiProperties::GetNuclearMass(1.,1.)/GeV;
  //static const G4double mDeut = G4NucleiProperties::GetNuclearMass(2.,1.)/GeV/2.;
  static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
  static const G4double dmN=mN+mN;    // Doubled nucleon mass (2*AtomicMassUnit, GeV)
  static const G4double mmu=.105658369; // Mass of a muon in GeV
  static const G4double mmu2=mmu*mmu;   // Squared mass of a muon in GeV^2
  static const G4double thresh=mmu+mmu2/dmN; // Universal threshold in GeV
  // ---------
  //static const G4double infEn = 9.e27;
  G4double dN=0.;
  if(Z>0||N>0) dN=thresh*GeV; // @@ if upgraded, change it in a total cross section
  //@@ "dN=mmu+mmu2/G4NucleiProperties::GetNuclearMass(<G4double>(Z+N),<G4double>(Z)/GeV"
  return dN;
}

// The main member function giving the gamma-A cross section (E_kin in MeV, CS in mb)
G4double G4QANuMuNuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I,
                                       G4int , G4int targZ, G4int targN, G4double Momentum)
{
  static const G4double mb38=1.E-11;// Conversion 10^-38 cm^2 to mb=10^-27 cm^2
  static const G4int nE=65;   // !! If change this, change it in GetFunctions() (*.hh) !!
  static const G4int mL=nE-1;
  static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
  static const G4double dmN=mN+mN;      // Doubled nucleon mass (2*AtomicMassUnit, GeV)
  static const G4double mmu=.105658369; // Mass of a muon in GeV
  static const G4double mmu2=mmu*mmu;   // Squared mass of a muon in GeV^2
  static const G4double EMi=mmu+mmu2/dmN; // Universal threshold of the reaction in GeV
  static const G4double EMa=300.;       // Maximum tabulated Energy of nu_mu in GeV 
  // *** Begin of the Associative memory for acceleration of the cross section calculations
  static std::vector <G4double> colH;   //?? Vector of HighEnergyCoefficients (functional)
  static std::vector <G4double*> TX;    // Vector of pointers to the TX tabulated functions
  static std::vector <G4double*> QE;    // Vector of pointers to the QE tabulated functions
  static G4bool first=true;             // Flag of initialization of the energy axis
  // *** End of Static Definitions (Associative Memory) ***
  //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
  //G4double TotEnergy2=Momentum;
  onlyCS=CS;                            // Flag to calculate only CS (not TX & QE)
  lastE=Momentum/GeV;                   // Kinetic energy of the muon neutrino (in GeV!)
  if (lastE<=EMi)                       // Energy is below the minimum energy in the table
  {
    lastE=0.;
    lastSig=0.;
    return 0.;
  }
  G4int Z=targZ;                        // New Z, which can change the sign
  if(F<=0)                              // This isotope was not the last used isotop
  {
    if(F<0)                          // This isotope was found in DAMDB =========> RETRIEVE
				{
      lastTX =TX[I];                 // Pointer to the prepared TX function (same isotope)
      lastQE =QE[I];                 // Pointer to the prepared QE function (same isotope)
	  }
	  else                              // This isotope wasn't calculated previously => CREATE
	  {
      if(first)
      {
        lastEN = new G4double[nE];   // This must be done only once!
        Z=-Z;                        // To explain GetFunctions that E-axis must be filled
        first=false;                 // To make it only once
      }
      lastTX = new G4double[nE];     // Allocate memory for the new TX function
      lastQE = new G4double[nE];     // Allocate memory for the new QE function
      G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
      if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
      // *** The synchronization check ***
      G4int sync=TX.size();
      if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
      TX.push_back(lastTX);
      QE.push_back(lastQE);
	   } // End of creation of the new set of parameters
  } // End of parameters udate
  // ============================== NOW Calculate the Cross Section =====================
  if (lastE<=EMi)                    // Check that antiNuEnergy is higher than ThreshE
  {
    lastE=0.;
    lastSig=0.;
    return 0.;
  }
  if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
  {
    G4int chk=1;
    G4int ran=mL/2;
    G4int sep=ran;  // as a result = an index of the left edge of the interval
    while(ran>=2)
				{
      G4int newran=ran/2;
      if(lastE<=lastEN[sep]) sep-=newran;
      else                   sep+=newran;
      ran=newran;
      chk=chk+chk; 
    }
    if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
    G4double lowE=lastEN[sep];
    G4double highE=lastEN[sep+1];
    G4double lowTX=lastTX[sep];
    if(lastE<lowE||sep>=mL||lastE>highE)
      G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
            <<", sep="<<sep<<", mL="<<mL<<G4endl;
    lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
    if(!onlyCS)                       // Skip the differential cross-section parameters
    {
      G4double lowQE=lastQE[sep];
      lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
#ifdef pdebug
      G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
#endif
    }
  }
  else
  {
    lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
    lastQEL=lastQE[mL];
  }
  if(lastQEL<0.) lastQEL = 0.;
  if(lastSig<0.) lastSig = 0.;
  // The cross-sections are expected to be in mb
  lastSig*=mb38;
  if(!onlyCS) lastQEL*=mb38;
  return lastSig;
}

// Calculate the cros-section functions
// ****************************************************************************************
// *** This tables are the same for all lepto-nuclear reactions, only mass is different ***
// ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ ***
// ****************************************************************************************
G4int G4QANuMuNuclearCrossSection::GetFunctions (G4int z, G4int n,
                                                     G4double* t, G4double* q, G4double* e)
{
  static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
  static const G4double dmN=mN+mN;    // Doubled nucleon mass (2*AtomicMassUnit, GeV)
  static const G4double mmu=.105658369; // Mass of a muon in GeV
  static const G4double mmu2=mmu*mmu;   // Squared mass of a muon in GeV^2
  static const G4double thresh=mmu+mmu2/dmN; // Universal threshold in GeV
  static const G4int nE=65; // !! If change this, change it in GetCrossSection() (*.cc) !!
  static const G4double nuEn[nE]={thresh,
.112039,.116079,.120416,.125076,.130090,.135494,.141324,.147626,.154445,.161838,
.169864,.178594,.188105,.198485,.209836,.222272,.235923,.250941,.267497,.285789,
.306045,.328530,.353552,.381466,.412689,.447710,.487101,.531538,.581820,.638893,
.703886,.778147,.863293,.961275,1.07445,1.20567,1.35843,1.53701,1.74667,1.99390,
2.28679,2.63542,3.05245,3.55386,4.15990,4.89644,5.79665,6.90336,8.27224,9.97606,
12.1106,14.8029,18.2223,22.5968,28.2351,35.5587,45.1481,57.8086,74.6682,97.3201,
128.036,170.085,228.220,309.420};
  static const G4double TOTX[nE]={0.,
.077498,.247583,.329691,.386384,.429087,.462699,.489899,.512316,.530996,.546614,
.559616,.570292,.578840,.585395,.590053,.593083,.594197,.593614,.591396,.587611,
.582335,.575653,.567667,.558490,.548417,.537270,.525352,.512825,.499857,.486620,
.473283,.460014,.446970,.434294,.422116,.410656,.399782,.389665,.380349,.371860,
.364207,.357387,.351388,.346192,.341778,.338122,.335198,.332980,.331439,.330544,
.330263,.330558,.331391,.332718,.334494,.336667,.339182,.341697,.344470,.348125,
.351322,.354481,.357507,.359239};
  static const G4double QELX[nE]={0.,
.008683,.028739,.039700,.048327,.055820,.062693,.069235,.075631,.082010,.088463,
.095059,.101851,.108883,.116192,.123814,.131826,.140185,.148962,.158197,.167933,
.178221,.189119,.200700,.213045,.226326,.240454,.255277,.270612,.286388,.302608,
.319318,.336582,.354468,.373031,.392427,.412445,.433146,.454448,.476222,.498289,
.520430,.542558,.564130,.585003,.604928,.623680,.641266,.657255,.671704,.684586,
.696111,.706028,.714553,.721951,.728085,.733182,.737348,.740958,.743716,.746059,
.747806,.749129,.750331,.751100};

  // --------------------------------
  G4int first=0;
  if(z<0.)
		{
    first=1;
    z=-z;
  }
  if(z<1 || z>92)             // neutron & plutonium are forbidden
  {
    G4cout<<"**G4QANuMuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
    return -1;
  }
  for(G4int k=0; k<nE; k++)
  {
    G4double a=n+z;
    G4double za=z+a;
    G4double dz=z+z;
    G4double da=a+a;
    G4double ta=da+a;
    if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
    t[k]=TOTX[k]*nuEn[k]*(za+za)/ta+QELX[k]*(dz+dz-da)/ta; // TotalCrossSection
    q[k]=QELX[k]*dz/a;                                     // QuasiElasticCrossSection
	 }
  return first;
}

// Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
G4double G4QANuMuNuclearCrossSection::GetQEL_ExchangeQ2()
{
  static const G4double mmu=.105658369;// Mass of muon in GeV
  static const G4double mmu2=mmu*mmu;  // Squared Mass of muon in GeV^2
  static const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
  static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
  static const double MN2=MN*MN;       // M_N^2 in GeV^2
  static const G4double power=-3.5;    // direct power for the magic variable
  static const G4double pconv=1./power;// conversion power for the magic variable
  static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
  static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
  static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
  // Reversed table
  static const G4double Xl[nQ2]={5.20224e-16,
 .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898,
	.093878,	.103332,	.112844,	.122410,	.132023,	.141680,	.151376,	.161109,	.170875,	.180672,
	.190499,	.200352,	.210230,	.220131,	.230055,	.239999,	.249963,	.259945,	.269944,	.279960,
	.289992,	.300039,	.310099,	.320173,	.330260,	.340359,	.350470,	.360592,	.370724,	.380867,
	.391019,	.401181,	.411352,	.421531,	.431719,	.441915,	.452118,	.462329,	.472547,	.482771,
	.493003,	.503240,	.513484,	.523734,	.533989,	.544250,	.554517,	.564788,	.575065,	.585346,
	.595632,	.605923,	.616218,	.626517,	.636820,	.647127,	.657438,	.667753,	.678072,	.688394,
	.698719,	.709048,	.719380,	.729715,	.740053,	.750394,	.760738,	.771085,	.781434,	.791786,
	.802140,	.812497,	.822857,	.833219,	.843582,	.853949,	.864317,	.874687,	.885060,	.895434,
	.905810,	.916188,	.926568,	.936950,	.947333,	.957719,	.968105,	.978493,	.988883,	.999275};
	  // Direct table
  static const G4double Xmax=Xl[lQ2];
  static const G4double Xmin=Xl[0];
  static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
  static const G4double inl[nQ2]={0,
	1.52225,	2.77846,	3.96651,	5.11612,	6.23990,	7.34467,	8.43466,	9.51272,	10.5809,	11.6406,
	12.6932,	13.7394,	14.7801,	15.8158,	16.8471,	17.8743,	18.8979,	19.9181,	20.9353,	21.9496,
	22.9614,	23.9707,	24.9777,	25.9826,	26.9855,	27.9866,	28.9860,	29.9837,	30.9798,	31.9745,
	32.9678,	33.9598,	34.9505,	35.9400,	36.9284,	37.9158,	38.9021,	39.8874,	40.8718,	41.8553,
	42.8379,	43.8197,	44.8007,	45.7810,	46.7605,	47.7393,	48.7174,	49.6950,	50.6718,	51.6481,
	52.6238,	53.5990,	54.5736,	55.5476,	56.5212,	57.4943,	58.4670,	59.4391,	60.4109,	61.3822,
	62.3531,	63.3236,	64.2937,	65.2635,	66.2329,	67.2019,	68.1707,	69.1390,	70.1071,	71.0748,
	72.0423,	73.0095,	73.9763,	74.9429,	75.9093,	76.8754,	77.8412,	78.8068,	79.7721,	80.7373,
	81.7022,	82.6668,	83.6313,	84.5956,	85.5596,	86.5235,	87.4872,	88.4507,	89.4140,	90.3771,
	91.3401,	92.3029,	93.2656,	94.2281,	95.1904,	96.1526,	97.1147,	98.0766,	99.0384,	100.000};
  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
  G4double ME=Enu*MN;                 // M*E
  G4double dME=ME+ME;                 // 2*M*E
  G4double dEMN=(dEnu+MN)*ME;
  G4double MEm=ME-hmmu2;
  G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
  G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
  G4double ymax=(E2M+sqE)/dEMN;
  G4double ymin=(E2M-sqE)/dEMN;
  G4double rmin=1.-ymin;
  G4double rhm2E=hmmu2/Enu2;
  G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
  G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
  G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
  G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
  G4double rXi=(Xmi-Xmin)/dX;
  G4int    iXi=static_cast<int>(rXi);
  if(iXi<0) iXi=0;
  if(iXi>bQ2) iXi=bQ2;
  G4double dXi=rXi-iXi;
  G4double bnti=inl[iXi];
  G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
  //
  G4double rXa=(Xma-Xmin)/dX;
  G4int    iXa=static_cast<int>(rXa);
  if(iXa<0) iXa=0;
  if(iXa>bQ2) iXa=bQ2;
  G4double dXa=rXa-iXa;
  G4double bnta=inl[iXa];
  G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
  // *** Find X using the reversed table ***
  G4double intx=inti+(inta-inti)*G4UniformRand();
  G4int    intc=static_cast<int>(intx);
  if(intc<0) intc=0;
  if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
  G4double dint=intx-intc;
  G4double mX=Xl[intc];
  G4double X=mX+dint*(Xl[intc+1]-mX);
  G4double Q2=std::pow(X,pconv)-1.;
  return Q2*GeV*GeV;
}

// Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
G4double G4QANuMuNuclearCrossSection::GetNQE_ExchangeQ2()
{
  static const double mpi=.13957018;    // charged pi meson mass in GeV
  static const G4double mmu=.105658369; // Mass of muon in GeV
  static const G4double mmu2=mmu*mmu;   // Squared Mass of muon in GeV^2
  static const double hmmu2=mmu2/2;     // .5*m_mu^2 in GeV^2
  static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
  static const double MN2=MN*MN;        // M_N^2 in GeV^2
  static const double dMN=MN+MN;        // 2*M_N in GeV
  static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
  static const G4int power=7;           // direct power for the magic variable
  static const G4double pconv=1./power; // conversion power for the magic variable
  static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
  static const G4int lX=nX-1;           // index of the last in the Xl table
  static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
  static const G4int nE=20;             // #Of point in the El table (in GeV^2)
  static const G4int bE=nE-1;           // index of the last in the El table
  static const G4int pE=bE-1;           // index of the before last in the El table
  // Reversed table
  static const G4double X0[nX]={5.21412e-05,
	.437860,	.681908,	.891529,	1.08434,	1.26751,	1.44494,	1.61915,	1.79198,	1.96493,	2.13937,
	2.31664,	2.49816,	2.68559,	2.88097,	3.08705,	3.30774,	3.54917,	3.82233,	4.15131, 4.62182};
  static const G4double X1[nX]={.00102591,
	1.00443,	1.55828,	2.03126,	2.46406,	2.87311,	3.26723,	3.65199,	4.03134,	4.40835,	4.78561,
	5.16549,	5.55031,	5.94252,	6.34484,	6.76049,	7.19349,	7.64917,	8.13502,	8.66246, 9.25086};
  static const G4double X2[nX]={.0120304,
	2.59903,	3.98637,	5.15131,	6.20159,	7.18024,	8.10986,	9.00426,	9.87265,	10.7217,	11.5564,
	12.3808,	13.1983,	14.0116,	14.8234,	15.6359,	16.4515,	17.2723,	18.1006,	18.9386, 19.7892};
  static const G4double X3[nX]={.060124,
	5.73857,	8.62595,	10.9849,	13.0644,	14.9636,	16.7340,	18.4066,	20.0019,	21.5342,	23.0142,
	24.4497,	25.8471,	27.2114,	28.5467,	29.8564,	31.1434,	32.4102,	33.6589,	34.8912, 36.1095};
  static const G4double X4[nX]={.0992363,
	8.23746,	12.1036,	15.1740,	17.8231,	20.1992,	22.3792,	24.4092,	26.3198,	28.1320,	29.8615,
	31.5200,	33.1169,	34.6594,	36.1536,	37.6044,	39.0160,	40.3920,	41.7353,	43.0485, 44.3354};
  static const G4double X5[nX]={.0561127,
	7.33661,	10.5694,	13.0778,	15.2061,	17.0893,	18.7973,	20.3717,	21.8400,	23.2211,	24.5291,
	25.7745,	26.9655,	28.1087,	29.2094,	30.2721,	31.3003,	32.2972,	33.2656,	34.2076, 35.1265};
  static const G4double X6[nX]={.0145859,
	4.81774,	6.83565,	8.37399,	9.66291,	10.7920,	11.8075,	12.7366,	13.5975,	14.4025,	15.1608,
	15.8791,	16.5628,	17.2162,	17.8427,	18.4451,	19.0259,	19.5869,	20.1300,	20.6566, 21.1706};
  static const G4double X7[nX]={.00241155,
	2.87095,	4.02492,	4.89243,	5.61207,	6.23747,	6.79613,	7.30433,	7.77270,	8.20858,	8.61732,
	9.00296,	9.36863,	9.71682,	10.0495,	10.3684,	10.6749,	10.9701,	11.2550,	11.5306, 11.7982};
  static const G4double X8[nX]={.000316863,
	1.76189,	2.44632,	2.95477,	3.37292,	3.73378,	4.05420,	4.34415,	4.61009,	4.85651,	5.08666,
	5.30299,	5.50738,	5.70134,	5.88609,	6.06262,	6.23178,	6.39425,	6.55065,	6.70149, 6.84742};
  static const G4double X9[nX]={3.73544e-05,
	1.17106,	1.61289,	1.93763,	2.20259,	2.42976,	2.63034,	2.81094,	2.97582,	3.12796,	3.26949,
	3.40202,	3.52680,	3.64482,	3.75687,	3.86360,	3.96557,	4.06323,	4.15697,	4.24713, 4.33413};
  static const G4double XA[nX]={4.19131e-06,
	.849573,	1.16208,	1.38955,	1.57379,	1.73079,	1.86867,	1.99221,	2.10451,	2.20770,	2.30332,
	2.39252,	2.47622,	2.55511,	2.62977,	2.70066,	2.76818,	2.83265,	2.89437,	2.95355, 3.01051};
  static const G4double XB[nX]={4.59981e-07,
	.666131,	.905836,	1.07880,	1.21796,	1.33587,	1.43890,	1.53080,	1.61399,	1.69011,	1.76040,
	1.82573,	1.88682,	1.94421,	1.99834,	2.04959,	2.09824,	2.14457,	2.18878,	2.23107, 2.27162};
  static const G4double XC[nX]={4.99861e-08,
	.556280,	.752730,	.893387,	1.00587,	1.10070,	1.18317,	1.25643,	1.32247,	1.38269,	1.43809,
	1.48941,	1.53724,	1.58203,	1.62416,	1.66391,	1.70155,	1.73728,	1.77128,	1.80371, 1.83473};
  static const G4double XD[nX]={5.40832e-09,
	.488069,	.657650,	.778236,	.874148,	.954621,	1.02432,	1.08599,	1.14138,	1.19172,	1.23787,
	1.28049,	1.32008,	1.35705,	1.39172,	1.42434,	1.45514,	1.48429,	1.51197,	1.53829, 1.56339};
  static const G4double XE[nX]={5.84029e-10,
	.445057,	.597434,	.705099,	.790298,	.861468,	.922865,	.976982,	1.02542,	1.06930,	1.10939,
	1.14630,	1.18050,	1.21233,	1.24208,	1.27001,	1.29630,	1.32113,	1.34462,	1.36691, 1.38812};
  static const G4double XF[nX]={6.30137e-11,
	.418735,	.560003,	.659168,	.737230,	.802138,	.857898,	.906854,	.950515,	.989915,	1.02580,
	1.05873,	1.08913,	1.11734,	1.14364,	1.16824,	1.19133,	1.21306,	1.23358,	1.25298, 1.27139};
  static const G4double XG[nX]={6.79627e-12,
	.405286,	.539651,	.633227,	.706417,	.766929,	.818642,	.863824,	.903931,	.939963,	.972639,
	1.00250,	1.02995,	1.05532,	1.07887,	1.10082,	1.12134,	1.14058,	1.15867,	1.17572, 1.19183};
  static const G4double XH[nX]={7.32882e-13,
	.404391,	.535199,	.625259,	.695036,	.752243,	.800752,	.842823,	.879906,	.912994,	.942802,
	.969862,	.994583,	1.01729,	1.03823,	1.05763,	1.07566,	1.09246,	1.10816,	1.12286, 1.13667};
  static const G4double XI[nX]={7.90251e-14,
	.418084,	.548382,	.636489,	.703728,	.758106,	.803630,	.842633,	.876608,	.906576,	.933269,
	.957233,	.978886,	.998556,	1.01651,	1.03295,	1.04807,	1.06201,	1.07489,	1.08683, 1.09792};
  static const G4double XJ[nX]={8.52083e-15,
	.447299,	.579635,	.666780,	.731788,	.783268,	.825512,	.861013,	.891356,	.917626,	.940597,
	.960842,	.978802,	.994820,	1.00917,	1.02208,	1.03373,	1.04427,	1.05383,	1.06253, 1.07046};
  // Direct table
  static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
                        X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
  static const G4double dX[nE]={
    (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
    (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
    (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
    (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
    (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
  static const G4double* Xl[nE]=
                             {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
  static const G4double I0[nX]={0,
	.354631,	1.08972,	2.05138,	3.16564,	4.38343,	5.66828,	6.99127,	8.32858,	9.65998,	10.9680,
	12.2371,	13.4536,	14.6050,	15.6802,	16.6686,	17.5609,	18.3482,	19.0221,	19.5752,	20.0000};
  static const G4double I1[nX]={0,
	.281625,	.877354,	1.67084,	2.60566,	3.64420,	4.75838,	5.92589,	7.12829,	8.34989,	9.57708,
	10.7978,	12.0014,	13.1781,	14.3190,	15.4162,	16.4620,	17.4496,	18.3724,	19.2245,	20.0000};
  static const G4double I2[nX]={0,
	.201909,	.642991,	1.24946,	1.98463,	2.82370,	3.74802,	4.74263,	5.79509,	6.89474,	8.03228,
	9.19947,	10.3889,	11.5938,	12.8082,	14.0262,	15.2427,	16.4527,	17.6518,	18.8356,	20.0000};
  static const G4double I3[nX]={0,
	.140937,	.461189,	.920216,	1.49706,	2.17728,	2.94985,	3.80580,	4.73758,	5.73867,	6.80331,
	7.92637,	9.10316,	10.3294,	11.6013,	12.9150,	14.2672,	15.6548,	17.0746,	18.5239,	20.0000};
  static const G4double I4[nX]={0,
	.099161,	.337358,	.694560,	1.16037,	1.72761,	2.39078,	3.14540,	3.98768,	4.91433,	5.92245,
	7.00942,	8.17287,	9.41060,	10.7206,	12.1010,	13.5500,	15.0659,	16.6472,	18.2924,	20.0000};
  static const G4double I5[nX]={0,
	.071131,	.255084,	.543312,	.932025,	1.41892,	2.00243,	2.68144,	3.45512,	4.32283,	5.28411,
	6.33859,	7.48602,	8.72621,	10.0590,	11.4844,	13.0023,	14.6128,	16.3158,	18.1115,	20.0000};
  static const G4double I6[nX]={0,
	.053692,	.202354,	.443946,	.778765,	1.20774,	1.73208,	2.35319,	3.07256,	3.89177,	4.81249,
	5.83641,	6.96528,	8.20092,	9.54516,	10.9999,	12.5670,	14.2486,	16.0466,	17.9630,	20.0000};
  static const G4double I7[nX]={0,
	.043065,	.168099,	.376879,	.672273,	1.05738,	1.53543,	2.10973,	2.78364,	3.56065,	4.44429,
	5.43819,	6.54610,	7.77186,	9.11940,	10.5928,	12.1963,	13.9342,	15.8110,	17.8313,	20.0000};
  static const G4double I8[nX]={0,
	.036051,	.143997,	.327877,	.592202,	.941572,	1.38068,	1.91433,	2.54746,	3.28517,	4.13277,
	5.09574,	6.17984,	7.39106,	8.73568,	10.2203,	11.8519,	13.6377,	15.5854,	17.7033,	20.0000};
  static const G4double I9[nX]={0,
	.030977,	.125727,	.289605,	.528146,	.846967,	1.25183,	1.74871,	2.34384,	3.04376,	3.85535,
	4.78594,	5.84329,	7.03567,	8.37194,	9.86163,	11.5150,	13.3430,	15.3576,	17.5719,	20.0000};
  static const G4double IA[nX]={0,
	.027129,	.111420,	.258935,	.475812,	.768320,	1.14297,	1.60661,	2.16648,	2.83034,	3.60650,
	4.50394,	5.53238,	6.70244,	8.02569,	9.51488,	11.1841,	13.0488,	15.1264,	17.4362,	20.0000};
  static const G4double IB[nX]={0,
	.024170,	.100153,	.234345,	.433198,	.703363,	1.05184,	1.48607,	2.01409,	2.64459,	3.38708,
	4.25198,	5.25084,	6.39647,	7.70319,	9.18708,	10.8663,	12.7617,	14.8968,	17.2990,	20.0000};
  static const G4double IC[nX]={0,
	.021877,	.091263,	.214670,	.398677,	.650133,	.976322,	1.38510,	1.88504,	2.48555,	3.19709,
	4.03129,	5.00127,	6.12184,	7.40989,	8.88482,	10.5690,	12.4888,	14.6748,	17.1638,	20.0000};
  static const G4double ID[nX]={0,
	.020062,	.084127,	.198702,	.370384,	.606100,	.913288,	1.30006,	1.77535,	2.34912,	3.03253,
	3.83822,	4.78063,	5.87634,	7.14459,	8.60791,	10.2929,	12.2315,	14.4621,	17.0320,	20.0000};
  static const G4double IE[nX]={0,
	.018547,	.078104,	.185102,	.346090,	.567998,	.858331,	1.22535,	1.67824,	2.22735,	2.88443,
	3.66294,	4.57845,	5.64911,	6.89637,	8.34578,	10.0282,	11.9812,	14.2519,	16.8993,	20.0000};
  static const G4double IF[nX]={0,
	.017143,	.072466,	.172271,	.323007,	.531545,	.805393,	1.15288,	1.58338,	2.10754,	2.73758,
	3.48769,	4.37450,	5.41770,	6.64092,	8.07288,	9.74894,	11.7135,	14.0232,	16.7522,	20.0000};
  static const G4double IG[nX]={0,
	.015618,	.066285,	.158094,	.297316,	.490692,	.745653,	1.07053,	1.47479,	1.96931,	2.56677,
	3.28205,	4.13289,	5.14068,	6.33158,	7.73808,	9.40133,	11.3745,	13.7279,	16.5577,	20.0000};
  static const G4double IH[nX]={0,
	.013702,	.058434,	.139923,	.264115,	.437466,	.667179,	.961433,	1.32965,	1.78283,	2.33399,
	2.99871,	3.79596,	4.74916,	5.88771,	7.24937,	8.88367,	10.8576,	13.2646,	16.2417,	20.0000};
  static const G4double II[nX]={0,
	.011264,	.048311,	.116235,	.220381,	.366634,	.561656,	.813132,	1.13008,	1.52322,	2.00554,
	2.59296,	3.30542,	4.16834,	5.21490,	6.48964,	8.05434,	9.99835,	12.4580,	15.6567,	20.0000};
  static const G4double IJ[nX]={0,
	.008628,	.037206,	.089928,	.171242,	.286114,	.440251,	.640343,	.894382,	1.21208,	1.60544,
	2.08962,	2.68414,	3.41486,	4.31700,	5.44048,	6.85936,	8.69067,	11.1358,	14.5885,	20.0000};
  static const G4double* Il[nE]=
                             {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
  static const G4double lE[nE]={
-1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215,	.459141,	.867068,	1.27499,	1.68292,
 2.09085,	2.49877,	2.90670,	3.31463,	3.72255,	4.13048,	4.53840,	4.94633,	5.35426,	5.76218};
  static const G4double lEmi=lE[0];
  static const G4double lEma=lE[nE-1];
  static const G4double dlE=(lEma-lEmi)/bE;
	 //***************************************************************************************
  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
  G4double lEn=std::log(Enu);         // log(E) for interpolation
  G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
  G4int fE=static_cast<int>(rE);      // Left bin for interpolation
  if(fE<0) fE=0;
  if(fE>pE)fE=pE;
  G4int    sE=fE+1;                   // Right bin for interpolation
  G4double dE=rE-fE;                  // relative log shift from the left bin
  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
  G4double Emu=Enu-mmu;               // Free Energy of neutrino/anti-neutrino
  G4double ME=Enu*MN;                 // M*E
  G4double dME=ME+ME;                 // 2*M*E
  G4double dEMN=(dEnu+MN)*ME;
  G4double MEm=ME-hmmu2;
  G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
  G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
  G4double ymax=(E2M+sqE)/dEMN;
  G4double ymin=(E2M-sqE)/dEMN;
  G4double rmin=1.-ymin;
  G4double rhm2E=hmmu2/Enu2;
  G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
  G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
  G4double Q2nq=Emu*dMN-mcV;
  if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
  // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
  G4double Rmi=Q2mi/Q2ma;
  G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78);
  // --- E-interpolation must be done in a log scale ---
  G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
  G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
  G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
  G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
  G4double rXi=(Xmi-iXmi)/idX;
  G4int    iXi=static_cast<int>(rXi);
  if(iXi<0) iXi=0;
  if(iXi>bX) iXi=bX;
  G4double dXi=rXi-iXi;
  G4double bntil=Il[fE][iXi];
  G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
  G4double bntir=Il[sE][iXi];
  G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
  G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
  //
  G4double rXa=(Xma-iXmi)/idX;
  G4int    iXa=static_cast<int>(rXa);
  if(iXa<0) iXa=0;
  if(iXa>bX) iXa=bX;
  G4double dXa=rXa-iXa;
  G4double bntal=Il[fE][iXa];
  G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
  G4double bntar=Il[sE][iXa];
  G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
  G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
  //
  // *** Find X using the reversed table ***
  G4double intx=inti+(inta-inti)*G4UniformRand(); 
  G4int    intc=static_cast<int>(intx);
  if(intc<0) intc=0;
  if(intc>bX) intc=bX;
  G4double dint=intx-intc;
  G4double mXl=Xl[fE][intc];
  G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
  G4double mXr=Xl[sE][intc];
  G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
  G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
  G4double R=shift-std::pow(X,pconv);
  G4double Q2=R*Q2ma;
  return Q2*GeV*GeV;
}

// It returns a fraction of the direct interaction of the neutrino with quark-partons
G4double G4QANuMuNuclearCrossSection::GetDirectPart(G4double Q2)
{
  G4double f=Q2/4.62;
  G4double ff=f*f;
  G4double r=ff*ff;
  G4double s=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
  //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
  return 1.-s*(1.-s/2);
}

// #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
G4double G4QANuMuNuclearCrossSection::GetNPartons(G4double Q2)
{
  return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
}

// This class can provide only virtual exchange pi+ (a substitute for W+ boson)
G4int G4QANuMuNuclearCrossSection::GetExchangePDGCode() {return -211;}
