G4ComponentGGHadronNucleusXsc Class Reference

#include <G4ComponentGGHadronNucleusXsc.hh>

Inheritance diagram for G4ComponentGGHadronNucleusXsc:

Public Member Functions
	G4ComponentGGHadronNucleusXsc ()
virtual	~G4ComponentGGHadronNucleusXsc ()
virtual G4double	GetTotalIsotopeCrossSection (const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
virtual G4double	GetTotalElementCrossSection (const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
virtual G4double	GetInelasticIsotopeCrossSection (const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
virtual G4double	GetInelasticElementCrossSection (const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
virtual G4double	GetElasticElementCrossSection (const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
virtual G4double	GetElasticIsotopeCrossSection (const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
virtual G4double	ComputeQuasiElasticRatio (const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
G4bool	IsIsoApplicable (const G4DynamicParticle *aDP, G4int Z, G4int A, const G4Element *elm=0, const G4Material *mat=0)
G4double	GetIsoCrossSection (const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
G4double	GetRatioSD (const G4DynamicParticle *, G4int At, G4int Zt)
G4double	GetRatioQE (const G4DynamicParticle *, G4int At, G4int Zt)
G4double	GetHadronNucleonXsc (const G4DynamicParticle *, const G4Element *)
G4double	GetHadronNucleonXsc (const G4DynamicParticle *, G4int At, G4int Zt)
G4double	GetHadronNucleonXscPDG (const G4DynamicParticle *, const G4Element *)
G4double	GetHadronNucleonXscPDG (const G4DynamicParticle *, G4int At, G4int Zt)
G4double	GetHadronNucleonXscNS (const G4DynamicParticle *, const G4Element *)
G4double	GetHadronNucleonXscNS (const G4DynamicParticle *, G4int At, G4int Zt)
G4double	GetKaonNucleonXscVector (const G4DynamicParticle *, G4int At, G4int Zt)
G4double	GetHNinelasticXsc (const G4DynamicParticle *, const G4Element *)
G4double	GetHNinelasticXsc (const G4DynamicParticle *, G4int At, G4int Zt)
G4double	GetHNinelasticXscVU (const G4DynamicParticle *, G4int At, G4int Zt)
G4double	CalculateEcmValue (const G4double, const G4double, const G4double)
G4double	CalcMandelstamS (const G4double, const G4double, const G4double)
G4double	GetNucleusRadius (const G4DynamicParticle *, const G4Element *)
G4double	GetNucleusRadius (G4int At)
virtual void	CrossSectionDescription (std::ostream &) const
G4double	GetElasticGlauberGribov (const G4DynamicParticle *, G4int Z, G4int A)
G4double	GetInelasticGlauberGribov (const G4DynamicParticle *, G4int Z, G4int A)
G4double	GetTotalGlauberGribovXsc ()
G4double	GetElasticGlauberGribovXsc ()
G4double	GetInelasticGlauberGribovXsc ()
G4double	GetProductionGlauberGribovXsc ()
G4double	GetDiffractionGlauberGribovXsc ()
G4double	GetRadiusConst ()
G4double	GetParticleBarCorTot (const G4ParticleDefinition *theParticle, G4int Z)
G4double	GetParticleBarCorIn (const G4ParticleDefinition *theParticle, G4int Z)
void	SetEnergyLowerLimit (G4double E)

---

Detailed Description

Definition at line 51 of file G4ComponentGGHadronNucleusXsc.hh.

---

Constructor & Destructor Documentation

G4ComponentGGHadronNucleusXsc::G4ComponentGGHadronNucleusXsc

(

)

Definition at line 44 of file G4ComponentGGHadronNucleusXsc.cc.

References G4Alpha::Alpha(), G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), G4Deuteron::Deuteron(), G4Gamma::Gamma(), G4He3::He3(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4Neutron::Neutron(), G4OmegaMinus::OmegaMinus(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4Proton::Proton(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), G4SigmaZero::SigmaZero(), G4Triton::Triton(), G4XiMinus::XiMinus(), and G4XiZero::XiZero().

 : G4VComponentCrossSection("Glauber-Gribov"),
   fUpperLimit(100000*GeV), fLowerLimit(10.*MeV),// fLowerLimit(3*GeV),
   fRadiusConst(1.08*fermi),  // 1.1, 1.3 ?
   fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
   fDiffractionXsc(0.0), fHadronNucleonXsc(0.0)
{
  theGamma    = G4Gamma::Gamma();
  theProton   = G4Proton::Proton();
  theNeutron  = G4Neutron::Neutron();
  theAProton  = G4AntiProton::AntiProton();
  theANeutron = G4AntiNeutron::AntiNeutron();
  thePiPlus   = G4PionPlus::PionPlus();
  thePiMinus  = G4PionMinus::PionMinus();
  thePiZero   = G4PionZero::PionZero();
  theKPlus    = G4KaonPlus::KaonPlus();
  theKMinus   = G4KaonMinus::KaonMinus();
  theK0S      = G4KaonZeroShort::KaonZeroShort();
  theK0L      = G4KaonZeroLong::KaonZeroLong();
  theL        = G4Lambda::Lambda();
  theAntiL    = G4AntiLambda::AntiLambda();
  theSPlus    = G4SigmaPlus::SigmaPlus();
  theASPlus   = G4AntiSigmaPlus::AntiSigmaPlus();
  theSMinus   = G4SigmaMinus::SigmaMinus();
  theASMinus  = G4AntiSigmaMinus::AntiSigmaMinus();
  theS0       = G4SigmaZero::SigmaZero();
  theAS0      = G4AntiSigmaZero::AntiSigmaZero();
  theXiMinus  = G4XiMinus::XiMinus();
  theXi0      = G4XiZero::XiZero();
  theAXiMinus = G4AntiXiMinus::AntiXiMinus();
  theAXi0     = G4AntiXiZero::AntiXiZero();
  theOmega    = G4OmegaMinus::OmegaMinus();
  theAOmega   = G4AntiOmegaMinus::AntiOmegaMinus();
  theD        = G4Deuteron::Deuteron();
  theT        = G4Triton::Triton();
  theA        = G4Alpha::Alpha();
  theHe3      = G4He3::He3();

  hnXsc = new G4HadronNucleonXsc();
}

G4ComponentGGHadronNucleusXsc::~G4ComponentGGHadronNucleusXsc

(

)

[virtual]

Definition at line 89 of file G4ComponentGGHadronNucleusXsc.cc.

{
  if (hnXsc) delete hnXsc;
}

---

Member Function Documentation

G4double G4ComponentGGHadronNucleusXsc::CalcMandelstamS	(	const	G4double,
		const	G4double,
		const	G4double
	)

Definition at line 1428 of file G4ComponentGGHadronNucleusXsc.cc.

Referenced by GetHadronNucleonXsc(), GetHadronNucleonXscNS(), and GetHadronNucleonXscPDG().

{
  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
  G4double sMand  = mp*mp + mt*mt + 2*Elab*mt ;

  return sMand;
}

G4double G4ComponentGGHadronNucleusXsc::CalculateEcmValue	(	const	G4double,
		const	G4double,
		const	G4double
	)

Definition at line 1412 of file G4ComponentGGHadronNucleusXsc.cc.

{
  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
  G4double Ecm  = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
  // G4double Pcm  = Plab * mt / Ecm;
  // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;

  return Ecm ; // KEcm;
}

G4double G4ComponentGGHadronNucleusXsc::ComputeQuasiElasticRatio	(	const G4ParticleDefinition *	aParticle,
		G4double	kinEnergy,
		G4int	Z,
		G4int	A
	)			[virtual]

Reimplemented from G4VComponentCrossSection.

Definition at line 180 of file G4ComponentGGHadronNucleusXsc.cc.

References GetIsoCrossSection().

{
  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.), 
                                                kinEnergy);
  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
  delete aDP;
  G4double ratio = 0.;

  if(fInelasticXsc > 0.)
  {
    ratio = (fInelasticXsc - fProductionXsc)/fInelasticXsc;
    if(ratio < 0.) ratio = 0.;
  }
  return ratio;
}

void G4ComponentGGHadronNucleusXsc::CrossSectionDescription

(

std::ostream &

)

const [virtual]

Definition at line 1442 of file G4ComponentGGHadronNucleusXsc.cc.

{
  outFile << "G4ComponentGGHadronNucleusXsc calculates total, inelastic and\n"
          << "elastic cross sections for hadron-nucleus cross sections using\n"
          << "the Glauber model with Gribov corrections.  It is valid for all\n"
          << "targets except hydrogen, and for incident p, pbar, n, sigma-,\n"
          << "pi+, pi-, K+, K- and gammas with energies above 3 GeV.  This is\n"
          << "a cross section component which is to be used to build a cross\n"
          << "data set.\n";
}

G4double G4ComponentGGHadronNucleusXsc::GetDiffractionGlauberGribovXsc

(

)

[inline]

Definition at line 141 of file G4ComponentGGHadronNucleusXsc.hh.

{ return fDiffractionXsc; }; 

G4double G4ComponentGGHadronNucleusXsc::GetElasticElementCrossSection	(	const G4ParticleDefinition *	aParticle,
		G4double	kinEnergy,
		G4int	Z,
		G4double	A
	)			[virtual]

Implements G4VComponentCrossSection.

Definition at line 152 of file G4ComponentGGHadronNucleusXsc.cc.

References GetIsoCrossSection().

{
  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.), 
                                                kinEnergy);
  fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
  delete aDP;

  return fElasticXsc;
}

G4double G4ComponentGGHadronNucleusXsc::GetElasticGlauberGribov	(	const G4DynamicParticle *	,
		G4int	Z,
		G4int	A
	)			[inline]

Definition at line 211 of file G4ComponentGGHadronNucleusXsc.hh.

References GetIsoCrossSection().

{
  GetIsoCrossSection(dp, Z, A);
  return fElasticXsc;
}

G4double G4ComponentGGHadronNucleusXsc::GetElasticGlauberGribovXsc

(

)

[inline]

Definition at line 138 of file G4ComponentGGHadronNucleusXsc.hh.

{ return fElasticXsc;   }; 

G4double G4ComponentGGHadronNucleusXsc::GetElasticIsotopeCrossSection	(	const G4ParticleDefinition *	aParticle,
		G4double	kinEnergy,
		G4int	Z,
		G4int	A
	)			[virtual]

Implements G4VComponentCrossSection.

Definition at line 166 of file G4ComponentGGHadronNucleusXsc.cc.

References GetIsoCrossSection().

{
  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.), 
                                                kinEnergy);
  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
  delete aDP;

  return fElasticXsc;
}

G4double G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc	(	const G4DynamicParticle *	,
		G4int	At,
		G4int	Zt
	)

Definition at line 460 of file G4ComponentGGHadronNucleusXsc.cc.

References CalcMandelstamS(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMass(), and G4DynamicParticle::GetMomentum().

{
  G4double xsection;

  //G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);

  G4double targ_mass = 0.939*GeV;  // ~mean neutron and proton ???

  G4double proj_mass     = aParticle->GetMass();
  G4double proj_momentum = aParticle->GetMomentum().mag();
  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );

  sMand /= GeV*GeV;  // in GeV for parametrisation
  proj_momentum /= GeV;

  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
  
  G4double aa = At;

  if(theParticle == theGamma) 
  {
    xsection = aa*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
  } 
  else if(theParticle == theNeutron) // as proton ??? 
  {
    xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
  } 
  else if(theParticle == theProton) 
  {
    xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
    // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
  } 
  else if(theParticle == theAProton) 
  {
    xsection = aa*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
  } 
  else if(theParticle == thePiPlus) 
  {
    xsection = aa*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
  } 
  else if(theParticle == thePiMinus) 
  {
    // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
    xsection = aa*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
  } 
  else if(theParticle == theKPlus) 
  {
    xsection = aa*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
  } 
  else if(theParticle == theKMinus) 
  {
    xsection = aa*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
  }
  else  // as proton ??? 
  {
    xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
  } 
  xsection *= millibarn;
  return xsection;
}

G4double G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc	(	const G4DynamicParticle *	,
		const G4Element *
	)

Definition at line 443 of file G4ComponentGGHadronNucleusXsc.cc.

References G4lrint(), G4Element::GetN(), and G4Element::GetZ().

{
  G4int At = G4lrint(anElement->GetN());  // number of nucleons 
  G4int Zt = G4lrint(anElement->GetZ());  // number of protons

  return GetHadronNucleonXsc(aParticle, At, Zt);
}

G4double G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS	(	const G4DynamicParticle *	,
		G4int	At,
		G4int	Zt
	)

Definition at line 682 of file G4ComponentGGHadronNucleusXsc.cc.

References CalcMandelstamS(), G4DynamicParticle::GetDefinition(), GetHadronNucleonXscPDG(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), G4ParticleTable::GetParticleTable(), G4DynamicParticle::GetTotalEnergy(), G4InuclParticleNames::nn, and G4InuclParticleNames::s0.

{
  G4double xsection(0);
  // G4double Delta;   DHW 19 May 2011: variable set but not used
  G4double A0, B0;
  G4double hpXscv(0);
  G4double hnXscv(0);

  G4int Nt = At-Zt;              // number of neutrons
  if (Nt < 0) Nt = 0;  

  G4double aa = At;
  G4double zz = Zt;
  G4double nn = Nt;

  G4double targ_mass = G4ParticleTable::GetParticleTable()->
  GetIonTable()->GetIonMass(Zt, At);

  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???

  G4double proj_mass     = aParticle->GetMass();
  G4double proj_energy   = aParticle->GetTotalEnergy(); 
  G4double proj_momentum = aParticle->GetMomentum().mag();

  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );

  sMand         /= GeV*GeV;  // in GeV for parametrisation
  proj_momentum /= GeV;
  proj_energy   /= GeV;
  proj_mass     /= GeV;

  // General PDG fit constants

  G4double s0   = 5.38*5.38; // in Gev^2
  G4double eta1 = 0.458;
  G4double eta2 = 0.458;
  G4double B    = 0.308;


  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
  

  if(theParticle == theNeutron) 
  {
    if( proj_momentum >= 373.)
    {
      return GetHadronNucleonXscPDG(aParticle,At,Zt);
    }
    else if( proj_momentum >= 10.)
    // if( proj_momentum >= 2.)
    {
      //  Delta = 1.;  // DHW 19 May 2011: variable set but not used
      // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;

      if(proj_momentum >= 10.)
      {
        B0 = 7.5;
        A0 = 100. - B0*std::log(3.0e7);

        xsection = A0 + B0*std::log(proj_energy) - 11
                  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
                     0.93827*0.93827,-0.165);        //  mb
      }
      xsection *= zz + nn;
    }
    else
    {
      // nn to be pp

      if( proj_momentum < 0.73 )
      {
        hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
      }
      else if( proj_momentum < 1.05  )
      {
       hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
                         (std::log(proj_momentum/0.73));
      }
      else  // if( proj_momentum < 10.  )
      {
         hnXscv = 39.0+
              75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
      }
      // pn to be np

      if( proj_momentum < 0.8 )
      {
        hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
      }      
      else if( proj_momentum < 1.4 )
      {
        hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
      }
      else    // if( proj_momentum < 10.  )
      {
        hpXscv = 33.3+
              20.8*(std::pow(proj_momentum,2.0)-1.35)/
                 (std::pow(proj_momentum,2.50)+0.95);
      }
      xsection = hpXscv*zz + hnXscv*nn;
    }
  } 
  else if(theParticle == theProton) 
  {
    if( proj_momentum >= 373.)
    {
      return GetHadronNucleonXscPDG(aParticle,At,Zt);
    }
    else if( proj_momentum >= 10.)
    // if( proj_momentum >= 2.)
    {
      // Delta = 1.;  DHW 19 May 2011: variable set but not used
      // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;

      if(proj_momentum >= 10.)
      {
        B0 = 7.5;
        A0 = 100. - B0*std::log(3.0e7);

        xsection = A0 + B0*std::log(proj_energy) - 11
                  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
                     0.93827*0.93827,-0.165);        //  mb
      }
      xsection *= zz + nn;
    }
    else
    {
      // pp

      if( proj_momentum < 0.73 )
      {
        hpXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
      }
      else if( proj_momentum < 1.05  )
      {
       hpXscv = 23 + 40*(std::log(proj_momentum/0.73))*
                         (std::log(proj_momentum/0.73));
      }
      else    // if( proj_momentum < 10.  )
      {
         hpXscv = 39.0+
              75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
      }
      // pn to be np

      if( proj_momentum < 0.8 )
      {
        hnXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
      }      
      else if( proj_momentum < 1.4 )
      {
        hnXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
      }
      else   // if( proj_momentum < 10.  )
      {
        hnXscv = 33.3+
              20.8*(std::pow(proj_momentum,2.0)-1.35)/
                 (std::pow(proj_momentum,2.50)+0.95);
      }
      xsection = hpXscv*zz + hnXscv*nn;
      // xsection = hpXscv*(Zt + Nt);
      // xsection = hnXscv*(Zt + Nt);
    }    
    // xsection *= 0.95;
  } 
  else if( theParticle == theAProton ) 
  {
    // xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
    //                       + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));

    // xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
    //                    + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));

    G4double logP = std::log(proj_momentum);

    if( proj_momentum <= 1.0 )
    {
      xsection  = zz*(65.55 + 53.84/(proj_momentum+1.e-6)  );
    }
    else
    {
      xsection  = zz*( 41.1 + 77.2*std::pow( proj_momentum, -0.68) 
                       + 0.293*logP*logP - 1.82*logP );
    }
    if ( nn > 0.)  
    {
      xsection += nn*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP);
    }
    else // H
    {
      fInelasticXsc =   38.0 + 38.0*std::pow( proj_momentum, -0.96) 
                        - 0.169*logP*logP;
      fInelasticXsc *=  millibarn;
    }    
  } 
  else if( theParticle == thePiPlus ) 
  {
    if(proj_momentum < 0.4)
    {
      G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
      hpXscv      = Ex3+20.0;
    }
    else if( proj_momentum < 1.15 )
    {
      G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
      hpXscv = Ex4+14.0;
    }
    else if(proj_momentum < 3.5)
    {
      G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
      G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
      hpXscv = Ex1+Ex2+27.5;
    }
    else //  if(proj_momentum > 3.5) // mb
    {
      hpXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
    }
    // pi+n = pi-p??

    if(proj_momentum < 0.37)
    {
      hnXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
    }
    else if(proj_momentum<0.65)
    {
       hnXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
    }
    else if(proj_momentum<1.3)
    {
      hnXscv = 36.1+
                10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
                24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
    }
    else if(proj_momentum<3.0)
    {
      hnXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
                3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
                1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
    }
    else   // mb
    {
      hnXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); 
    }
    xsection = hpXscv*zz + hnXscv*nn;
  } 
  else if(theParticle == thePiMinus) 
  {
    // pi-n = pi+p??

    if(proj_momentum < 0.4)
    {
      G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
      hnXscv      = Ex3+20.0;
    }
    else if(proj_momentum < 1.15)
    {
      G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
      hnXscv = Ex4+14.0;
    }
    else if(proj_momentum < 3.5)
    {
      G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
      G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
      hnXscv = Ex1+Ex2+27.5;
    }
    else //  if(proj_momentum > 3.5) // mb
    {
      hnXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
    }
    // pi-p

    if(proj_momentum < 0.37)
    {
      hpXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
    }
    else if(proj_momentum<0.65)
    {
       hpXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
    }
    else if(proj_momentum<1.3)
    {
      hpXscv = 36.1+
                10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
                24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
    }
    else if(proj_momentum<3.0)
    {
      hpXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
                3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
                1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
    }
    else   // mb
    {
      hpXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); 
    }
    xsection = hpXscv*zz + hnXscv*nn;
  } 
  else if(theParticle == theKPlus) 
  {
    xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
                          + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));

    xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
                          + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
  } 
  else if(theParticle == theKMinus) 
  {
    xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
                          + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));

    xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
                          + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
  }
  else if(theParticle == theSMinus) 
  {
    xsection  = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.) 
                          - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
  } 
  else if(theParticle == theGamma) // modify later on
  {
    xsection  = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.) 
                          + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
   
  } 
  else  // as proton ??? 
  {
    xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
                          + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));

    xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
  } 
  xsection *= millibarn; // parametrised in mb
  return xsection;
}

