G4QPhotonNuclearCrossSection Class Reference

#include <G4QPhotonNuclearCrossSection.hh>

Inheritance diagram for G4QPhotonNuclearCrossSection:

Public Member Functions
	~G4QPhotonNuclearCrossSection ()
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)
G4double	ThresholdEnergy (G4int Z, G4int N, G4int PDG=22)
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
	G4QPhotonNuclearCrossSection ()
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 51 of file G4QPhotonNuclearCrossSection.hh.

Constructor & Destructor Documentation

G4QPhotonNuclearCrossSection()

G4QPhotonNuclearCrossSection::G4QPhotonNuclearCrossSection ( )

inlineprotected

Definition at line 55 of file G4QPhotonNuclearCrossSection.hh.

{}

~G4QPhotonNuclearCrossSection()

G4QPhotonNuclearCrossSection::~G4QPhotonNuclearCrossSection

(

)

Definition at line 76 of file G4QPhotonNuclearCrossSection.cc.

{
  G4int lens=GDR->size();
  for(G4int i=0; i<lens; ++i) delete[] (*GDR)[i];
  delete GDR;
  G4int hens=HEN->size();
  for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
  delete HEN;
}

Member Function Documentation

CalculateCrossSection()

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

virtual

Implements G4VQCrossSection.

Definition at line 308 of file G4QPhotonNuclearCrossSection.cc.

{
#ifdef debug
  G4cout<<"G4QPhotonNucCrossSection::CalculateCrossSection: ***Called***"<<G4endl;
#endif
  static const G4double THmin=2.;    // minimum Energy Threshold
  static const G4double dE=1.;       // step for the GDR table
  static const G4int    nL=105;      // A#of GDResonance points in E (1 MeV from 2 to 106)
  static const G4double Emin=THmin+(nL-1)*dE; // minE for the HighE part
  static const G4double Emax=50000.;       // maxE for the HighE part
  static const G4int    nH=224;            // A#of HResonance points in lnE
  static const G4double milE=std::log(Emin);    // Low logarithm energy for the HighE part
  static const G4double malE=std::log(Emax);    // High logarithm energy (each 2.75 %)
  static const G4double dlE=(malE-milE)/(nH-1); // Step in log energy in the HighE part
  //
  //static const G4double shd=1.075-.0023*log(2.);  // HE PomShadowing(D)
  static const G4double shd=1.0734;  // HE PomShadowing(D)
  static const G4double shc=0.072;   // HE Shadowing constant
  static const G4double poc=0.0375;  // HE Pomeron coefficient
  static const G4double pos=16.5;    // HE Pomeron shift
  static const G4double reg=.11;     // HE Reggeon slope
  static const G4double shp=1.075;   // HE PomShadowing(P)
  //
  // Associative memory for acceleration
  static std::vector <G4double> spA; // shadowing coefficients (A-dependent)
  //
  onlyCS=CS;                         // Flag to calculate only CS (not Si/Bi)
#ifdef pdebug
  G4cout<<"G4QPhotonNucCS::CalcCS: P="<<Energy<<", F="<<F<<", I="<<I<<", Z="<<targZ
        <<", N="<<targN<<", onlyCS="<<CS<<",E="<<Energy<<",th="<<THmin<<G4endl;
  if(F==-27) return 0.;
#endif
  //if (Energy<THmin)
  //{
  //  lastE=0.;
  //  lastSig=0.;
#ifdef debug
  //  G4cout<<"---> G4QMuonNucCS::CalcCS: CS=0  as E="<<Energy<<" < "<<THmin<<G4endl;
#endif
  //  return 0.;                       // @@ This can be dangerouse for the heaviest nuc.!
  //}
  G4double sigma=0.;
  G4double A=targN+targZ;
  if(F<=0)                           // This isotope was not the last used isotop
  {
    if(F<0)                          // This isotope was found in DAMDB =-------=> RETRIEVE
    {
      lastGDR=(*GDR)[I];             // Pointer to prepared GDR cross sections
      lastHEN=(*HEN)[I];             // Pointer to prepared High Energy cross sections
      lastSP =spA[I];                // Shadowing coefficient for UHE
    }
    else                             // This isotope wasn't calculated previously => CREATE
    {
      G4double lnA=std::log(A);      // The nucleus is not found in DB. It is new.
      if(A==1.) lastSP=1.;           // The Reggeon shadowing (A=1)
      else      lastSP=A*(1.-shc*lnA); // The Reggeon shadowing
#ifdef debug
      G4cout<<">>>G4QPhotonNuclearCrossSect::CalcCS:lnA="<<lnA<<",lastSP="<<lastSP<<G4endl;
#endif
#ifdef debug3
      if(A==3) G4cout<<"G4QPhotonNuclearCrossSection::CalcCS: lastSP="<<lastSP<<G4endl;
#endif
      lastGDR = new G4double[nL];    // Allocate memory for the new GDR cross sections
      lastHEN = new G4double[nH];    // Allocate memory for the new HEN cross sections
      G4int er=GetFunctions(A,lastGDR,lastHEN);// set newZeroPosition and fill theFunctions
      if(er<1) G4cerr<<"***G4QPhotNucCrosSec::CalcCrossSection: A="<<A<<" failed"<<G4endl;
#ifdef pdebug
      G4cout<<">>G4QPhotNucCS::CalcCS:**GDR/HEN're made for A="<<A<<"** GFEr="<<er<<G4endl;
#endif
      // *** The synchronization check ***
      G4int sync=GDR->size();
      if(sync!=I) G4cerr<<"***G4QPhotoNuclearCS::CalcCS: PDG=22, S="<<sync<<"#"<<I<<G4endl;
      GDR->push_back(lastGDR);       // added GDR, found by AH 10/7/02
      HEN->push_back(lastHEN);       // added HEN, found by AH 10/7/02
      spA.push_back(lastSP);         // Pomeron Shadowing
    } // End of creation of the new set of parameters
  } // End of parameters udate
  // =-------------------= NOW the Magic Formula =--------------------------=
  if (Energy<lastTH) return 0.;             // It must be already checked in the interface
  else if (Energy<=Emin)                    // GDR region (approximated in E, not in lnE)
  {
#ifdef debug
    G4cout<<"G4QPhNCS::CalcCS:bGDR A="<<A<<", nL="<<nL<<",TH="<<THmin<<",dE="<<dE<<G4endl;
#endif
    if(A<=1. || dE <= THmin) sigma=0.;      // No GDR for A=1
    else      sigma=EquLinearFit(Energy,nL,THmin,dE,lastGDR);
#ifdef debugn
    if(sigma<0.)
      G4cout<<"G4QPhoNucCS::CalcCS:A="<<A<<",E="<<Energy<<",T="<<THmin<<",dE="<<dE<<G4endl;
#endif
  }
  else if (Energy<Emax)                     // High Energy region
  {
    G4double lE=std::log(Energy);
#ifdef pdebug
    G4cout<<"G4QPhotNucCS::CalcCS:lE="<<milE<<",n="<<nH<<",iE="<<milE<<",dE="<<dlE<<",HEN="
          <<lastHEN<<",A="<<A<<G4endl;
#endif
    if(lE > milE) sigma=EquLinearFit(lE,nH,milE,dlE,lastHEN);
    else          sigma=0.;
  }
  else                                      // UHE region (calculation, not frequent)
  {
    G4double lE=std::log(Energy);
    G4double sh=shd;
    if(A==1.)sh=shp;
    sigma=lastSP*(poc*(lE-pos)+sh*std::exp(-reg*lE));
  }
#ifdef debug
  G4cout<<"G4QPhotonNuclearCrossSection::CalcCS: sigma="<<sigma<<G4endl;
#endif
#ifdef debug
  if(Energy>45000.&&Energy<60000.)
    G4cout<<"G4QPhotoNucCS::GetCS: A="<<A<<", E="<<Energy<<",CS="<<sigma<<G4endl;
#endif
  if(sigma<0.) return 0.;
  return sigma;
}

