G4QNeutronElasticCrossSection Class Reference

#include <G4QNeutronElasticCrossSection.hh>

Inheritance diagram for G4QNeutronElasticCrossSection:

Public Member Functions
	~G4QNeutronElasticCrossSection ()
virtual G4double	GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=2112)
G4double	CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int pPDG, G4int Z, G4int N, G4double pP)
G4double	GetSlope (G4int tZ, G4int tN, G4int pPDG)
G4double	GetExchangeT (G4int tZ, G4int tN, G4int pPDG)
G4double	GetHMaxT ()
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
	G4QNeutronElasticCrossSection ()
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 47 of file G4QNeutronElasticCrossSection.hh.

Constructor & Destructor Documentation

G4QNeutronElasticCrossSection()

G4QNeutronElasticCrossSection::G4QNeutronElasticCrossSection ( )

protected

Definition at line 100 of file G4QNeutronElasticCrossSection.cc.

{
}

~G4QNeutronElasticCrossSection()

G4QNeutronElasticCrossSection::~G4QNeutronElasticCrossSection

(

)

Definition at line 104 of file G4QNeutronElasticCrossSection.cc.

{
  std::vector<G4double*>::iterator pos;
  for (pos=CST.begin(); pos<CST.end(); pos++)
  { delete [] *pos; }
  CST.clear();
  for (pos=PAR.begin(); pos<PAR.end(); pos++)
  { delete [] *pos; }
  PAR.clear();
  for (pos=SST.begin(); pos<SST.end(); pos++)
  { delete [] *pos; }
  SST.clear();
  for (pos=S1T.begin(); pos<S1T.end(); pos++)
  { delete [] *pos; }
  S1T.clear();
  for (pos=B1T.begin(); pos<B1T.end(); pos++)
  { delete [] *pos; }
  B1T.clear();
  for (pos=S2T.begin(); pos<S2T.end(); pos++)
  { delete [] *pos; }
  S2T.clear();
  for (pos=B2T.begin(); pos<B2T.end(); pos++)
  { delete [] *pos; }
  B2T.clear();
  for (pos=S3T.begin(); pos<S3T.end(); pos++)
  { delete [] *pos; }
  S3T.clear();
  for (pos=B3T.begin(); pos<B3T.end(); pos++)
  { delete [] *pos; }
  B3T.clear();
  for (pos=S4T.begin(); pos<S4T.end(); pos++)
  { delete [] *pos; }
  S4T.clear();
  for (pos=B4T.begin(); pos<B4T.end(); pos++)
  { delete [] *pos; }
  B4T.clear();
}

Member Function Documentation

CalculateCrossSection()

G4double G4QNeutronElasticCrossSection::CalculateCrossSection ( G4bool  CS, G4int  F, G4int  I, G4int  pPDG, G4int  Z, G4int  N, G4double  pP  )

virtual

Implements G4VQCrossSection.

Definition at line 286 of file G4QNeutronElasticCrossSection.cc.

