G4QMuonNuclearCrossSection Class Reference

#include <G4QMuonNuclearCrossSection.hh>

Inheritance diagram for G4QMuonNuclearCrossSection:

Public Member Functions
	~G4QMuonNuclearCrossSection ()
G4double	ThresholdEnergy (G4int Z, G4int N, G4int PDG=13)
virtual G4double	GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=0)
G4double	CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum)
G4int	GetExchangePDGCode ()
G4double	GetExchangeEnergy ()
G4double	GetVirtualFactor (G4double nu, G4double Q2)
G4double	GetExchangeQ2 (G4double nu)
Public Member Functions inherited from G4VQCrossSection
virtual	~G4VQCrossSection ()
virtual G4double	GetCrossSection (G4bool, G4double, G4int, G4int, G4int pPDG=0)
virtual G4double	ThresholdEnergy (G4int Z, G4int N, G4int PDG=0)
virtual G4double	CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int PDG, G4int tgZ, G4int tgN, G4double pMom)=0
virtual G4double	GetLastTOTCS ()
virtual G4double	GetLastQELCS ()
virtual G4double	GetDirectPart (G4double Q2)
virtual G4double	GetNPartons (G4double Q2)
virtual G4double	GetExchangeEnergy ()
virtual G4double	GetExchangeT (G4int tZ, G4int tN, G4int pPDG)
virtual G4double	GetSlope (G4int tZ, G4int tN, G4int pPDG)
virtual G4double	GetHMaxT ()
virtual G4double	GetExchangeQ2 (G4double nu=0)
virtual G4double	GetVirtualFactor (G4double nu, G4double Q2)
virtual G4double	GetQEL_ExchangeQ2 ()
virtual G4double	GetNQE_ExchangeQ2 ()
virtual G4int	GetExchangePDGCode ()

Static Public Member Functions
static G4VQCrossSection *	GetPointer ()
Static Public Member Functions inherited from G4VQCrossSection
static void	setTolerance (G4double tol)

Protected Member Functions
	G4QMuonNuclearCrossSection ()
Protected Member Functions inherited from G4VQCrossSection
	G4VQCrossSection ()
G4double	LinearFit (G4double X, G4int N, G4double *XN, G4double *YN)
G4double	EquLinearFit (G4double X, G4int N, G4double X0, G4double DX, G4double *Y)

Additional Inherited Members
Static Protected Attributes inherited from G4VQCrossSection
static G4double	tolerance =.001

Detailed Description

Definition at line 56 of file G4QMuonNuclearCrossSection.hh.

Constructor & Destructor Documentation

G4QMuonNuclearCrossSection()

G4QMuonNuclearCrossSection::G4QMuonNuclearCrossSection ( )

inlineprotected

Definition at line 60 of file G4QMuonNuclearCrossSection.hh.

{};

~G4QMuonNuclearCrossSection()

G4QMuonNuclearCrossSection::~G4QMuonNuclearCrossSection

(

)

Definition at line 83 of file G4QMuonNuclearCrossSection.cc.

{
  G4int lens=J1->size();
  for(G4int i=0; i<lens; ++i) delete[] (*J1)[i];
  delete J1;
  G4int hens=J2->size();
  for(G4int i=0; i<hens; ++i) delete[] (*J2)[i];
  delete J2;
  G4int pens=J3->size();
  for(G4int i=0; i<pens; ++i) delete[] (*J3)[i];
  delete J3;
}

Member Function Documentation

CalculateCrossSection()

G4double G4QMuonNuclearCrossSection::CalculateCrossSection ( G4bool  CS, G4int  F, G4int  I, G4int  PDG, G4int  Z, G4int  N, G4double  Momentum  )

virtual

Implements G4VQCrossSection.

Definition at line 322 of file G4QMuonNuclearCrossSection.cc.

