G4QIonIonCrossSection Class Reference

#include <G4QIonIonCrossSection.hh>

Inheritance diagram for G4QIonIonCrossSection:

Public Member Functions
	~G4QIonIonCrossSection ()
virtual G4double	GetCrossSection (G4bool fCS, G4double pMom, G4int Z, G4int N, G4int PDG)
G4double	CalculateCrossSection (G4bool fCS, G4int F, G4int I, G4int PDG, G4int tZ, G4int tN, G4double Momentum)
G4double	ThresholdMomentum (G4int pZ, G4int pN, G4int tZ, G4int tN)
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
	G4QIonIonCrossSection ()
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 59 of file G4QIonIonCrossSection.hh.

Constructor & Destructor Documentation

G4QIonIonCrossSection()

G4QIonIonCrossSection::G4QIonIonCrossSection ( )

inlineprotected

Definition at line 63 of file G4QIonIonCrossSection.hh.

{}

~G4QIonIonCrossSection()

G4QIonIonCrossSection::~G4QIonIonCrossSection ( )

inline

Definition at line 67 of file G4QIonIonCrossSection.hh.

{}

Member Function Documentation

CalculateCrossSection()

G4double G4QIonIonCrossSection::CalculateCrossSection ( G4bool  fCS, G4int  F, G4int  I, G4int  PDG, G4int  tZ, G4int  tN, G4double  Momentum  )

virtual

Implements G4VQCrossSection.

Definition at line 252 of file G4QIonIonCrossSection.cc.

{
  //static const G4double third=1./3.; // power for A^P->R conversion [R=1.1*A^(1/3)]
  //static const G4double conv=38.; // coeff. R2->sig=c*(pR+tR)^2, c=pi*10(mb/fm^2)*1.21
  // If change the following, please change in ::GetFunctions:
  static const G4double THmin=0.;  // @@ start from threshold (?) minimum Energy Threshold
  static const G4double dP=10.;    // step for the LEN table
  static const G4int    nL=100;    // A#of LENesonance points in E (each MeV from 2 to 106)
  static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part
  static const G4double Pmax=300000.;   // maxE for the HighE part
  static const G4int    nH=100;         // A#of HResonance points in lnE
  static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part
  static const G4double malP=std::log(Pmax); // High logarithm energy (each 2.75 percent)
  static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part
  //
  // Associative memory for acceleration
  static std::vector <G4double*> LENI;   // Vector of pointers: LowEnIneIonIonCrossSection
  static std::vector <G4double*> HENI;   // Vector of pointers: HighEnIneIonIonCrossSection
  static std::vector <G4double*> LENE;   // Vector of pointers: LowEnElaIonIonCrossSection
  static std::vector <G4double*> HENE;   // Vector of pointers: HighEnElaIonIonCrossSection
#ifdef debug
  G4cout<<"G4QIonIonCrossSection::CalcCS: Z="<<tZ<<", N="<<tN<<", P="<<TotMom<<G4endl;
#endif
  G4int dPDG=pPDG/10;                // 10SZZZAAA
  G4int zPDG=dPDG/1000;              // 10SZZZ (?)
  G4int zA=dPDG%1000;                // proj A
  G4int pZ=zPDG%1000;                // proj Z (?)
  G4int pN=zA-pZ;                    // proj N (?)
  G4double Momentum=TotMom/zA;       // Momentum per nucleon
  if (Momentum<THmin) return 0.;     // @@ This can be dangerouse for the heaviest nuc.!
  G4double sigma=0.;
  if(F&&I) sigma=0.;                 // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
  G4double tA=tN+tZ;                 // Target weight
  G4double pA=zA;                    // Projectile weight
  if(F<=0)                           // This isotope was not the last used isotop
  {
    if(F<0 || !XS)                   // This isotope was found in DAMDB or Elast =>RETRIEVE
    {
      lastLENI=LENI[I];              // Pointer to Low Energy inelastic cross sections
      lastHENI=HENI[I];              // Pointer to High Energy inelastic cross sections
      lastLENE=LENE[I];              // Pointer to Low Energy inelastic cross sections
      lastHENE=HENE[I];              // Pointer to High Energy inelastic cross sections
    }
    else                             // This isotope wasn't calculated previously => CREATE
    {
      lastLENI = new G4double[nL];   // Allocate memory for the new LEN cross sections
      lastHENI = new G4double[nH];   // Allocate memory for the new HEN cross sections
      lastLENE = new G4double[nL];   // Allocate memory for the new LEN cross sections
      lastHENE = new G4double[nH];   // Allocate memory for the new HEN cross sections
      G4int er=GetFunctions(pZ,pN,tZ,tN,lastLENI,lastHENI,lastLENE,lastHENE);
      if(er<1) G4cerr<<"*W*G4QIonIonCroSec::CalcCrossSection: pA="<<tA<<",tA="<<tA<<G4endl;
#ifdef debug
      G4cout<<"G4QIonIonCrossSection::CalcCS: GetFunctions er="<<er<<",pA="<<pA<<",tA="<<tA
            <<G4endl;
#endif
      // *** The synchronization check ***
      G4int sync=LENI.size();
      if(sync!=I) G4cout<<"*W*G4IonIonCrossSec::CalcCrossSect:Sync="<<sync<<"#"<<I<<G4endl;
      LENI.push_back(lastLENI);      // added LEN Inelastic
      HENI.push_back(lastHENI);      // added HEN Inelastic
      LENE.push_back(lastLENE);      // added LEN Elastic
      HENE.push_back(lastHENE);      // added HEN Elastic
    } // End of creation of the new set of parameters
  } // End of parameters udate
  // =------------= NOW the Magic Formula =--------------------------=
  if (Momentum<lastTH) return 0.;    // It must be already checked in the interface class
  else if (Momentum<Pmin)            // LEN region (approximated in E, not in lnE)
  {
#ifdef debug
    G4cout<<"G4QIICS::CalCS:p="<<pA<<",t="<<tA<<",n="<<nL<<",T="<<THmin<<",d="<<dP<<G4endl;
#endif
    if(tA<1. || pA<1.)
    {
      G4cout<<"-Warning-G4QIICS::CalcCS: pA="<<pA<<" or tA="<<tA<<" aren't nuclei"<<G4endl;
      sigma=0.;
    }
    else
    {
      G4double dPp=dP*pA;
      if(XS) sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENI);
      else   sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENE);
    }