{
  // *** Begin of Associative Memory DB for acceleration of the cross section calculations
  static std::vector <G4double>  PIN;   // Vector of max initialized log(P) in the table
  // *** End of Static Definitions (Associative Memory Data Base) ***
  G4double pMom=pIU/GeV;                // All calculations are in GeV
  onlyCS=CS;                            // Flag to calculate only CS (not Si/Bi)
#ifdef debug
  G4cout<<"G4QNeutronElasticCrosS::CalcCS:->onlyCS="<<onlyCS<<",F="<<F<<",p="<<pIU<<G4endl;
#endif
  lastLP=std::log(pMom);                // Make a logarithm of the momentum for calculation
  if(F)                                 // This isotope was found in AMDB =>RETRIEVE/UPDATE
  {
    if(F<0)                             // the AMDB must be loded
    {
      lastPIN = PIN[I];                 // Max log(P) initialised for this table set
      lastPAR = PAR[I];                 // Pointer to the parameter set
      lastCST = CST[I];                 // Pointer to the total sross-section table
      lastSST = SST[I];                 // Pointer to the first squared slope
      lastS1T = S1T[I];                 // Pointer to the first mantissa
      lastB1T = B1T[I];                 // Pointer to the first slope
      lastS2T = S2T[I];                 // Pointer to the second mantissa
      lastB2T = B2T[I];                 // Pointer to the second slope
      lastS3T = S3T[I];                 // Pointer to the third mantissa
      lastB3T = B3T[I];                 // Pointer to the rhird slope
      lastS4T = S4T[I];                 // Pointer to the 4-th mantissa
      lastB4T = B4T[I];                 // Pointer to the 4-th slope
#ifdef debug
      G4cout<<"G4QNElasticCS::CalcCS: DB is updated for I="<<I<<",*,PIN4="<<PIN[4]<<G4endl;
#endif
    }
#ifdef debug
    G4cout<<"G4QNeutronElasticCrosS::CalcCS:*read*, LP="<<lastLP<<",PIN="<<lastPIN<<G4endl;
#endif
    if(lastLP>lastPIN && lastLP<lPMax)
    {
      lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);// Can update upper logP-Limit in tabs
#ifdef debug
      G4cout<<"G4QNElCrS::CalcCS:updated(I),LP="<<lastLP<<"<IN["<<I<<"]="<<lastPIN<<G4endl;
#endif
      PIN[I]=lastPIN;                   // Remember the new P-Limit of the tables
    }
  }
  else                                  // This isotope wasn't initialized => CREATE
  {
    lastPAR = new G4double[nPoints];    // Allocate memory for parameters of CS function
    lastPAR[nLast]=0;                   // Initialization for VALGRIND
    lastCST = new G4double[nPoints];    // Allocate memory for Tabulated CS function    
    lastSST = new G4double[nPoints];    // Allocate memory for Tabulated first sqaredSlope 
    lastS1T = new G4double[nPoints];    // Allocate memory for Tabulated first mantissa 
    lastB1T = new G4double[nPoints];    // Allocate memory for Tabulated first slope    
    lastS2T = new G4double[nPoints];    // Allocate memory for Tabulated second mantissa
    lastB2T = new G4double[nPoints];    // Allocate memory for Tabulated second slope   
    lastS3T = new G4double[nPoints];    // Allocate memory for Tabulated third mantissa 
    lastB3T = new G4double[nPoints];    // Allocate memory for Tabulated third slope    
    lastS4T = new G4double[nPoints];    // Allocate memory for Tabulated 4-th mantissa 
    lastB4T = new G4double[nPoints];    // Allocate memory for Tabulated 4-th slope    
#ifdef debug
    G4cout<<"G4QNeutronElasticCrosS::CalcCS:*ini*,lastLP="<<lastLP<<",min="<<lPMin<<G4endl;
#endif
    lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN); // Returns the new P-limit for tables
#ifdef debug
    G4cout<<"G4QNElCS::CalcCS:i,Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<",LP"<<lastPIN<<G4endl;
#endif
    PIN.push_back(lastPIN);             // Fill parameters of CS function to AMDB
    PAR.push_back(lastPAR);             // Fill parameters of CS function to AMDB
    CST.push_back(lastCST);             // Fill Tabulated CS function to AMDB    
    SST.push_back(lastSST);             // Fill Tabulated first sq.slope to AMDB 
    S1T.push_back(lastS1T);             // Fill Tabulated first mantissa to AMDB 
    B1T.push_back(lastB1T);             // Fill Tabulated first slope to AMDB    
    S2T.push_back(lastS2T);             // Fill Tabulated second mantissa to AMDB 
    B2T.push_back(lastB2T);             // Fill Tabulated second slope to AMDB    
    S3T.push_back(lastS3T);             // Fill Tabulated third mantissa to AMDB 
    B3T.push_back(lastB3T);             // Fill Tabulated third slope to AMDB    
    S4T.push_back(lastS4T);             // Fill Tabulated 4-th mantissa to AMDB 
    B4T.push_back(lastB4T);             // Fill Tabulated 4-th slope to AMDB    
  } // End of creation/update of the new set of parameters and tables
  // =-------= NOW Update (if necessary) and Calculate the Cross Section =---------=
#ifdef debug
  G4cout<<"G4QNElCrS::CalcCS:?update?,LP="<<lastLP<<",IN="<<lastPIN<<",ML="<<lPMax<<G4endl;
#endif
  if(lastLP>lastPIN && lastLP<lPMax)
  {
    lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
#ifdef debug
    G4cout<<"G4QNeutronElCrS::CalcCS:*updated(O)*, LP="<<lastLP<<" < IN="<<lastPIN<<G4endl;
#endif
  }
#ifdef debug
  G4cout<<"G4QNEltCrS::CalcCS: lastLP="<<lastLP<<",lPM="<<lPMin<<",lPIN="<<lastPIN<<G4endl;