{
  static const G4int nE=336;    // !! If change this, change it in GetFunctions() (*.hh) !!
  static const G4int mL=nE-1;
  static const G4double EMi=2.0612;          // Minimum tabulated KineticEnergy of Muon
  static const G4double EMa=50000.;          // Maximum tabulated Energy of the Muon 
  static const G4double lEMi=std::log(EMi);  // Minimum tabulatedLogarithmKinEnergy of Muon
  static const G4double lEMa=std::log(EMa);  // Maximum tabulatedLogarithmKinEnergy of Muon
  static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic table step for MuonKinEnergy
  static const G4double alop=1./137.036/3.14159265; //coeffitient for  E>50000 calculations
  static const G4double mmu=105.65839;       // Mass of the muon in MeV
  static const G4double mmu2=mmu*mmu;        // Squared Mass of muon in MeV^2
  static const G4double lmmu=std::log(mmu);       // Log of the muon mass
  // *** Begin of the Associative memory for acceleration of the cross section calculations
  static std::vector <G4int> colF;     // Vector of LastStartPosition in Ji-function tables
  static std::vector <G4double> colH;  // Vector of HighEnergyCoefficients (functional)
#ifdef pdebug
  G4cout<<"G4QMuonNucCrossSection::CalculateCrossSection: ***Called*** "<<J3->size();
  if(J3->size()) G4cout<<", p="<<(*J3)[0];
  G4cout<<G4endl;
#endif
  // *** End of Static Definitions (Associative Memory) ***
  //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
  onlyCS=CS;                           // Flag to calculate only CS (not Si/Bi)
  G4double TotEnergy2=Momentum*Momentum+mmu2;
  G4double TotEnergy=std::sqrt(TotEnergy2); // Total energy of the muon
  lastE=TotEnergy-mmu;               // Kinetic energy of the muon
#ifdef pdebug
  G4cout<<"G4QMuonNucCS::CalcCS: P="<<Momentum<<", F="<<F<<", I="<<I<<", Z="<<targZ
        <<", N="<<targN<<", onlyCS="<<CS<<",E="<<lastE<<",th="<<EMi<<G4endl;
#endif
  G4double A=targN+targZ;            // New A (can be different from targetAtomicNumber)
  if(F<=0)                           // >>>This isotope was not the last used isotop<<<
  {
    if(F<0)                          // This isotope was found in DAMDB =------=> RETRIEVE
    {                                // ...........................................=------=
      if (lastE<=EMi)                // Energy is below the minimum energy in the table
      {
        lastE=0.;
        lastG=0.;
        lastSig=0.;
#ifdef pdebug
        G4cout<<"--> G4QMuonNucCS::CalcCS: Old CS=0  as lastE="<<lastE<<" < "<<EMi<<G4endl;
#endif
        return 0.;
      }
      lastJ1 =(*J1)[I];              // Pointer to the prepared J1 function
      lastJ2 =(*J2)[I];              // Pointer to the prepared J2 function
      lastJ3 =(*J3)[I];              // Pointer to the prepared J3 function
      lastF  =colF[I];               // Last ZeroPosition in the J-functions
      lastH  =colH[I];               // Last High Energy Coefficient (A-dependent)
    }
    else                             // This isotope wasn't calculated previously => CREATE
    {                                // ............................................=-----=
      lastJ1 = new G4double[nE];     // Allocate memory for the new J1 function
      lastJ2 = new G4double[nE];     // Allocate memory for the new J2 function
      lastJ3 = new G4double[nE];     // Allocate memory for the new J3 function
      lastF   = GetFunctions(A,lastJ1,lastJ2,lastJ3);//newZeroPos and J-functions filling
      lastH   = alop*A*(1.-.072*std::log(A)); // like lastSP of G4PhotonuclearCrossSection
#ifdef pdebug
      G4cout<<"==>G4QMuNCS::CalcCS:lJ1="<<lastJ1<<",lJ2="<<lastJ2<<",lJ3="<<lastJ3<<G4endl;
#endif
      // *** The synchronization check ***
      G4int sync=J1->size();
      if(sync!=I) G4cerr<<"***G4QMuonNuclearCS::CalcCrS: PDG=13, S="<<sync<<"#"<<I<<G4endl;
      J1->push_back(lastJ1);
      J2->push_back(lastJ2);
      J3->push_back(lastJ3);
      colF.push_back(lastF);
      colH.push_back(lastH);
    } // End of creation of the new set of parameters
  } // End of parameters udate
  // =-------------------= NOW Calculate the Cross Section =--------------------------=
  if (lastE<=lastTH || lastE<=EMi)   // Check that muKiE is higher than ThreshE
  {
    lastE=0.;
    lastG=0.;
    lastSig=0.;
#ifdef pdebug
    G4cout<<"--> G4QMuonNucCS::CalcCS:CS=0 as T="<<lastE<<"<"<<EMi<<" || "<<lastTH<<G4endl;
#endif
    return 0.;
  }
  G4double lE=std::log(lastE);       // log(muE) (it is necessary for the fit)
  lastG=lE-lmmu;                     // Gamma of the muon (used to recover log(muE))
  G4double dlg1=lastG+lastG-1.;
  G4double lgoe=lastG/lastE;
  if(lE<lEMa) // Log fit is made explicitly to fix the last bin for the randomization
  {
    G4double shift=(lE-lEMi)/dlnE;
    G4int    blast=static_cast<int>(shift);
#ifdef pdebug
    G4cout<<"-->G4QMuonNuclearCS::CalcCrossSect:LOGfit b="<<blast<<",max="<<mL<<",lJ1="
          <<lastJ1<<",lJ2="<<lastJ2<<",lJ3="<<lastJ3<<G4endl;
#endif
    if(blast<0)   blast=0;
    if(blast>=mL) blast=mL-1;
    shift-=blast;
    lastL=blast+1;
    G4double YNi=dlg1*lastJ1[blast]-lgoe*(lastJ2[blast]+lastJ2[blast]-lastJ3[blast]/lastE);
    G4double YNj=dlg1*lastJ1[lastL]-lgoe*(lastJ2[lastL]+lastJ2[lastL]-lastJ3[lastL]/lastE);
    lastSig = YNi+shift*(YNj-YNi);
    if(lastSig>YNj)lastSig=YNj;
#ifdef pdebug
    G4cout<<"G4QMuNucCS::CalcCS:S="<<lastSig<<",E="<<lE<<",Yi="<<YNi<<",Yj="<<YNj<<",M="
          <<lEMa<<G4endl;
    G4cout<<"G4QMuNucCS::CalcCS:s="<<shift<<",Jb="<<lastJ1[blast]<<",J="<<lastJ1[lastL]
          <<",b="<<blast<<G4endl;
#endif
  }
  else
  {
#ifdef pdebug
    G4cout<<"->G4QMuonNucCS::CCS:LOGex="<<lastJ1<<",lJ2="<<lastJ2<<",lJ3="<<lastJ3<<G4endl;
#endif
    lastL=mL;
    G4double term1=lastJ1[mL]+lastH*HighEnergyJ1(lE);
    G4double term2=lastJ2[mL]+lastH*HighEnergyJ2(lE);
    G4double term3=lastJ3[mL]+lastH*HighEnergyJ3(lE);
    lastSig=dlg1*term1-lgoe*(term2+term2-term3/lastE);
#ifdef pdebug
    G4cout<<"G4QMuNucCS::CalculateCrossSection:S="<<lastSig<<",lE="<<lE<<",J1="
          <<lastH*HighEnergyJ1(lE)<<",Pm="<<lastJ1[mL]<<",Fm="<<lastJ2[mL]<<",Fh="
          <<lastH*HighEnergyJ2(lE)<<",EM="<<lEMa<<G4endl;
#endif
  }
  if(lastSig<0.) lastSig = 0.;
  return lastSig;
}

