//
// ********************************************************************
// * 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: G4QNuENuclearCrossSection.cc,v 1.4 2009/05/08 15:16:26 mkossov Exp $
// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
//
//
// G4 Physics class: G4QNuENuclearCrossSection for A(nu_e,e-) cross sections
// Created: M.V. Kossov, CERN/ITEP(Moscow), 24-SEP-2007
// The last update: M.V. Kossov, CERN/ITEP (Moscow) 24-SEP-2007
// 
// ****************************************************************************************
// ********** 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 "G4QNuENuclearCrossSection.hh"

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

// Returns Pointer to the G4VQCrossSection class
G4VQCrossSection* G4QNuENuclearCrossSection::GetPointer()
{
  static G4QNuENuclearCrossSection 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 G4QNuENuclearCrossSection::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<<"G4QNENCS::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!=12)
  {
#ifdef debug
    G4cout<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pEn;
        }
        break;                           // Go out of the LOOP
      }
#ifdef pdebug
      G4cout<<"---G4QNENCrossSec::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<<"G4QNENCS::GetCrSec: CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<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<<"G4QNENCrossSection::GetCrossSect:NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
#endif
        if(pEn>lastTH)
        {
#ifdef pdebug
          G4cout<<"G4QNENCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pEn;
        }
      }
#ifdef pdebug
      G4cout<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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<<"G4QNENCS::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 G4QNuENuclearCrossSection::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 me=.00051099892; // electron mass in GeV
  static const G4double me2=me*me;    // Squared mass of an electron in GeV^2
  static const G4double thresh=me+me2/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=me+me2/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 G4QNuENuclearCrossSection::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()33/65!!
  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 me=.00051099892;// electron mass in GeV
  static const G4double me2=me*me;      // Squared mass of an electron in GeV^2
  static const G4double EMi=me+me2/dmN; // Universal threshold of the reaction in GeV
  static const G4double EMa=300.;       // Maximum tabulated Energy of nu_e 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 Electron
  //G4double TotEnergy2=Momentum;
  onlyCS=CS;                            // Flag to calculate only CS (not TX & QE)
  lastE=Momentum/GeV;                   // Kinetic energy of the electron 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*G4NuENuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
      // *** The synchronization check ***
      G4int sync=TX.size();
      if(sync!=I) G4cerr<<"***G4NuENuclearCS::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 the neutrinoEnergy 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*G4NuENuclearCS::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*G4NuENuclearCS::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<<"G4NuENuclearCS::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 G4QNuENuclearCrossSection::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 me=.00051099892; // electron mass in GeV
  static const G4double me2=me*me;    // Squared mass of an electron in GeV^2
  static const G4double thresh=me+me2/dmN; // Universal threshold in GeV
  static const G4int nE=65; // !! If change this, change it in CalculateCrossSection() !!
  static const G4double nuEn[nE]={thresh,
    .00051331,.00053602,.00056078,.00058783,.00061743,.00064990,.00068559,.00072492,
    .00076834,.00081641,.00086975,.00092912,.00099536,.00106950,.00115273,.00124646,
    .00135235,.00147241,.00160901,.00176503,.00194392,.00214986,.00238797,.00266448,
    .00298709,.00336531,.00381094,.00433879,.00496745,.00572047,.00662785,.00772806,
    .00907075,.01072050,.01276190,.01530660,.01850330,.02255110,.02771990,.03437780,
    .04303240,.05438970,.06944210,.08959920,.11688400,.15423600,.20597200,.27851200,
    .38153100,.52979600,.74616300,1.0665200,1.5480900,2.2834800,3.4251100,5.2281000,
    8.1270200,12.875900,20.808500,34.331200,57.877800,99.796200,176.16300,318.68200};
  static const G4double TOTX[nE]={0.,
    .00047551,.00162896,.00232785,.00292938,.00349456,.00404939,.00460908,.00518455,
    .00578488,.00641848,.00709376,.00781964,.00860585,.00946334,.01040460,.01144420,
    .01259910,.01388920,.01533830,.01697480,.01883260,.02095280,.02338500,.02618960,
    .02944060,.03322870,.03766580,.04289050,.04907540,.06123530,.07521120,.09034730,
    .10803800,.12856800,.15277600,.18174900,.21764300,.26083100,.31424900,.37935600,
    .45871900,.55375100,.66506600,.79039400,.92276600,1.0489000,1.1500300,1.2071700,
    1.2096800,1.1612200,1.0782900,.98251100,.89137400,.81472700,.75500100,.71061200,
    .67873900,.65619000,.64098400,.63085100,.62389900,.61664900,.61261000,.60635700};
  static const G4double QELX[nE]={0.,
  2.44084e-7,8.73147e-7,1.30540e-6,1.72196e-6,2.15765e-6,2.63171e-6,3.15996e-6,3.75836e-6,
  4.44474e-6,5.24008e-6,6.16982e-6,7.26537e-6,8.56595e-6,1.01211e-5,1.19937e-5,1.42647e-5,
  1.70384e-5,2.04506e-5,2.46796e-5,2.99611e-5,3.66090e-5,4.50456e-5,5.58425e-5,6.97817e-5,
  8.79417e-5,.000111825,.000143542,.000186093,.000243780,.000350295,.000498489,.000698209,
  .000979989,.001378320,.001949710,.002781960,.004027110,.005882030,.008710940,.013041400,
  .019739800,.030118300,.046183600,.070818700,.107857000,.161778000,.236873000,.336212000,
  .455841000,.565128000,.647837000,.701208000,.729735000,.742062000,.746495000,.748182000,
  .749481000,.750637000,.751471000,.752237000,.752763000,.753103000,.753159000,.753315000};

  // --------------------------------
  G4int first=0;
  if(z<0.)
  {
    first=1;
    z=-z;
  }
  if(z<1 || z>92)             // neutron and plutonium are forbidden
  {
    G4cout<<"***G4QNuENuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
    return -1;
  }
  for(G4int k=0; k<nE; k++)
  {
    G4double a=n+z;
    G4double na=n+a;
    G4double dn=n+n;
    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]*(na+na)/ta+QELX[k]*(dn+dn-da)/ta; // TotalCrossSection
    q[k]=QELX[k]*dn/a;                                     // QuasiElasticCrossSection
  }
  return first;
}