G4double G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS	(	const G4DynamicParticle *	,
		const G4Element *
	)

Definition at line 664 of file G4ComponentGGHadronNucleusXsc.cc.

References G4lrint(), G4Element::GetN(), and G4Element::GetZ().

Referenced by GetHNinelasticXsc(), GetIsoCrossSection(), GetRatioQE(), and GetRatioSD().

{
  G4int At = G4lrint(anElement->GetN());  // number of nucleons 
  G4int Zt = G4lrint(anElement->GetZ());  // number of protons

  return GetHadronNucleonXscNS(aParticle, At, Zt);
}

G4double G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG	(	const G4DynamicParticle *	,
		G4int	At,
		G4int	Zt
	)

Definition at line 549 of file G4ComponentGGHadronNucleusXsc.cc.

References CalcMandelstamS(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), G4ParticleTable::GetParticleTable(), G4InuclParticleNames::nn, and G4InuclParticleNames::s0.

{
  G4double xsection;

  G4int Nt = At-Zt;              // number of neutrons
  if (Nt < 0) Nt = 0;
  
  G4double zz = Zt;
  G4double aa = At;
  G4double nn = Nt;

  G4double targ_mass = G4ParticleTable::GetParticleTable()->
    GetIonTable()->GetIonMass(Zt, At);

  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???

  G4double proj_mass     = aParticle->GetMass(); 
  G4double proj_momentum = aParticle->GetMomentum().mag();

  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );

  sMand         /= GeV*GeV;  // in GeV for parametrisation

  // General PDG fit constants

  G4double s0   = 5.38*5.38; // in Gev^2
  G4double eta1 = 0.458;
  G4double eta2 = 0.458;
  G4double B    = 0.308;


  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
  

  if(theParticle == theNeutron) // proton-neutron fit 
  {
    xsection = zz*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    xsection  += nn*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
                      + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
  } 
  else if(theParticle == theProton) 
  {
      
      xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
                          + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));

      xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
  } 
  else if(theParticle == theAProton) 
  {
    xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
                          + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));

    xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
                          + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
  } 
  else if(theParticle == thePiPlus) 
  {
    xsection  = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.) 
                          + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
  } 
  else if(theParticle == thePiMinus) 
  {
    xsection  = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.) 
                          + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
  } 
  else if(theParticle == theKPlus || theParticle == theK0L ) 
  {
    xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
                          + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));

    xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
                          + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
  } 
  else if(theParticle == theKMinus || theParticle == theK0S ) 
  {
    xsection  = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 
                          + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));

    xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 
                          + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
  }
  else if(theParticle == theSMinus) 
  {
    xsection  = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.) 
                          - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
  } 
  else if(theParticle == theGamma) // modify later on
  {
    xsection  = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.) 
                          + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
   
  } 
  else  // as proton ??? 
  {
    xsection  = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
                          + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));

    xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
  } 
  xsection *= millibarn; // parametrised in mb
  return xsection;
}