Referenced by GetCrossSection().

GetCrossSection()

G4double G4QMuonNuclearCrossSection::GetCrossSection ( G4bool  fCS, G4double  pMom, G4int  tgZ, G4int  tgN, G4int  pPDG = 0  )

virtual

!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)

Reimplemented from G4VQCrossSection.

Definition at line 98 of file G4QMuonNuclearCrossSection.cc.

{
  static const G4double mmu=105.65839; // Mass of the muon in MeV
  static const G4double mmu2=mmu*mmu;  // Squared Mass of muon in MeV^2
  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=std::sqrt(pMom*pMom+mmu2)-mmu; // ==> mu-/mu+ kinEnergy
#ifdef pdebug
  G4cout<<"G4QMNCS::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(std::abs(pPDG)!=13)
  {
#ifdef pdebug
    G4cout<<"G4QMNCS::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<<"G4QMNCS::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<<"G4QMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=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)
        if(lastP==pMom)                  // VI do not use tolerance
        {
#ifdef pdebug
          G4cout<<"G4QMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
#endif
          CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only
          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<<"G4QMNCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
#endif
        lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
#ifdef pdebug
        G4cout<<"G4QMNCS::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<<"G4QMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pEn;
        }
        break;                           // Go out of the LOOP
      }
#ifdef pdebug
      G4cout<<"---G4QMNCrossSec::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<<"G4QMNCS::GetCrosSec: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<<"G4QMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
#endif
        if(pEn>lastTH)
        {
#ifdef pdebug
          G4cout<<"G4QMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pEn;
        }
      }