#endif
  if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom); // Calculate (-t)_max=Q2_max (GeV2)
#ifdef debug
  G4cout<<"G4QNeutElastCroSec::CalcCS:oCS="<<onlyCS<<",-t="<<lastTM<<",p="<<lastLP<<G4endl;
#endif
  if(lastLP>lPMin && lastLP<=lastPIN)   // Linear fit is made using precalculated tables
  {
    if(lastLP==lastPIN)
    {
      G4double shift=(lastLP-lPMin)/dlnP+.000001; // Log distance from lPMin
      G4int    blast=static_cast<int>(shift); // this is a bin number of the lower edge (0)
      if(blast<0 || blast>=nLast) G4cout<<"G4QNeutElCS::CCS:b="<<blast<<","<<nLast<<G4endl;
      lastSIG = lastCST[blast];
      if(!onlyCS)                       // Skip the differential cross-section parameters
      {
        theSS  = lastSST[blast];
        theS1  = lastS1T[blast];
        theB1  = lastB1T[blast];
        theS2  = lastS2T[blast];
        theB2  = lastB2T[blast];
        theS3  = lastS3T[blast];
        theB3  = lastB3T[blast];
        theS4  = lastS4T[blast];
        theB4  = lastB4T[blast];
      }
#ifdef debug
      G4cout<<"G4QNeutronElasticCrossSection::CalcCS:(E)S1="<<theS1<<",B1="<<theB1<<G4endl;
#endif
    }
    else
    {
      G4double shift=(lastLP-lPMin)/dlnP;        // a shift from the beginning of the table
      G4int    blast=static_cast<int>(shift);    // the lower bin number
      if(blast<0)   blast=0;
      if(blast>=nLast) blast=nLast-1;            // low edge of the last bin
      shift-=blast;                              // step inside the unit bin
      G4int lastL=blast+1;                       // the upper bin number
      G4double SIGL=lastCST[blast];              // the basic value of the cross-section
      lastSIG= SIGL+shift*(lastCST[lastL]-SIGL); // calculated total elastic cross-section
#ifdef debug
      G4cout<<"G4QNeutronElCS::CalcCrossSection: Sig="<<lastSIG<<", P="<<pMom<<", Z="<<tgZ
            <<", N="<<tgN<<", PDG="<<PDG<<", onlyCS="<<onlyCS<<G4endl;
#endif
      if(!onlyCS)                       // Skip the differential cross-section parameters
      {
        G4double SSTL=lastSST[blast];           // the low bin of the first squared slope
        theSS=SSTL+shift*(lastSST[lastL]-SSTL); // the basic value of the first sq.slope
        G4double S1TL=lastS1T[blast];           // the low bin of the first mantissa
        theS1=S1TL+shift*(lastS1T[lastL]-S1TL); // the basic value of the first mantissa
        G4double B1TL=lastB1T[blast];           // the low bin of the first slope
#ifdef debug
        G4cout<<"G4QNeutronElCrS::CalcCrossSection:bl="<<blast<<",ls="<<lastL<<",SL="<<S1TL
              <<",SU="<<lastS1T[lastL]<<",BL="<<B1TL<<",BU="<<lastB1T[lastL]<<G4endl;
#endif
        theB1=B1TL+shift*(lastB1T[lastL]-B1TL); // the basic value of the first slope
        G4double S2TL=lastS2T[blast];           // the low bin of the second mantissa
        theS2=S2TL+shift*(lastS2T[lastL]-S2TL); // the basic value of the second mantissa
        G4double B2TL=lastB2T[blast];           // the low bin of the second slope
        theB2=B2TL+shift*(lastB2T[lastL]-B2TL); // the basic value of the second slope
        G4double S3TL=lastS3T[blast];           // the low bin of the third mantissa
        theS3=S3TL+shift*(lastS3T[lastL]-S3TL); // the basic value of the third mantissa
#ifdef debug
        G4cout<<"G4QNElCS::CCS: s3l="<<S3TL<<",sh3="<<shift<<",s3h="<<lastS3T[lastL]<<",b="
              <<blast<<",l="<<lastL<<G4endl;
#endif
        G4double B3TL=lastB3T[blast];           // the low bin of the third slope
        theB3=B3TL+shift*(lastB3T[lastL]-B3TL); // the basic value of the third slope
        G4double S4TL=lastS4T[blast];           // the low bin of the 4-th mantissa
        theS4=S4TL+shift*(lastS4T[lastL]-S4TL); // the basic value of the 4-th mantissa
#ifdef debug
        G4cout<<"G4QNElCS::CCS: s4l="<<S4TL<<",sh4="<<shift<<",s4h="<<lastS4T[lastL]<<",b="
              <<blast<<",l="<<lastL<<G4endl;
#endif
        G4double B4TL=lastB4T[blast];           // the low bin of the 4-th slope
        theB4=B4TL+shift*(lastB4T[lastL]-B4TL); // the basic value of the 4-th slope
      }
#ifdef debug
      G4cout<<"G4QNeutronElasticCrossSection::CalcCS:(I)S1="<<theS1<<",B1="<<theB1<<G4endl;
#endif
    }
  }
  else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN); // Direct calculation beyond the table
  if(lastSIG<0.) lastSIG = 0.;                   // @@ a Warning print can be added