G4double G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG	(	const G4DynamicParticle *	,
		const G4Element *
	)

Definition at line 530 of file G4ComponentGGHadronNucleusXsc.cc.

References G4lrint(), G4Element::GetN(), and G4Element::GetZ().

Referenced by GetHadronNucleonXscNS(), and GetKaonNucleonXscVector().

{
  G4int At = G4lrint(anElement->GetN());  // number of nucleons 
  G4int Zt = G4lrint(anElement->GetZ());  // number of protons

  return GetHadronNucleonXscPDG(aParticle, At, Zt);
}

G4double G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc	(	const G4DynamicParticle *	,
		G4int	At,
		G4int	Zt
	)

Definition at line 1072 of file G4ComponentGGHadronNucleusXsc.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), and GetHNinelasticXscVU().

{
  G4ParticleDefinition* hadron = aParticle->GetDefinition();
  G4double sumInelastic;
  G4int Nt = At - Zt;
  if(Nt < 0) Nt = 0;
  
  if( hadron == theKPlus )
  {
    sumInelastic =  GetHNinelasticXscVU(aParticle, At, Zt);
  }
  else
  {
    //sumInelastic  = Zt*GetHadronNucleonXscMK(aParticle, theProton);
    // sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);    
    sumInelastic  = G4double(Zt)*GetHadronNucleonXscNS(aParticle, 1, 1);
    sumInelastic += G4double(Nt)*GetHadronNucleonXscNS(aParticle, 1, 0);    
  } 
  return sumInelastic;
}