#ifdef pdebug
      G4cout<<"G4QMNCS::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<<"G4QMNCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<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<<"G4QMNCS::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<<"G4QMNCS::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)
  else if(lastP==pMom)                 // VI do not use tolerance
  {
#ifdef pdebug
    G4cout<<"G4QMNCS::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<<"G4QMNCS::GetCS:UpdatCur 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<<"G4QMNCS::GetCroSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
  //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
  return lastCS*millibarn;
}

GetExchangeEnergy()

G4double G4QMuonNuclearCrossSection::GetExchangeEnergy ( )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 2618 of file G4QMuonNuclearCrossSection.cc.

{
  // @@ All constants are copy of that from GetCrossSection funct. => Make them general.
  static const G4int nE=336; // !!  If change this, change it in GetFunctions() (*.hh) !!
  static const G4int mL=nE-1;
  static const G4double EMi=2.0612;          // Minimum Energy
  static const G4double EMa=50000.;          // Maximum Energy
  static const G4double lEMi=std::log(EMi);       // Minimum logarithmic Energy
  static const G4double lEMa=std::log(EMa);       // Maximum logarithmic Energy
  static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic step in Energy
  static const G4double mmu=105.65839;       // Mass of muon in MeV
  static const G4double lmmu=std::log(mmu);       // Log of muon mass
  G4double phLE=0.;                     // Prototype of the log(nu=E_gamma)
  G4double Y[nE];                       // Prepare the array for randomization
#ifdef debug
  G4cout<<"G4QMuonNuclCrossSect::GetExchanEn:B="<<lastF<<",l="<<lastL<<",1="<<lastJ1[lastL]
        <<",2="<<lastJ2[lastL]<<",3="<<lastJ3[lastL]<<",S="<<lastSig<<",E="<<lastE<<G4endl;
#endif
  G4double lastLE=lastG+lmmu;           // recover log(eE) from the gamma (lastG)
  G4double dlg1=lastG+lastG-1.;
  G4double lgoe=lastG/lastE;
  for(G4int i=lastF;i<=lastL;i++)
    Y[i]=dlg1*lastJ1[i]-lgoe*(lastJ2[i]+lastJ2[i]-lastJ3[i]/lastE);
  G4double ris=lastSig*G4UniformRand(); // Sig can be > Y[lastL=mL], then it is funct. reg.
#ifdef debug
  G4cout<<"G4QMuonNuclearCrossSection::GetExchangeEnergy: "<<ris<<",Y="<<Y[lastL]<<G4endl;
#endif
  if(ris<Y[lastL])                      // Search in the table
  {
    G4int j=lastF;
    G4double Yj=Y[j];                   // It mast be 0 (some times just very small)
    while (ris>Yj && j<lastL)           // Associative search
    {
      j++;
      Yj=Y[j];                          // High value
    }
    G4int j1=j-1;
    G4double Yi=Y[j1];                  // Low value
    phLE=lEMi+(j1+(ris-Yi)/(Yj-Yi))*dlnE;
#ifdef debug
    G4cout<<"G4QMuNuclearCS::E="<<phLE<<",l="<<lEMi<<",j="<<j<<",ris="<<ris<<",Yi="<<Yi
          <<",Y="<<Yj<<G4endl;
#endif
  }
  else                                  // Search with the function
  {
    if(lastL<mL)
      G4cerr<<"**G4QMuonNucCS::GetExEn:L="<<lastL<<",S="<<lastSig<<",Y="<<Y[lastL]<<G4endl;
    G4double f=(ris-Y[lastL])/lastH;    // ScaledResidualValue of the cross-sec. integral
#ifdef pdebug
    G4cout<<"G4QMuNucCS::GetExEn:HighEnergy f="<<f<<",ris="<<ris<<",lastH="<<lastH<<G4endl;
#endif
    phLE=SolveTheEquation(f);           // Solve equation to find theLog(phE) (comp lastLE)
#ifdef pdebug
    G4cout<<"G4QMuonNuclearCrossSection::GetExchangeEnergy:HighEnergy lphE="<<phLE<<G4endl;
#endif
  }
  if(phLE>lastLE)
  {
    G4cerr<<"***G4QMuNuclearCS::GetExEnergy: N="<<lastN<<",Z="<<lastZ<<",lpE"<<phLE<<">leE"
          <<lastLE<<",Sig="<<lastSig<<",rndSig="<<ris<<",Beg="<<lastF<<",End="<<lastL
          <<",Y="<<Y[lastL]<<G4endl;
    if(lastLE<7.2) phLE=std::log(std::exp(lastLE)-mmu);
    else phLE=7.;
  }
  return std::exp(phLE);
}