#ifdef debug
  G4cout<<"G4QNeutronElasticCrossSection::CalculateCS: END, onlyCS="<<onlyCS<<G4endl;
#endif
  return lastSIG;
}

Referenced by GetCrossSection().

GetCrossSection()

G4double G4QNeutronElasticCrossSection::GetCrossSection ( G4bool  fCS, G4double  pMom, G4int  tgZ, G4int  tgN, G4int  pPDG = 2112  )

virtual

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

Reimplemented from G4VQCrossSection.

Definition at line 151 of file G4QNeutronElasticCrossSection.cc.

{
  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;
  onlyCS=fCS;
#ifdef debug
  G4cout<<"G4QNElCS::GetCS:>>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="
        <<tgN<<"("<<lastN<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="<<colN.size()<<G4endl;
  //CalculateCrossSection(fCS,-27,j,pPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
  if(pPDG!=2112)
  {
    G4cout<<"G4QNeutronElCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
    //CalculateCrossSection(fCS,-27,j,pPDG,lastZ,lastN,pMom); // DUMMY TEST
    return 0.;                         // projectile PDG=0 is a mistake (?!) @@
  }
  G4bool in=false;                     // By default the isotope must be found in the AMDB
  lastP   = 0.;                      // New momentum history (nothing to compare with)
  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
  if(lastI) for(G4int i=0; i<lastI; i++) // Loop over proj/tgZ/tgN lines of DB
  {                                  // The nucleus with projPDG is found in AMDB
    if(colN[i]==tgN && colZ[i]==tgZ) // Isotope is foind in AMDB
    {
      lastI=i;
      lastTH =colTH[i];              // Last THreshold (A-dependent)
#ifdef debug
      G4cout<<"G4QNeutElCS::GetCS:Found,P="<<pMom<<",Threshold="<<lastTH<<",i="<<i<<G4endl;
      //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
      if(pEn<=lastTH)
      {
#ifdef debug
        G4cout<<"G4QNeutElCS::GetCS:Found,T="<<pEn<<"<Threshold="<<lastTH<<",CS=0"<<G4endl;
        //CalculateCrossSection(fCS,-27,i,pPDG,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) //VI (do not use tolerance)
      if(lastP == pMom)              // Do not recalculate
      {
#ifdef debug
        G4cout<<"G4QNeutronElasticCS::GetCrosS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
#endif
        CalculateCrossSection(fCS,-1,i,pPDG,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 debug
      G4cout<<"G4QNElCrS::G:UpdateDB,P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",i="<<i<<G4endl;
#endif
      lastCS=CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom); // read & update
#ifdef debug
      G4cout<<"G4QNeutronElCS::GetCrosSec:*****>New (inDB) Calculated CS="<<lastCS<<G4endl;
      //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
      if(lastCS<=0. && pEn>lastTH)    // Correct the threshold
      {
#ifdef debug
        G4cout<<"G4QNeutronElCS::GetCS:New,T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
        lastTH=pEn;
      }
      break;                           // Go out of the LOOP with found lastI
    }
#ifdef debug
    G4cout<<"---G4QNeutronElasticCrossSection::GetCrosSec:pPDG="<<pPDG<<",i="<<i<<",N="
          <<colN[i]<<",Z["<<i<<"]="<<colZ[i]<<G4endl;
    //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
  }
  if(!in)                            // This nucleus has not been calculated previously
  {
#ifdef debug
    G4cout<<"G4QNElCrS::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,lastI,pPDG,lastZ,lastN,pMom);//calculate&create
    if(lastCS<=0.)
    {
      lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
#ifdef debug
      G4cout<<"G4QNeutronElCrosSect::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
#endif
      if(pEn>lastTH)
      {
#ifdef debug
        G4cout<<"G4QNeutElCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
#endif
        lastTH=pEn;
      }
    }
#ifdef debug
    G4cout<<"G4QNElCrS::GetCrosSec: New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
    //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
    colN.push_back(tgN);
    colZ.push_back(tgZ);
    colP.push_back(pMom);
    colTH.push_back(lastTH);
    colCS.push_back(lastCS);
#ifdef debug
    G4cout<<"G4QNElCrS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
    //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
#endif
    return lastCS*millibarn;
  } // End of creation of the new set of parameters
  else
  {
#ifdef debug
    G4cout<<"G4QNeutronElasticCrossSection::GetCS: Update lastI="<<lastI<<G4endl;
#endif
    colP[lastI]=pMom;
    colCS[lastI]=lastCS;
  }
#ifdef debug
  G4cout<<"G4QNElCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
  //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
  G4cout<<"G4QNeutronElasticCrossSection::GetCrSec:***End***, onlyCS="<<onlyCS<<G4endl;
#endif
  return lastCS*millibarn;
}

GetExchangeT()

G4double G4QNeutronElasticCrossSection::GetExchangeT ( G4int  tZ, G4int  tN, G4int  pPDG  )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 1997 of file G4QNeutronElasticCrossSection.cc.

{
  static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
  static const G4double third=1./3.;
  static const G4double fifth=1./5.;
  static const G4double sevth=1./7.;
#ifdef tdebug
  G4cout<<"G4QNeutElCrS::GetExcT:F="<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
#endif
  if(PDG!=2112) G4cout<<"*Warning*G4QNeutronElasticCrossSection::GetExT:PDG="<<PDG<<G4endl;
  if(onlyCS) G4cout<<"*Warning*G4QNeutronElasticCrossSection::GetExchangeT:onCS=1"<<G4endl;
  if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();// S-wave for p<14 MeV/c (kinE<.1MeV)
  G4double q2=0.;
  if(tgZ==1 && tgN==0)                // ===> n+p=n+p
  {
#ifdef tdebug
    G4cout<<"G4QNeutronElasticCrS::GetExchangeT: TM="<<lastTM<<",S1="<<theS1<<",B1="<<theB1
          <<",S2="<<theS2<<",B2="<<theB2<<",GeV2="<<GeVSQ<<G4endl;
#endif
    G4double E1=lastTM*theB1;
    G4double R1=(1.-std::exp(-E1));
#ifdef tdebug
    G4double ts1=-std::log(1.-R1)/theB1;
    G4double ds1=std::fabs(ts1-lastTM)/lastTM;
    if(ds1>.0001)
      G4cout<<"*Warning*G4QNeutronElasticCrossSection::GetExT:1n "<<ts1<<"#"<<lastTM<<",d="
            <<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl;
#endif
    G4double E2=lastTM*theB2;
    G4double R2=(1.-std::exp(-E2));
#ifdef tdebug
    G4double ts2=-std::log(1.-R2)/theB2;
    G4double ds2=std::fabs(ts2-lastTM)/lastTM;
    if(ds2>.0001)
      G4cout<<"*Warning*G4QNeutronElasticCrossSection::GetExT:2n "<<ts2<<"#"<<lastTM<<",d="
            <<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl;
#endif
    //G4double E3=lastTM*theB3;
    //G4double R3=(1.-std::exp(-E3));
#ifdef tdebug
    //G4double ts3=-std::log(1.-R3)/theB3;
    //G4double ds3=std::fabs(ts3-lastTM)/lastTM;
    //if(ds3>.01)G4cout<<"Warn*G4QNElCS::GetExT:3n="<<ts3<<"#"<<lastTM<<",d="<<ds3<<G4endl;
#endif
    G4double I1=R1*theS1;
    G4double I2=R2*theS2/theB2;
    //G4double I3=R3*theS3/theB3;
    G4double I12=I1+I2;
    //G4double rand=(I12+I3)*G4UniformRand();
    G4double rand=I12*G4UniformRand();
    if     (rand<I1 )
    {
      G4double ran=R1*G4UniformRand();
      if(ran>1.) ran=1.;
      q2=-std::log(1.-ran)/theB1;       // t-chan
    }
    else
    {
      G4double ran=R2*G4UniformRand();
      if(ran>1.) ran=1.;
      q2=lastTM+std::log(1.-ran)/theB2; // u-chan (ChEx)
    }
  }
  else
  {
    G4double a=tgZ+tgN;
#ifdef tdebug
    G4cout<<"G4QNeutronElCroSec::GetExT:a="<<a<<",t="<<lastTM<<",S1="<<theS1<<",B1="<<theB1
          <<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3
          <<",S4="<<theS4<<",B4="<<theB4<<G4endl;
#endif
    G4double E1=lastTM*(theB1+lastTM*theSS);
    G4double R1=(1.-std::exp(-E1));
    G4double tss=theSS+theSS; // for future solution of quadratic equation (imediate check)
#ifdef tdebug
    G4double ts1=-std::log(1.-R1)/theB1;
    if(std::fabs(tss)>1.e-7) ts1=(std::sqrt(theB1*(theB1+(tss+tss)*ts1))-theB1)/tss;
    G4double ds1=(ts1-lastTM)/lastTM;
    if(ds1>.0001)
      G4cout<<"*Warning*G4QNeutronElasticCrossSection::GetExT:1a "<<ts1<<"#"<<lastTM<<",d="
            <<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl;
#endif
    G4double tm2=lastTM*lastTM;
    G4double E2=lastTM*tm2*theB2;                   // power 3 for lowA, 5 for HighA (1st)
    if(a>6.5)E2*=tm2;                               // for heavy nuclei
    G4double R2=(1.-std::exp(-E2));
#ifdef tdebug
    G4double ts2=-std::log(1.-R2)/theB2;
    if(a<6.5)ts2=std::pow(ts2,third);
    else     ts2=std::pow(ts2,fifth);
    G4double ds2=std::fabs(ts2-lastTM)/lastTM;
    if(ds2>.0001)
      G4cout<<"*Warning*G4QNeutronElasticCrossSection::GetExT:2a "<<ts2<<"#"<<lastTM<<",d="
            <<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl;
#endif
    G4double E3=lastTM*theB3;
    if(a>6.5)E3*=tm2*tm2*tm2;                       // power 1 for lowA, 7 (2nd) for HighA
    G4double R3=(1.-std::exp(-E3));
#ifdef tdebug
    G4double ts3=-std::log(1.-R3)/theB3;
    if(a>6.5)ts3=std::pow(ts3,sevth);
    G4double ds3=std::fabs(ts3-lastTM)/lastTM;
    if(ds3>.0001)
      G4cout<<"*Warning*G4QNeutronElasticCrossSection::GetExT:3a "<<ts3<<"#"<<lastTM<<",d="
            <<ds3<<",R3="<<R3<<",E3="<<E3<<G4endl;
#endif
    G4double E4=lastTM*theB4;
    G4double R4=(1.-std::exp(-E4));
#ifdef tdebug
    G4double ts4=-std::log(1.-R4)/theB4;
    G4double ds4=std::fabs(ts4-lastTM)/lastTM;
    if(ds4>.0001)
      G4cout<<"*Warning*G4QNeutronElasticCrissSection::GetExT:4a "<<ts4<<"#"<<lastTM<<",d="
            <<ds4<<",R4="<<R4<<",E4="<<E4<<G4endl;
#endif
    G4double I1=R1*theS1;
    G4double I2=R2*theS2;
    G4double I3=R3*theS3;
    G4double I4=R4*theS4;
    G4double I12=I1+I2;
    G4double I13=I12+I3;
    G4double rand=(I13+I4)*G4UniformRand();
#ifdef tdebug
    G4cout<<"G4QNElCS::GtExT:1="<<I1<<",2="<<I2<<",3="<<I3<<",4="<<I4<<",r="<<rand<<G4endl;
#endif
    if(rand<I1)
    {
      G4double ran=R1*G4UniformRand();
      if(ran>1.) ran=1.;
      q2=-std::log(1.-ran)/theB1;
      if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
#ifdef tdebug
      G4cout<<"G4QNElCS::GetET:Q2="<<q2<<",ss="<<tss/2<<",b1="<<theB1<<",t1="<<ts1<<G4endl;
#endif
    }
    else if(rand<I12)
    {
      G4double ran=R2*G4UniformRand();
      if(ran>1.) ran=1.;
      q2=-std::log(1.-ran)/theB2;
      if(q2<0.) q2=0.;
      if(a<6.5) q2=std::pow(q2,third);
      else      q2=std::pow(q2,fifth);
#ifdef tdebug
      G4cout<<"G4QNElCS::GetExT: Q2="<<q2<<",r2="<<R2<<", b2="<<theB2<<",t2="<<ts2<<G4endl;
#endif
    }
    else if(rand<I13)
    {
      G4double ran=R3*G4UniformRand();
      if(ran>1.) ran=1.;
      q2=-std::log(1.-ran)/theB3;
      if(q2<0.) q2=0.;
      if(a>6.5) q2=std::pow(q2,sevth);
#ifdef tdebug
      G4cout<<"G4QNElCS::GetExT:Q2="<<q2<<", r3="<<R2<<", b3="<<theB2<<",t3="<<ts2<<G4endl;
#endif
    }
    else
    {
      G4double ran=R4*G4UniformRand();
      if(ran>1.) ran=1.;
      q2=-std::log(1.-ran)/theB4;
      if(a<6.5) q2=lastTM-q2;                    // u reduced for lightA (starts from 0)
#ifdef tdebug
      G4cout<<"G4QNElCS::GetET:Q2="<<q2<<",m="<<lastTM<<",b4="<<theB3<<",t4="<<ts3<<G4endl;
#endif
    }
  }
  if(q2<0.) q2=0.;
  if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QNeutronElCroSect::GetExchangeT: -t="<<q2<<G4endl;
  if(q2>lastTM)
  {
#ifdef tdebug
    G4cout<<"*Warning*G4QNeutronElasticCrossSection::GetExT:-t="<<q2<<" >"<<lastTM<<G4endl;
#endif
    q2=lastTM;
  }
  return q2*GeVSQ;
}

GetHMaxT()

G4double G4QNeutronElasticCrossSection::GetHMaxT ( )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 2204 of file G4QNeutronElasticCrossSection.cc.

{
  static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
  return lastTM*HGeVSQ;
}

GetPointer()

G4VQCrossSection * G4QNeutronElasticCrossSection::GetPointer ( )

static

Definition at line 143 of file G4QNeutronElasticCrossSection.cc.

{
  static G4QNeutronElasticCrossSection theCrossSection;//*StatBody of the QEl CrossSection*
  return &theCrossSection;
}

Referenced by G4QuasiFreeRatios::ChExer(), G4CHIPSElastic::G4CHIPSElastic(), G4CHIPSElasticXS::G4CHIPSElasticXS(), G4QHadronElasticDataSet::GetIsoCrossSection(), G4QElastic::GetMeanFreePath(), G4QElastic::PostStepDoIt(), G4QIonIonElastic::PostStepDoIt(), G4QLowEnergy::PostStepDoIt(), and G4QuasiFreeRatios::Scatter().

GetSlope()

G4double G4QNeutronElasticCrossSection::GetSlope ( G4int  tZ, G4int  tN, G4int  pPDG  )

virtual

Reimplemented from G4VQCrossSection.

Definition at line 2179 of file G4QNeutronElasticCrossSection.cc.

{
  static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
#ifdef tdebug
  G4cout<<"G4QNeutrElCS::GetSlope:"<<onlyCS<<", Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
#endif
  if(onlyCS) G4cout<<"Warning*G4QNeutronElasticCrossSection::GetSlope:onlyCS=true"<<G4endl;
  if(lastLP<-4.3) return 0.;          // S-wave for p<14 MeV/c (kinE<.1MeV)
  if(PDG!=2112)
  {
    // G4cout<<"*Error*G4QNeutronElasticCrossSection::GetSlope: PDG="<<PDG<<", Z="<<tgZ
    //       <<", N="<<tgN<<", while it is defined only for PDG=2112"<<G4endl;
    // throw G4QException("G4QNeutronElasticCrossSection::GetSlope: only nA are implemented");
    G4ExceptionDescription ed;
    ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
       <<", while it is defined only for PDG=2112 (n) " << G4endl;
    G4Exception("G4QNeutronElasticCrossSection::GetSlope()", "HAD_CHPS_0000",
                FatalException, ed);
  }
  if(theB1<0.) theB1=0.;
  if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QNeutElasticCrosS::Getslope:"<<theB1<<G4endl;
  return theB1/GeVSQ;
}

---

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

• G4QNeutronElasticCrossSection.hh
• G4QNeutronElasticCrossSection.cc