Referenced by GetCrossSection().

GetCrossSection()

G4double G4QPhotonNuclearCrossSection::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 88 of file G4QPhotonNuclearCrossSection.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> 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<<"G4QPhCS::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)
  {
#ifdef debug
    G4cout<<"G4QPhCS::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 (?!) @@
  }
  if(pPDG==22 && tgZ==1 && !tgN && pMom<150.*MeV) return 0.; // Examle of pre-threshold (A)
  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 debug
        G4cout<<"G4QPhCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
        //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
        if(pEn<=lastTH)
        {
#ifdef debug
          G4cout<<"G4QPhCS::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(lastP == pMom)
        {
#ifdef debug
          G4cout<<"G4QPhCS::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<<"G4QPhCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
#endif
        lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
#ifdef debug
        G4cout<<"G4QPhCS::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 debug
          G4cout<<"G4QPhCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pEn;
        }
        break;                         // Go out of the LOOP
      }
#ifdef debug
      G4cout<<"---G4QPhCrossSec::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 debug
      G4cout<<"G4QPhCS::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 debug
        G4cout<<"G4QPhCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
#endif
        if(pEn>lastTH)
        {
#ifdef debug
          G4cout<<"G4QPhCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
          lastTH=pEn;
        }
      }
#ifdef debug
      G4cout<<"G4QPhCS::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 debug
      G4cout<<"G4QPhCS::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 debug
      G4cout<<"G4QPrCS::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 debug
    G4cout<<"G4QPhCS::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(lastP == pMom)
  {
#ifdef debug
    G4cout<<"G4QPhCS::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 debug
    G4cout<<"G4QPhCS::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 debug
  G4cout<<"G4QPhCS::GetCroSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
  //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
  return lastCS*millibarn;
}

GetPointer()

G4VQCrossSection * G4QPhotonNuclearCrossSection::GetPointer ( )

static

Definition at line 70 of file G4QPhotonNuclearCrossSection.cc.

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

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

ThresholdEnergy()

G4double G4QPhotonNuclearCrossSection::ThresholdEnergy ( G4int  Z, G4int  N, G4int  PDG = 22  )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 263 of file G4QPhotonNuclearCrossSection.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 134.9766; // Pi0 threshold for the nucleon
  // 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
  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 debug
  G4cout<<"G4QPhotoNucCS::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:

• G4QPhotonNuclearCrossSection.hh
• G4QPhotonNuclearCrossSection.cc