#ifdef debugn
    if(sigma<0.) G4cout<<"-Warning-G4QIICS::CalcCS:pA="<<pA<<",tA="<<tA<<",XS="<<XS<<",P="
                       <<Momentum<<", Th="<<THmin<<", dP="<<dP<<G4endl;
#endif
  }
  else if (Momentum<Pmax*pA)                     // High Energy region
  {
    G4double lP=std::log(Momentum);
#ifdef debug
    G4cout<<"G4QIonIonCS::CalcCS:before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl;
#endif
    if(tA<1. || pA<1.)
    {
      G4cout<<"-Warning-G4QIICS::CalCS:pA="<<pA<<" or tA="<<tA<<" aren't composit"<<G4endl;
      sigma=0.;
    }
    else
    {
      G4double milPp=milP+std::log(pA);
      if(XS) sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENI);
      else   sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENE);
    }
  }
  else                                      // UltraHighE region (not frequent)
  {
    std::pair<G4double, G4double> inelel = CalculateXS(pZ, pN, tZ, tN, Momentum);
    if(XS) sigma=inelel.first;
    else   sigma=inelel.second;
  }
#ifdef debug
  G4cout<<"G4IonIonCrossSection::CalculateCrossSection: sigma="<<sigma<<G4endl;
#endif
  if(sigma<0.) return 0.;
  return sigma;
}

Referenced by GetCrossSection().

GetCrossSection()

G4double G4QIonIonCrossSection::GetCrossSection ( G4bool  fCS, G4double  pMom, G4int  Z, G4int  N, G4int  PDG  )

virtual

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

Reimplemented from G4VQCrossSection.

Definition at line 73 of file G4QIonIonCrossSection.cc.