G4double G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc	(	const G4DynamicParticle *	,
		const G4Element *
	)

Definition at line 1058 of file G4ComponentGGHadronNucleusXsc.cc.

References G4lrint(), G4Element::GetN(), and G4Element::GetZ().

Referenced by GetIsoCrossSection(), and GetRatioQE().

{
  G4int At = G4lrint(anElement->GetN());  // number of nucleons 
  G4int Zt = G4lrint(anElement->GetZ());  // number of protons

  return GetHNinelasticXsc(aParticle, At, Zt);
}

G4double G4ComponentGGHadronNucleusXsc::GetHNinelasticXscVU	(	const G4DynamicParticle *	,
		G4int	At,
		G4int	Zt
	)

Definition at line 1100 of file G4ComponentGGHadronNucleusXsc.cc.

References G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetMomentum(), G4ParticleDefinition::GetPDGEncoding(), and G4DynamicParticle::GetTotalEnergy().

Referenced by GetHNinelasticXsc().

{
  G4int PDGcode    = aParticle->GetDefinition()->GetPDGEncoding();
  G4int absPDGcode = std::abs(PDGcode);

  G4double Elab = aParticle->GetTotalEnergy();              
                          // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
  G4double Plab = aParticle->GetMomentum().mag();            
                          // std::sqrt(Elab * Elab - 0.88);

  Elab /= GeV;
  Plab /= GeV;

  G4double LogPlab    = std::log( Plab );
  G4double sqrLogPlab = LogPlab * LogPlab;

  //G4cout<<"Plab = "<<Plab<<G4endl;

  G4double NumberOfTargetProtons = G4double(Zt); 
  G4double NumberOfTargetNucleons = G4double(At);
  G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;

  if(NumberOfTargetNeutrons < 0.0) NumberOfTargetNeutrons = 0.0;

  G4double Xtotal, Xelastic, Xinelastic;

  if( absPDGcode > 1000 )  //------Projectile is baryon --------
  {
       G4double XtotPP = 48.0 +  0. *std::pow(Plab, 0.  ) +
                         0.522*sqrLogPlab - 4.51*LogPlab;

       G4double XtotPN = 47.3 +  0. *std::pow(Plab, 0.  ) +
                         0.513*sqrLogPlab - 4.27*LogPlab;

       G4double XelPP  = 11.9 + 26.9*std::pow(Plab,-1.21) +
                         0.169*sqrLogPlab - 1.85*LogPlab;

       G4double XelPN  = 11.9 + 26.9*std::pow(Plab,-1.21) +
                         0.169*sqrLogPlab - 1.85*LogPlab;

       Xtotal          = (NumberOfTargetProtons * XtotPP +
                          NumberOfTargetNeutrons * XtotPN);

       Xelastic        = (NumberOfTargetProtons * XelPP +
                          NumberOfTargetNeutrons * XelPN);
  }
  else if( PDGcode ==  211 ) //------Projectile is PionPlus -------
  {
       G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) +
                          0.19 *sqrLogPlab - 0.0 *LogPlab;

       G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) +
                          0.456*sqrLogPlab - 4.03*LogPlab;

       G4double XelPiP  =  0.0 + 11.4*std::pow(Plab,-0.40) +
                           0.079*sqrLogPlab - 0.0 *LogPlab;

       G4double XelPiN  = 1.76 + 11.2*std::pow(Plab,-0.64) +
                          0.043*sqrLogPlab - 0.0 *LogPlab;

       Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
                            NumberOfTargetNeutrons * XtotPiN  );

       Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
                            NumberOfTargetNeutrons * XelPiN   );
  }
  else if( PDGcode == -211 ) //------Projectile is PionMinus -------
  {
       G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) +
                          0.456*sqrLogPlab - 4.03*LogPlab;

       G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) +
                          0.19 *sqrLogPlab - 0.0 *LogPlab;

       G4double XelPiP  = 1.76 + 11.2*std::pow(Plab,-0.64) +
                          0.043*sqrLogPlab - 0.0 *LogPlab;

       G4double XelPiN  =  0.0 + 11.4*std::pow(Plab,-0.40) +
                           0.079*sqrLogPlab - 0.0 *LogPlab;

       Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
                            NumberOfTargetNeutrons * XtotPiN  );

       Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
                            NumberOfTargetNeutrons * XelPiN   );
  }
  else if( PDGcode ==  111 )  //------Projectile is PionZero  -------
  {
       G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) +
                          0.19 *sqrLogPlab - 0.0 *LogPlab +   //Pi+
                          33.0 + 14.0 *std::pow(Plab,-1.36) +
                          0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-

       G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) +
                          0.456*sqrLogPlab - 4.03*LogPlab +   //Pi+
                          16.4 + 19.3 *std::pow(Plab,-0.42) +
                          0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-

       G4double XelPiP  =( 0.0 + 11.4*std::pow(Plab,-0.40) +
                           0.079*sqrLogPlab - 0.0 *LogPlab +    //Pi+
                           1.76 + 11.2*std::pow(Plab,-0.64) +
                           0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-

       G4double XelPiN  =( 1.76 + 11.2*std::pow(Plab,-0.64) +
                           0.043*sqrLogPlab - 0.0 *LogPlab +   //Pi+
                           0.0  + 11.4*std::pow(Plab,-0.40) +
                           0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-

       Xtotal           = ( NumberOfTargetProtons  * XtotPiP +
                            NumberOfTargetNeutrons * XtotPiN  );

       Xelastic         = ( NumberOfTargetProtons  * XelPiP  +
                            NumberOfTargetNeutrons * XelPiN   );
  }
  else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
  {
       G4double XtotKP = 18.1 +  0. *std::pow(Plab, 0.  ) +
                         0.26 *sqrLogPlab - 1.0 *LogPlab;
       G4double XtotKN = 18.7 +  0. *std::pow(Plab, 0.  ) +
                         0.21 *sqrLogPlab - 0.89*LogPlab;

       G4double XelKP  =  5.0 +  8.1*std::pow(Plab,-1.8 ) +
                          0.16 *sqrLogPlab - 1.3 *LogPlab;

       G4double XelKN  =  7.3 +  0. *std::pow(Plab,-0.  ) +
                          0.29 *sqrLogPlab - 2.4 *LogPlab;

       Xtotal          = ( NumberOfTargetProtons  * XtotKP +
                           NumberOfTargetNeutrons * XtotKN  );

       Xelastic        = ( NumberOfTargetProtons  * XelKP  +
                           NumberOfTargetNeutrons * XelKN   );
  }
  else if( PDGcode ==-321 )  //------Projectile is KaonMinus ------
  {
       G4double XtotKP = 32.1 +  0. *std::pow(Plab, 0.  ) +
                         0.66 *sqrLogPlab - 5.6 *LogPlab;
       G4double XtotKN = 25.2 +  0. *std::pow(Plab, 0.  ) +
                         0.38 *sqrLogPlab - 2.9 *LogPlab;

       G4double XelKP  =  7.3 +  0. *std::pow(Plab,-0.  ) +
                          0.29 *sqrLogPlab - 2.4 *LogPlab;

       G4double XelKN  =  5.0 +  8.1*std::pow(Plab,-1.8 ) +
                          0.16 *sqrLogPlab - 1.3 *LogPlab;

       Xtotal          = ( NumberOfTargetProtons  * XtotKP +
                           NumberOfTargetNeutrons * XtotKN  );

       Xelastic        = ( NumberOfTargetProtons  * XelKP  +
                           NumberOfTargetNeutrons * XelKN   );
  }
  else if( PDGcode == 311 ) //------Projectile is KaonZero ------
  {
       G4double XtotKP = ( 18.1 +  0. *std::pow(Plab, 0.  ) +
                          0.26 *sqrLogPlab - 1.0 *LogPlab +   //K+
                          32.1 +  0. *std::pow(Plab, 0.  ) +
                          0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-

       G4double XtotKN = ( 18.7 +  0. *std::pow(Plab, 0.  ) +
                          0.21 *sqrLogPlab - 0.89*LogPlab +   //K+
                          25.2 +  0. *std::pow(Plab, 0.  ) +
                          0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-

       G4double XelKP  = (  5.0 +  8.1*std::pow(Plab,-1.8 )
                           + 0.16 *sqrLogPlab - 1.3 *LogPlab +   //K+
                           7.3 +  0. *std::pow(Plab,-0.  ) +
                           0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-

       G4double XelKN  = (  7.3 +  0. *std::pow(Plab,-0.  ) +
                           0.29 *sqrLogPlab - 2.4 *LogPlab +   //K+
                           5.0 +  8.1*std::pow(Plab,-1.8 ) +
                           0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-

       Xtotal          = ( NumberOfTargetProtons  * XtotKP +
                           NumberOfTargetNeutrons * XtotKN  );

       Xelastic        = ( NumberOfTargetProtons  * XelKP  +
                           NumberOfTargetNeutrons * XelKN   );
  }
  else  //------Projectile is undefined, Nucleon assumed
  {
       G4double XtotPP = 48.0 +  0. *std::pow(Plab, 0.  ) +
                         0.522*sqrLogPlab - 4.51*LogPlab;

       G4double XtotPN = 47.3 +  0. *std::pow(Plab, 0.  ) +
                         0.513*sqrLogPlab - 4.27*LogPlab;

       G4double XelPP  = 11.9 + 26.9*std::pow(Plab,-1.21) +
                         0.169*sqrLogPlab - 1.85*LogPlab;
       G4double XelPN  = 11.9 + 26.9*std::pow(Plab,-1.21) +
                         0.169*sqrLogPlab - 1.85*LogPlab;

       Xtotal          = ( NumberOfTargetProtons  * XtotPP +
                           NumberOfTargetNeutrons * XtotPN  );

       Xelastic        = ( NumberOfTargetProtons  * XelPP  +
                           NumberOfTargetNeutrons * XelPN   );
  }
  Xinelastic = Xtotal - Xelastic;

  if( Xinelastic < 0.) Xinelastic = 0.;

  return Xinelastic*= millibarn;
}

G4double G4ComponentGGHadronNucleusXsc::GetInelasticElementCrossSection	(	const G4ParticleDefinition *	aParticle,
		G4double	kinEnergy,
		G4int	Z,
		G4double	A
	)			[virtual]

Implements G4VComponentCrossSection.

Definition at line 138 of file G4ComponentGGHadronNucleusXsc.cc.

References GetIsoCrossSection().

{
  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.), 
                                                kinEnergy);
  fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
  delete aDP;

  return fInelasticXsc;
}

G4double G4ComponentGGHadronNucleusXsc::GetInelasticGlauberGribov	(	const G4DynamicParticle *	,
		G4int	Z,
		G4int	A
	)			[inline]

Definition at line 222 of file G4ComponentGGHadronNucleusXsc.hh.

References GetIsoCrossSection().

{
  GetIsoCrossSection(dp, Z, A);
  return fInelasticXsc;
}

G4double G4ComponentGGHadronNucleusXsc::GetInelasticGlauberGribovXsc

(

)

[inline]

Definition at line 139 of file G4ComponentGGHadronNucleusXsc.hh.

{ return fInelasticXsc; }; 

G4double G4ComponentGGHadronNucleusXsc::GetInelasticIsotopeCrossSection	(	const G4ParticleDefinition *	aParticle,
		G4double	kinEnergy,
		G4int	Z,
		G4int	A
	)			[virtual]

Implements G4VComponentCrossSection.

Definition at line 124 of file G4ComponentGGHadronNucleusXsc.cc.

References GetIsoCrossSection().

{
  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.), 
                                                kinEnergy);
  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
  delete aDP;

  return fInelasticXsc;
}