// Randomize Q2 from neutrino to the scattered electron when scattering is quasi-elastic
G4double G4QNuENuclearCrossSection::GetQEL_ExchangeQ2()
{
  static const G4double me=.00051099892; // electron mass in GeV
  static const G4double me2=me*me;     // Squared mass of an electron in GeV^2
  static const G4double hme2=me2/2;    // .5*m_e^2 in GeV^2
  static const G4double 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]={1.87905e-10,
 .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374,
 .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069,
 .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718,
 .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169,
 .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670,
 .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720,
 .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457,
 .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366,
 .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128,
 .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539};
  // 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.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811,
 15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039,
 26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876,
 36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950,
 46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377,
 55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162,
 65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852,
 74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767,
 83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106,
 92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 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-hme2;
  G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
  G4double E2M=MN*Enu2-(Enu+MN)*hme2;
  G4double ymax=(E2M+sqE)/dEMN;
  G4double ymin=(E2M-sqE)/dEMN;
  G4double rmin=1.-ymin;
  G4double rhm2E=hme2/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 electron when scattering is not quasiElastic
G4double G4QNuENuclearCrossSection::GetNQE_ExchangeQ2()
{
  static const double mpi=.13957018;    // charged pi meson mass in GeV
  static const G4double me=.00051099892;// electron mass in GeV
  static const G4double me2=me*me;      // Squared mass of an electron in GeV^2
  static const G4double hme2=me2/2;     // .5*m_e^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]={6.14081e-05,
 .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244,
 2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419};
  static const G4double X1[nX]={.00125268,
 .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049,
 4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330};
  static const G4double X2[nX]={.015694,
 1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728,
 10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276};
  static const G4double X3[nX]={.0866877,
 4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302,
 21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115};
  static const G4double X4[nX]={.160483,
 5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127,
 29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265};
  static const G4double X5[nX]={.0999307,
 5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398,
 26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649};
  static const G4double X6[nX]={.0276367,
 3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519,
 16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148};
  static const G4double X7[nX]={.00472383,
 2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473,
 9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057};
  static const G4double X8[nX]={.000630783,
 1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480,
 5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037};
  static const G4double X9[nX]={7.49179e-05,
 .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823,
 3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028};
  static const G4double XA[nX]={8.43437e-06,
 .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704,
 2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823};
  static const G4double XB[nX]={9.27028e-07,
 .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059,
 1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341};
  static const G4double XC[nX]={1.00807e-07,
 .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668,
 1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014};
  static const G4double XD[nX]={1.09102e-08,
 .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878,
 1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617};
  static const G4double XE[nX]={1.17831e-09,
 .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374,
 .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571};
  static const G4double XF[nX]={1.27141e-10,
 .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046,
 .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338};
  static const G4double XG[nX]={1.3713e-11,
 .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870,
 .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379};
  static const G4double XH[nX]={1.47877e-12,
 .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421,
 .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197};
  static const G4double XI[nX]={1.59454e-13,
 .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806,
 .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872};
  static const G4double XJ[nX]={1.71931e-14,
 .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147,
 .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813};
  // 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,
 .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555,
 13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000};
  static const G4double I1[nX]={0,
 .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365,
 13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000};
  static const G4double I2[nX]={0,
 .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586,
 12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000};
  static const G4double I3[nX]={0,
 .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461,
 12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000};
  static const G4double I4[nX]={0,
 .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515,
 11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000};
  static const G4double I5[nX]={0,
 .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040,
 10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000};
  static const G4double I6[nX]={0,
 .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243,
 9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000};
  static const G4double I7[nX]={0,
 .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717,
 8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000};
  static const G4double I8[nX]={0,
 .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415,
 8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000};
  static const G4double I9[nX]={0,
 .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672,
 8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000};
  static const G4double IA[nX]={0,
 .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265,
 7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000};
  static const G4double IB[nX]={0,
 .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298,
 7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000};
  static const G4double IC[nX]={0,
 .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696,
 7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000};
  static const G4double ID[nX]={0,
 .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212,
 7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000};
  static const G4double IE[nX]={0,
 .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949,
 7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000};
  static const G4double IF[nX]={0,
 .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164,
 7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000};
  static const G4double IG[nX]={0,
 .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249,
 7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000};
  static const G4double IH[nX]={0,
 .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804,
 7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000};
  static const G4double II[nX]={0,
 .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722,
 7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000};
  static const G4double IJ[nX]={0,
 .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196,
 6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 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 Ee=Enu-me;               // 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-hme2;
  G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
  G4double E2M=MN*Enu2-(Enu+MN)*hme2;
  G4double ymax=(E2M+sqE)/dEMN;
  G4double ymin=(E2M-sqE)/dEMN;
  G4double rmin=1.-ymin;
  G4double rhm2E=hme2/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=Ee*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=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu
  // --- 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 G4QNuENuclearCrossSection::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 G4QNuENuclearCrossSection::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 G4QNuENuclearCrossSection::GetExchangePDGCode() {return 211;}