{
  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> colICS;// Vector of last inelastic cross-sections
  static std::vector <G4double> colECS;// Vector of last elastic cross-sections
  // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
#ifdef pdebug
  G4cout<<"G4QIICS::GetCS:>>> f="<<fCS<<", Z="<<tZ<<"("<<lastZ<<"), N="<<tN<<"("<<lastN
        <<"),PDG="<<pPDG<<"("<<lastPDG<<"), p="<<pMom<<"("<<lastTH<<")"<<",Sz="
        <<colN.size()<<G4endl;
#endif
  if(!pPDG)
  {
#ifdef pdebug
    G4cout<<"G4QIonIonCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
#endif
    return 0.;                         // projectile PDG=0 is a mistake (?!) @@
  }
  G4bool in=false;                     // By default the isotope must be found in the AMDB
  if(tN!=lastN || tZ!=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   = tN;                      // The last N of the calculated nucleus
    lastZ   = tZ;                      // 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
#ifdef pdebug
    G4cout<<"G4QIICS::GetCS:FindI="<<lastI<<",pPDG="<<pPDG<<",tN="<<tN<<",tZ="<<tZ<<G4endl;
#endif
    if(lastI) for(G4int i=0; i<lastI; i++) // Loop over all DB
    {                                  // The nucleus with projPDG is found in AMDB
#ifdef pdebug
      G4cout<<"G4QII::GCS:P="<<colPDG[i]<<",N="<<colN[i]<<",Z="<<colZ[i]<<",j="<<j<<G4endl;
#endif
      if(colPDG[i]==pPDG && colN[i]==tN && colZ[i]==tZ)
      {
        lastI=i;
        lastTH =colTH[i];                // Last THreshold (A-dependent)
#ifdef pdebug
        G4cout<<"G4QIICS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
#endif
        if(pMom<=lastTH)
        {
#ifdef pdebug
          G4cout<<"G4QIICS::GetCS:Found P="<<pMom<<"<Threshold="<<lastTH<<"->XS=0"<<G4endl;
#endif
          return 0.;                     // Energy is below the Threshold value
        }
        lastP  =colP [i];                // Last Momentum  (A-dependent)
        lastICS=colICS[i];               // Last Inelastic Cross-Section (A-dependent)
        lastECS=colECS[i];               // Last Elastic Cross-Section (A-dependent)
        if(std::fabs(lastP/pMom-1.)<tolerance)
        {
#ifdef pdebug
          G4cout<<"G4QIonIonCS::GetCS:P="<<pMom<<",InXS="<<lastICS*millibarn<<",ElXS="
                <<lastECS*millibarn<<G4endl;
#endif
          CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only
          if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
          return         lastECS*millibarn;     // Use theLastElasticCS
        }
        in = true;                       // This is the case when the isotop is found in DB
        // Momentum pMom is in IU ! @@ Units
#ifdef pdebug
        G4cout<<"G4QIICS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
#endif
        lastICS=CalculateCrossSection( true,-1,j,lastPDG,lastZ,lastN,pMom);// read & update
        lastECS=CalculateCrossSection(false,-1,j,lastPDG,lastZ,lastN,pMom);// read & update
#ifdef pdebug
        G4cout<<"G4QIonIonCS::GetCS:=>New(inDB) InCS="<<lastICS<<",ElCS="<<lastECS<<G4endl;
#endif
        if((lastICS<=0. || lastECS<=0.) && pMom>lastTH) // Correct the threshold
        {
#ifdef pdebug
          G4cout<<"G4QIonIonCS::GetCS:New,T="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pMom;
        }
        break;                           // Go out of the LOOP
      }
#ifdef pdebug
      G4cout<<"--->G4QIonIonCrossSec::GetCrosSec: pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
            <<",Z["<<i<<"]="<<colZ[i]<<",PDG="<<colPDG[i]<<G4endl;
#endif
      j++;                             // Increment a#0f records found in DB for this pPDG
    }
    if(!in)                            // This nucleus has not been calculated previously
    {
#ifdef pdebug
      G4cout<<"G4QIICS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
#endif
      //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
      lastICS=CalculateCrossSection(true ,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
      lastECS=CalculateCrossSection(false,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
      if(lastICS<=0. || lastECS<=0.)
      {
        lastTH = ThresholdEnergy(tZ, tN); // Threshold Energy=Mom=0 which is now the last
#ifdef pdebug
        G4cout<<"G4QIonIonCrossSect::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
#endif
        if(pMom>lastTH)
        {
#ifdef pdebug
          G4cout<<"G4QIonIonCS::GetCS:1-st,P="<<pMom<<">Thresh="<<lastTH<<"->XS=0"<<G4endl;
#endif
          lastTH=pMom;
        }
      }
#ifdef pdebug
      G4cout<<"G4QIICS::GetCS: *New* ICS="<<lastICS<<", ECS="<<lastECS<<",N="<<lastN<<",Z="
            <<lastZ<<G4endl;
#endif
      colN.push_back(tN);
      colZ.push_back(tZ);
      colPDG.push_back(pPDG);
      colP.push_back(pMom);
      colTH.push_back(lastTH);
      colICS.push_back(lastICS);
      colECS.push_back(lastECS);
#ifdef pdebug
      G4cout<<"G4QIICS::GetCS:*1st*, P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn
            <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
#endif
      if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
      return         lastECS*millibarn;     // Use theLastElasticCS
    } // End of creation of the new set of parameters
    else
    {
#ifdef pdebug
      G4cout<<"G4QIICS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
#endif
      colP[lastI]=pMom;
      colPDG[lastI]=pPDG;
      colICS[lastI]=lastICS;
      colECS[lastI]=lastECS;
    }
  } // End of parameters udate
  else if(pMom<=lastTH)
  {
#ifdef pdebug
    G4cout<<"G4QIICS::GetCS: Current T="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
#endif
    return 0.;                         // Momentum is below the Threshold Value -> CS=0
  }
  else if(std::fabs(lastP/pMom-1.)<tolerance)
  {
#ifdef pdebug
    G4cout<<"G4QIICS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", InCS="<<lastICS*millibarn
          <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
#endif
    if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
    return         lastECS*millibarn;     // Use theLastElasticCS
  }
  else
  {
#ifdef pdebug
    G4cout<<"G4QIICS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
#endif
    lastICS=CalculateCrossSection( true,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
    lastECS=CalculateCrossSection(false,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
    lastP=pMom;
  }
#ifdef pdebug
  G4cout<<"G4QIICS::GetCroSec:*End*,P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn<<", ElCS="
        <<lastECS*millibarn<<"(mb)"<<G4endl;
#endif
    if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
    return         lastECS*millibarn;     // Use theLastElasticCS
}

GetPointer()

G4VQCrossSection * G4QIonIonCrossSection::GetPointer ( )

static

Definition at line 65 of file G4QIonIonCrossSection.cc.

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

Referenced by G4QIonIonElastic::GetMeanFreePath(), G4QLowEnergy::GetMeanFreePath(), and G4QIonIonElastic::PostStepDoIt().

ThresholdMomentum()

G4double G4QIonIonCrossSection::ThresholdMomentum	(	G4int	pZ,
		G4int	pN,
		G4int	tZ,
		G4int	tN
	)

Definition at line 519 of file G4QIonIonCrossSection.cc.

{
  static const G4double third=1./3.;
  static const G4double pM = G4QPDGCode(2212).GetMass(); // Proton mass in MeV
  static const G4double tpM= pM+pM;       // Doubled proton mass (MeV)
  if(pZ<.99 || pN<0. || tZ<.99 || tN<0.) return 0.;
  G4double tA=tZ+tN;
  G4double pA=pZ+pN;
  //G4double dE=1.263*tZ/(1.+std::pow(tA,third));
  G4double dE=pZ*tZ/(std::pow(pA,third)+std::pow(tA,third))/pA; // dE/pA (per projNucleon)
  return std::sqrt(dE*(tpM+dE));
}

---

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

• G4QIonIonCrossSection.hh
• G4QIonIonCrossSection.cc