G4double G4ComponentGGHadronNucleusXsc::GetIsoCrossSection	(	const G4DynamicParticle *	,
		G4int	Z,
		G4int	A,
		const G4Isotope *	iso = 0,
		const G4Element *	elm = 0,
		const G4Material *	mat = 0
	)

Definition at line 240 of file G4ComponentGGHadronNucleusXsc.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), G4HadronNucleonXsc::GetHadronNucleonXscNS(), GetHNinelasticXsc(), G4HadronNucleonXsc::GetInelasticHadronNucleonXsc(), GetKaonNucleonXscVector(), GetNucleusRadius(), GetParticleBarCorIn(), GetParticleBarCorTot(), and G4INCL::Math::pi.

Referenced by ComputeQuasiElasticRatio(), GetElasticElementCrossSection(), GetElasticGlauberGribov(), GetElasticIsotopeCrossSection(), GetInelasticElementCrossSection(), GetInelasticGlauberGribov(), GetInelasticIsotopeCrossSection(), GetTotalElementCrossSection(), and GetTotalIsotopeCrossSection().

{
  G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
  G4double hpInXsc(0.), hnInXsc(0.);
  G4double R             = GetNucleusRadius(A); 

  G4int N = A - Z;              // number of neutrons
  if (N < 0) N = 0;

  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();

  if( theParticle == theProton  || 
      theParticle == theNeutron ||
      theParticle == thePiPlus  || 
      theParticle == thePiMinus      )
  {
    // sigma        = GetHadronNucleonXscNS(aParticle, A, Z);

    sigma = Z*hnXsc->GetHadronNucleonXscNS(aParticle, theProton);

    hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();

    sigma += N*hnXsc->GetHadronNucleonXscNS(aParticle, theNeutron);

    hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();

    cofInelastic = 2.4;
    cofTotal     = 2.0;
  }
  else if( theParticle == theKPlus   || 
           theParticle == theKMinus  || 
           theParticle == theK0S     || 
           theParticle == theK0L        ) 
  {
    sigma        = GetKaonNucleonXscVector(aParticle, A, Z);
    cofInelastic = 2.2;
    cofTotal     = 2.0;
    R = 1.3*fermi;
    R *= std::pow(G4double(A), 0.3333);
  }
  else
  {
    sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
    cofInelastic = 2.2;
    cofTotal     = 2.0;
  }
  // cofInelastic = 2.0;

  if( A > 1 )
  { 
    nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
    ratio = sigma/nucleusSquare;

    xsection =  nucleusSquare*std::log( 1. + ratio );

    xsection *= GetParticleBarCorTot(theParticle, Z);

    fTotalXsc = xsection;

  

    fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;

    fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);

    fElasticXsc   = fTotalXsc - fInelasticXsc;

    if(fElasticXsc < 0.) fElasticXsc = 0.;
    
    G4double difratio = ratio/(1.+ratio);

    fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );


    // sigma = GetHNinelasticXsc(aParticle, A, Z);

    sigma = Z*hpInXsc + N*hnInXsc;

    ratio = sigma/nucleusSquare;

    fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;

    if (fElasticXsc < 0.) fElasticXsc = 0.;
  }
  else // H
  {
    fTotalXsc = sigma;
    xsection  = sigma;
    
    if ( theParticle != theAProton ) 
    {
      sigma         = GetHNinelasticXsc(aParticle, A, Z);
      fInelasticXsc = sigma;
      fElasticXsc   = fTotalXsc - fInelasticXsc;      
    }
    else
    {
      fElasticXsc   = fTotalXsc - fInelasticXsc;
    }
    if (fElasticXsc < 0.) fElasticXsc = 0.;
      
  }
  return xsection; 
}