GetExchangePDGCode()

G4int G4QMuonNuclearCrossSection::GetExchangePDGCode ( )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 2787 of file G4QMuonNuclearCrossSection.cc.

{return 22;}

GetExchangeQ2()

G4double G4QMuonNuclearCrossSection::GetExchangeQ2 ( G4double  nu )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 2733 of file G4QMuonNuclearCrossSection.cc.

{
  static const G4double mmu=105.65839;  // Mass of muon in MeV
  static const G4double mmu2=mmu*mmu;   // Squared Mass of muon in MeV
  G4double y=nu/lastE;                  // Part of energy carried by the equivalent pfoton
  if(y>=1.-1./(lastG+lastG)) return 0.; // The region where the method does not work
  G4double y2=y*y;                      // Squared photonic part of energy
  G4double ye=1.-y;                     // Part of energy carried by the secondary muon
  G4double Qi2=mmu2*y2/ye;              // Minimum Q2
  G4double Qa2=4*lastE*lastE*ye;        // Maximum Q2
  G4double iar=Qi2/Qa2;                 // Q2min/Q2max ratio
  G4double Dy=ye+.5*y2;                 // D(y) function
  G4double Py=ye/Dy;                    // P(y) function
  G4double ePy=1.-std::exp(Py);              // 1-exp(P(y)) part
  G4double Uy=Py*(1.-iar);              // U(y) function
  G4double Fy=(ye+ye)*(1.+ye)*iar/y2;   // F(y) function
  G4double fr=iar/(1.-ePy*iar);         // Q-fraction
  if(Fy<=-fr)
  {
#ifdef edebug
    G4cerr<<"***G4QMuonNucCS::GetExchQ2: Fy="<<Fy<<"+fr="<<fr<<" <0"<<",iar="<<iar<<G4endl;
#endif
    return 0.;
  }    
  G4double LyQa2=std::log(Fy+fr);            // L(y,Q2max) function
  G4bool cond=true;
  G4int maxTry=3;
  G4int cntTry=0;
  G4double Q2=Qi2;
  while(cond&&cntTry<maxTry)            // The loop to avoid x>1.
  {
    G4double R=G4UniformRand();         // Random number (0,1)
    Q2=Qi2*(ePy+1./(std::exp(R*LyQa2-(1.-R)*Uy)-Fy));
    cntTry++;
    cond = Q2>1878.*nu;
  }
  if(Q2<Qi2)
  {
#ifdef edebug
    G4cerr<<"***G4QMuonNuclearCrossSect::GetExchangeQ2: Q2="<<Q2<<" < Q2min="<<Qi2<<G4endl;
#endif
    return Qi2;
  }  
  if(Q2>Qa2)
  {
#ifdef edebug
    G4cerr<<"***G4QMuonNuclearCrossSect::GetExchangeQ2: Q2="<<Q2<<" > Q2max="<<Qi2<<G4endl;
#endif
    return Qa2;
  }  
  return Q2;
}

GetPointer()

G4VQCrossSection * G4QMuonNuclearCrossSection::GetPointer ( )

static

Definition at line 77 of file G4QMuonNuclearCrossSection.cc.

{
  static G4QMuonNuclearCrossSection theCrossSection; //**Static body of the Cross Section**
  return &theCrossSection;
}

Referenced by G4QInelastic::GetMeanFreePath(), G4QAtomicElectronScattering::PostStepDoIt(), and G4QInelastic::PostStepDoIt().

GetVirtualFactor()

G4double G4QMuonNuclearCrossSection::GetVirtualFactor ( G4double  nu, G4double  Q2  )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 2789 of file G4QMuonNuclearCrossSection.cc.

{
  static const G4double dM=938.27+939.57;// Mean double nucleon mass = m_n+m_p (no binding)
  static const G4double Q0=843.;         // Coefficient of the dipole nucleonic form-factor
  static const G4double Q02=Q0*Q0;       // Squared coefficient of theDipoleNuclFormFactor
  static const G4double blK0=std::log(185.);  // Coefficient of the b-function
  static const G4double bp=0.85;          // Power of the b-function
  static const G4double clK0=std::log(1390.);  // Coefficient of the c-function
  static const G4double cp=3.;            // Power of the c-function
  //G4double x=Q2/dM/nu;                  // Direct x definition
  G4double K=nu-Q2/dM;                    // K=nu*(1-x)
  if(K<0.)
  {
#ifdef edebug
    G4cerr<<"**G4QMuonNucCS::GetVirtFact:K="<<K<<",nu="<<nu<<",Q2="<<Q2<<",d="<<dM<<G4endl;
#endif
    return 0.;
  }
  G4double lK=std::log(K);                // ln(K)
  G4double x=1.-K/nu;                     // This definitin saves one div.
  G4double GD=1.+Q2/Q02;                  // Reversed nucleonic form-factor
  G4double b=std::exp(bp*(lK-blK0));      // b-factor
  G4double c=std::exp(cp*(lK-clK0));      // c-factor
  G4double r=.5*std::log(Q2+nu*nu)-lK;    // r=.5*log((Q^2+nu^2)/K^2)
  G4double ef=std::exp(r*(b-c*r*r));      // exponential factor
  return (1.-x)*ef/GD/GD;
}

ThresholdEnergy()

G4double G4QMuonNuclearCrossSection::ThresholdEnergy ( G4int  Z, G4int  N, G4int  PDG = 13  )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 276 of file G4QMuonNuclearCrossSection.cc.

{
  // CHIPS - Direct GEANT
  //static const G4double mNeut = G4QPDGCode(2112).GetMass();
  //static const G4double mProt = G4QPDGCode(2212).GetMass();
  static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/MeV;
  static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/MeV;
  static const G4double mAlph = G4NucleiProperties::GetNuclearMass(4,2)/MeV;
  // ---------
  static const G4double infEn = 9.e27;
 
  G4int A=Z+N;
  if(A<1) return infEn;
  else if(A==1) return 160.; // Pi0 threshold in MeV for the proton:T>m+(m^2+2lm)/2M,l=m_mu
  // CHIPS - Direct GEANT
  //G4double mT= G4QPDGCode(111).GetNuclMass(Z,N,0);
  G4double mT= 0.;
  if(Z&&G4NucleiProperties::IsInStableTable(A,Z))
                            mT = G4NucleiProperties::GetNuclearMass(A,Z)/MeV;
  else return 0.;       // If it is not in the Table of Stable Nuclei, then the Threshold=0
  // --------- Splitting thresholds
  G4double mP= infEn;
  if(A>1&&Z&&G4NucleiProperties::IsInStableTable(A-1,Z-1))
          mP = G4NucleiProperties::GetNuclearMass(A-1,Z-1)/MeV;// ResNucMass for a proton
 
  G4double mN= infEn;
  if(A>1&&G4NucleiProperties::IsInStableTable(A-1,Z))
          mN = G4NucleiProperties::GetNuclearMass(A-1,Z)/MeV;  // ResNucMass for a neutron
 
  G4double mA= infEn;
  if(A>4&&Z>1&&G4NucleiProperties::IsInStableTable(A-4,Z-2))
          mA=G4NucleiProperties::GetNuclearMass(A-4,Z-2)/MeV;  // ResNucMass for an alpha
 
  G4double dP= mP +mProt - mT;
  G4double dN= mN +mNeut - mT;
  G4double dA= mA +mAlph - mT;
#ifdef pdebug
  G4cout<<"G4QMuonNucCS::ThreshEn: mP="<<mP<<",dP="<<dP<<",mN="<<mN<<",dN="<<dN<<",mA="
        <<mA<<",dA="<<dA<<",mT="<<mT<<",A="<<A<<",Z="<<Z<<G4endl;
#endif
  if(dP<dN)dN=dP;
  if(dA<dN)dN=dA;
  return dN;
}

Referenced by GetCrossSection().

---

The documentation for this class was generated from the following files:

• G4QMuonNuclearCrossSection.hh
• G4QMuonNuclearCrossSection.cc