G4double G4ComponentGGHadronNucleusXsc::GetKaonNucleonXscVector	(	const G4DynamicParticle *	,
		G4int	At,
		G4int	Zt
	)

Definition at line 1020 of file G4ComponentGGHadronNucleusXsc.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscPDG(), G4DynamicParticle::GetKineticEnergy(), G4HadronNucleonXsc::GetKmNeutronTotXscVector(), G4HadronNucleonXsc::GetKmProtonTotXscVector(), G4HadronNucleonXsc::GetKpNeutronTotXscVector(), and G4HadronNucleonXsc::GetKpProtonTotXscVector().

Referenced by GetIsoCrossSection().

{
  G4double Tkin, logTkin, xsc, xscP, xscN;
  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();

  G4int Nt = At-Zt;              // number of neutrons
  if (Nt < 0) Nt = 0;  

  Tkin = aParticle->GetKineticEnergy(); // Tkin in MeV

  if( Tkin > 70*GeV ) return GetHadronNucleonXscPDG(aParticle,At,Zt);

  logTkin = std::log(Tkin); // Tkin in MeV!!!

 if( theParticle == theKPlus )
 {
   xscP = hnXsc->GetKpProtonTotXscVector(logTkin);
   xscN = hnXsc->GetKpNeutronTotXscVector(logTkin);
 }
 else if( theParticle == theKMinus )
 {
   xscP = hnXsc->GetKmProtonTotXscVector(logTkin);
   xscN = hnXsc->GetKmNeutronTotXscVector(logTkin);
 }
 else // K-zero as half of K+ and K-
 {
   xscP = (hnXsc->GetKpProtonTotXscVector(logTkin)+hnXsc->GetKmProtonTotXscVector(logTkin))*0.5;
   xscN = (hnXsc->GetKpNeutronTotXscVector(logTkin)+hnXsc->GetKmNeutronTotXscVector(logTkin))*0.5;
 }
 xsc = xscP*Zt + xscN*Nt;
  return xsc;
}

G4double G4ComponentGGHadronNucleusXsc::GetNucleusRadius

(

G4int

At

)

Definition at line 1369 of file G4ComponentGGHadronNucleusXsc.cc.

References G4INCL::Math::oneThird.

{
  G4double oneThird = 1.0/3.0;
  G4double cubicrAt = std::pow(G4double(At), oneThird); 

  G4double R;  // = fRadiusConst*cubicrAt;

  /*
  G4double tmp = std::pow( cubicrAt-1., 3.);
  tmp         += At;
  tmp         *= 0.5;

  if (At > 20.)
  {
    R = fRadiusConst*std::pow (tmp, oneThird); 
  }
  else
  {
    R = fRadiusConst*cubicrAt; 
  }
  */

  R = fRadiusConst*cubicrAt;

  G4double meanA = 20.;
  G4double tauA  = 20.; 

  if (At > 20)   // 20.
  {
    R *= ( 0.8 + 0.2*std::exp( -(G4double(At) - meanA)/tauA) ); 
  }
  else
  {
    R *= ( 1.0 + 0.1*( 1. - std::exp( (G4double(At) - meanA)/tauA) ) ); 
  }

  return R;
}

G4double G4ComponentGGHadronNucleusXsc::GetNucleusRadius	(	const G4DynamicParticle *	,
		const G4Element *
	)

Definition at line 1312 of file G4ComponentGGHadronNucleusXsc.cc.

References G4lrint(), G4Element::GetN(), and G4INCL::Math::oneThird.

Referenced by GetIsoCrossSection(), GetRatioQE(), and GetRatioSD().

{
  G4int At = G4lrint(anElement->GetN());
  G4double oneThird = 1.0/3.0;
  G4double cubicrAt = std::pow(G4double(At), oneThird); 

  G4double R;  // = fRadiusConst*cubicrAt;
  /*  
  G4double tmp = std::pow( cubicrAt-1., 3.);
  tmp         += At;
  tmp         *= 0.5;

  if (At > 20.)   // 20.
  {
    R = fRadiusConst*std::pow (tmp, oneThird); 
  }
  else
  {
    R = fRadiusConst*cubicrAt; 
  }
  */
  
  R = fRadiusConst*cubicrAt;

  G4double meanA  = 21.;

  G4double tauA1  = 40.; 
  G4double tauA2  = 10.; 
  G4double tauA3  = 5.; 

  G4double a1 = 0.85;
  G4double b1 = 1. - a1;

  G4double b2 = 0.3;
  G4double b3 = 4.;

  if (At > 20)   // 20.
  {
    R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) ); 
  }
  else if (At > 3)
  {
    R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) ); 
  }
  else 
  {
    R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) ); 
  }  
  return R;
 
}

G4double G4ComponentGGHadronNucleusXsc::GetParticleBarCorIn	(	const G4ParticleDefinition *	theParticle,
		G4int	Z
	)			[inline]

Definition at line 255 of file G4ComponentGGHadronNucleusXsc.hh.

Referenced by GetIsoCrossSection().

{
  if(Z >= 2 && Z <= 92)
  {
    if(      theParticle == theProton ) return fProtonBarCorrectionIn[Z]; 
    else if( theParticle == theNeutron) return fNeutronBarCorrectionIn[Z]; 
    else if( theParticle == thePiPlus ) return fPionPlusBarCorrectionIn[Z];
    else if( theParticle == thePiMinus) return fPionMinusBarCorrectionIn[Z];
    else return 1.0;
  }
  else return 1.0;
}

G4double G4ComponentGGHadronNucleusXsc::GetParticleBarCorTot	(	const G4ParticleDefinition *	theParticle,
		G4int	Z
	)			[inline]

Definition at line 235 of file G4ComponentGGHadronNucleusXsc.hh.

Referenced by GetIsoCrossSection().

{
  if(Z >= 2 && Z <= 92)
  {
    if(      theParticle == theProton ) return fProtonBarCorrectionTot[Z]; 
    else if( theParticle == theNeutron) return fNeutronBarCorrectionTot[Z]; 
    else if( theParticle == thePiPlus ) return fPionPlusBarCorrectionTot[Z];
    else if( theParticle == thePiMinus) return fPionMinusBarCorrectionTot[Z];
    else return 1.0;
  }
  else return 1.0;
}

G4double G4ComponentGGHadronNucleusXsc::GetProductionGlauberGribovXsc

(

)

[inline]

Definition at line 140 of file G4ComponentGGHadronNucleusXsc.hh.

{ return fProductionXsc; }; 

G4double G4ComponentGGHadronNucleusXsc::GetRadiusConst

(

)

[inline]

Definition at line 142 of file G4ComponentGGHadronNucleusXsc.hh.

{ return fRadiusConst;  }; 

G4double G4ComponentGGHadronNucleusXsc::GetRatioQE	(	const G4DynamicParticle *	,
		G4int	At,
		G4int	Zt
	)

Definition at line 396 of file G4ComponentGGHadronNucleusXsc.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), GetHNinelasticXsc(), GetNucleusRadius(), and G4INCL::Math::pi.

{
  G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
  G4double R             = GetNucleusRadius(A); 

  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();

  if( theParticle == theProton  || 
      theParticle == theNeutron ||
      theParticle == thePiPlus  || 
      theParticle == thePiMinus      )
  {
    sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
    cofInelastic = 2.4;
    cofTotal     = 2.0;
  }
  else
  {
    sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
    cofInelastic = 2.2;
    cofTotal     = 2.0;
  }
  nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
  ratio = sigma/nucleusSquare;

  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;

  sigma = GetHNinelasticXsc(aParticle, A, Z);
  ratio = sigma/nucleusSquare;

  fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;

  if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
  else                                ratio = 0.;
  if ( ratio < 0. )                   ratio = 0.;

  return ratio; 
}

G4double G4ComponentGGHadronNucleusXsc::GetRatioSD	(	const G4DynamicParticle *	,
		G4int	At,
		G4int	Zt
	)

Definition at line 354 of file G4ComponentGGHadronNucleusXsc.cc.

References G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), GetNucleusRadius(), and G4INCL::Math::pi.

{
  G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
  G4double R             = GetNucleusRadius(A); 

  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();

  if( theParticle == theProton  || 
      theParticle == theNeutron ||
      theParticle == thePiPlus  || 
      theParticle == thePiMinus      )
  {
    sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
    cofInelastic = 2.4;
    cofTotal     = 2.0;
  }
  else
  {
    sigma        = GetHadronNucleonXscNS(aParticle, A, Z);
    cofInelastic = 2.2;
    cofTotal     = 2.0;
  }
  nucleusSquare = cofTotal*pi*R*R;   // basically 2piRR
  ratio = sigma/nucleusSquare;

  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
   
  G4double difratio = ratio/(1.+ratio);

  fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );

  if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
  else                    ratio = 0.;

  return ratio; 
}

G4double G4ComponentGGHadronNucleusXsc::GetTotalElementCrossSection	(	const G4ParticleDefinition *	aParticle,
		G4double	kinEnergy,
		G4int	Z,
		G4double	A
	)			[virtual]

Implements G4VComponentCrossSection.

Definition at line 110 of file G4ComponentGGHadronNucleusXsc.cc.

References GetIsoCrossSection().

{
  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.), 
                                                kinEnergy);
  fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
  delete aDP;

  return fTotalXsc;
}

G4double G4ComponentGGHadronNucleusXsc::GetTotalGlauberGribovXsc

(

)

[inline]

Definition at line 137 of file G4ComponentGGHadronNucleusXsc.hh.

{ return fTotalXsc;     }; 

G4double G4ComponentGGHadronNucleusXsc::GetTotalIsotopeCrossSection	(	const G4ParticleDefinition *	aParticle,
		G4double	kinEnergy,
		G4int	Z,
		G4int	A
	)			[virtual]

Implements G4VComponentCrossSection.

Definition at line 96 of file G4ComponentGGHadronNucleusXsc.cc.

References GetIsoCrossSection().

{
  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.), 
                                                kinEnergy);
  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
  delete aDP;

  return fTotalXsc;
}

G4bool G4ComponentGGHadronNucleusXsc::IsIsoApplicable	(	const G4DynamicParticle *	aDP,
		G4int	Z,
		G4int	A,
		const G4Element *	elm = 0,
		const G4Material *	mat = 0
	)

Definition at line 204 of file G4ComponentGGHadronNucleusXsc.cc.

References G4DynamicParticle::GetDefinition(), and G4DynamicParticle::GetKineticEnergy().

{
  G4bool applicable      = false;
  // G4int baryonNumber     = aDP->GetDefinition()->GetBaryonNumber();
  G4double kineticEnergy = aDP->GetKineticEnergy();

  const G4ParticleDefinition* theParticle = aDP->GetDefinition();
 
  if ( ( kineticEnergy  >= fLowerLimit &&
         Z > 1 &&      // >=  He
       ( theParticle == theAProton   ||
         theParticle == theGamma     ||
         theParticle == theKPlus     ||
         theParticle == theKMinus    || 
         theParticle == theK0L     ||
         theParticle == theK0S    || 
         theParticle == theSMinus    ||  
         theParticle == theProton    ||
         theParticle == theNeutron   ||   
         theParticle == thePiPlus    ||
         theParticle == thePiMinus       ) )    ) applicable = true;

  return applicable;
}

void G4ComponentGGHadronNucleusXsc::SetEnergyLowerLimit

(

G4double

E

)

[inline]

Definition at line 147 of file G4ComponentGGHadronNucleusXsc.hh.

{fLowerLimit=E;};

---

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

• G4ComponentGGHadronNucleusXsc.hh
• G4ComponentGGHadronNucleusXsc.cc

---