#include <G4NeutronElectronElModel.hh>

Inheritance diagram for G4NeutronElectronElModel:

Public Member Functions
	G4NeutronElectronElModel (const G4String &name="n-e-elastic")

virtual	~G4NeutronElectronElModel ()

virtual G4bool	IsApplicable (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)

virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)

void	Initialise ()

G4double	XscIntegrand (G4double x)

G4double	CalculateAm (G4double momentum)

G4double	GetAm ()

G4double	SampleSin2HalfTheta (G4double Tkin)

G4double	GetTransfer (G4int iTkin, G4int iTransfer, G4double position)

void	SetCutEnergy (G4double ec)

G4double	GetCutEnergy ()

virtual void	ModelDescription (std::ostream &) const

Public Member Functions inherited from G4HadronElastic
	G4HadronElastic (const G4String &name="hElasticLHEP")

	~G4HadronElastic () override

G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override

G4double	GetSlopeCof (const G4int pdg)

void	SetLowestEnergyLimit (G4double value)

G4double	LowestEnergyLimit () const

G4double	ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)

Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")

virtual	~G4HadronicInteraction ()

G4double	GetMinEnergy () const

G4double	GetMinEnergy (const G4Material aMaterial, const G4Element anElement) const

void	SetMinEnergy (G4double anEnergy)

void	SetMinEnergy (G4double anEnergy, const G4Element *anElement)

void	SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)

G4double	GetMaxEnergy () const

G4double	GetMaxEnergy (const G4Material aMaterial, const G4Element anElement) const

void	SetMaxEnergy (const G4double anEnergy)

void	SetMaxEnergy (G4double anEnergy, const G4Element *anElement)

void	SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)

G4int	GetVerboseLevel () const

void	SetVerboseLevel (G4int value)

const G4String &	GetModelName () const

void	DeActivateFor (const G4Material *aMaterial)

void	ActivateFor (const G4Material *aMaterial)

void	DeActivateFor (const G4Element *anElement)

void	ActivateFor (const G4Element *anElement)

G4bool	IsBlocked (const G4Material *aMaterial) const

G4bool	IsBlocked (const G4Element *anElement) const

void	SetRecoilEnergyThreshold (G4double val)

G4double	GetRecoilEnergyThreshold () const

virtual const std::pair< G4double, G4double >	GetFatalEnergyCheckLevels () const

virtual std::pair< G4double, G4double >	GetEnergyMomentumCheckLevels () const

void	SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)

virtual void	BuildPhysicsTable (const G4ParticleDefinition &)

virtual void	InitialiseModel ()

	G4HadronicInteraction (const G4HadronicInteraction &right)=delete

const G4HadronicInteraction &	operator= (const G4HadronicInteraction &right)=delete

G4bool	operator== (const G4HadronicInteraction &right) const =delete

G4bool	operator!= (const G4HadronicInteraction &right) const =delete

Additional Inherited Members
Protected Member Functions inherited from G4HadronicInteraction
void	SetModelName (const G4String &nam)

G4bool	IsBlocked () const

void	Block ()

Protected Attributes inherited from G4HadronElastic
G4double	pLocalTmax

G4int	secID

Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange

G4int	verboseLevel

G4double	theMinEnergy

G4double	theMaxEnergy

G4bool	isBlocked

Detailed Description

Definition at line 51 of file G4NeutronElectronElModel.hh.

Constructor & Destructor Documentation

◆ G4NeutronElectronElModel()

G4NeutronElectronElModel::G4NeutronElectronElModel ( const G4String & name = "n-e-elastic" )

Definition at line 49 of file G4NeutronElectronElModel.cc.

  : G4HadronElastic(name)
{
  secID = G4PhysicsModelCatalog::GetModelID( "model_" + name );
  
 // neutron magneton squared
 
  fM   = neutron_mass_c2; // neutron mass
  fM2  = fM*fM;
  fme  = electron_mass_c2;
  fme2 = fme*fme;
  fMv2 = 0.7056*GeV*GeV;
 
  SetMinEnergy( 0.001*GeV );
  SetMaxEnergy( 10.*TeV );
  SetLowestEnergyLimit(1.e-6*eV);  
 
  theElectron = G4Electron::Electron();
  // PDG2016: sin^2 theta Weinberg
 
  fEnergyBin = 200;
  fMinEnergy = 1.*MeV;
  fMaxEnergy = 10000.*GeV;
  fEnergyVector = new G4PhysicsLogVector(fMinEnergy, fMaxEnergy, fEnergyBin, false);
 
  fAngleBin = 500;
  fAngleTable = 0;
 
  fCutEnergy = 0.; // default value
 
  Initialise();
}

Referenced by Initialise().

◆ ~G4NeutronElectronElModel()

G4NeutronElectronElModel::~G4NeutronElectronElModel ( )

virtual

Definition at line 84 of file G4NeutronElectronElModel.cc.

{
  if( fEnergyVector ) 
  {
    delete fEnergyVector;
    fEnergyVector = nullptr;
  }
  if( fAngleTable )
  {
    fAngleTable->clearAndDestroy();
    delete fAngleTable;
    fAngleTable = nullptr;
  }
}

Member Function Documentation

◆ ApplyYourself()

G4HadFinalState * G4NeutronElectronElModel::ApplyYourself	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	targetNucleus )

virtual

Reimplemented from G4HadronElastic.

Definition at line 262 of file G4NeutronElectronElModel.cc.

{
  theParticleChange.Clear();
 
  const G4HadProjectile* aParticle = &aTrack;
  G4double Tkin = aParticle->GetKineticEnergy();
  fAm = CalculateAm( Tkin);
  //   G4double En = aParticle->GetTotalEnergy();
 
  if( Tkin <= LowestEnergyLimit() ) 
  {
    theParticleChange.SetEnergyChange(Tkin);
    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
    return &theParticleChange;
  }
  // sample e-scattering angle and make final state in lab frame
 
  G4double sin2ht = SampleSin2HalfTheta( Tkin); // in n-rrest frame
 
  // G4cout<<"sin2ht = "<<sin2ht<<G4endl;
 
  G4double eTkin = fee; // fM;
 
  eTkin /= 1.+2.*fee*sin2ht/fM; // fme/En + 2*sin2ht;
 
  eTkin -= fme;
 
  // G4cout<<"eTkin = "<<eTkin<<G4endl;
 
  if( eTkin > fCutEnergy )
  {
    G4double ePlab = sqrt( eTkin*(eTkin + 2.*fme) );
 
    // G4cout<<"ePlab = "<<ePlab<<G4endl;
 
    G4double cost = 1. - 2*sin2ht;
 
    if( cost >  1. ) cost = 1.;
    if( cost < -1. ) cost = -1.;
 
    G4double sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
    G4double phi  = G4UniformRand()*CLHEP::twopi;
 
    G4ThreeVector eP( sint*std::cos(phi), sint*std::sin(phi), cost );
    eP *= ePlab;
    G4LorentzVector lvt2( eP, eTkin + electron_mass_c2 ); // recoil e- in n-rest frame
 
    G4LorentzVector lvp1 = aParticle->Get4Momentum();
    G4LorentzVector lvt1(0.,0.,0.,electron_mass_c2);
    G4LorentzVector lvsum = lvp1+lvt1;
 
    G4ThreeVector bst = lvp1.boostVector();
    lvt2.boost(bst);
 
    // G4cout<<"lvt2 = "<<lvt2<<G4endl;
 
    G4DynamicParticle * aSec = new G4DynamicParticle( theElectron, lvt2 );
    theParticleChange.AddSecondary( aSec, secID );
 
    G4LorentzVector lvp2 = lvsum-lvt2;
 
    // G4cout<<"lvp2 = "<<lvp2<<G4endl;
 
    G4double Tkin2 = lvp2.e()-aParticle->GetDefinition()->GetPDGMass();
    theParticleChange.SetEnergyChange(Tkin2);
    theParticleChange.SetMomentumChange(lvp2.vect().unit());
  }
  else if( eTkin > 0.0 ) 
  {
    theParticleChange.SetLocalEnergyDeposit( eTkin );
    Tkin -= eTkin;
 
    if( Tkin > 0. )
    {
      theParticleChange.SetEnergyChange( Tkin );
      theParticleChange.SetMomentumChange( aTrack.Get4Momentum().vect().unit() );
    }
  }
  else 
  {
    theParticleChange.SetEnergyChange( Tkin );
    theParticleChange.SetMomentumChange( aTrack.Get4Momentum().vect().unit() );
  }
  return &theParticleChange;
}

◆ CalculateAm()

G4double G4NeutronElectronElModel::CalculateAm ( G4double momentum )

inline

Definition at line 101 of file G4NeutronElectronElModel.hh.

{
  fee      = (Tkin+fM)*fme/fM;
    // G4cout<<"fee = "<<fee<<" MeV"<<G4endl;
  fee2     = fee*fee;
  G4double   momentum = std::sqrt( fee2 - fme2 );
  G4double k   = momentum/CLHEP::hbarc;
  G4double ch  = 1.13;
  G4double zn  = 1.77*k*CLHEP::Bohr_radius;
  G4double zn2 = zn*zn;
  fAm          = ch/zn2;
 
  return fAm;
}

Referenced by ApplyYourself(), and Initialise().

◆ GetAm()

G4double G4NeutronElectronElModel::GetAm ( )

inline

Definition at line 71 of file G4NeutronElectronElModel.hh.

71{return fAm;};

◆ GetCutEnergy()

G4double G4NeutronElectronElModel::GetCutEnergy ( )

inline

Definition at line 80 of file G4NeutronElectronElModel.hh.

80{return fCutEnergy;};

◆ GetTransfer()

G4double G4NeutronElectronElModel::GetTransfer	(	G4int	iTkin,
		G4int	iTransfer,
		G4double	position )

Definition at line 194 of file G4NeutronElectronElModel.cc.

{
  G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6;
 
  if( iTransfer == 0 ||  iTransfer == fAngleBin-1 )
  {
    randTransfer = (*fAngleTable)(iTkin)->Energy(iTransfer);
  }
  else
  {
    if ( iTransfer >= G4int((*fAngleTable)(iTkin)->GetVectorLength()) )
    {
      iTransfer = G4int((*fAngleTable)(iTkin)->GetVectorLength() - 1);
    }
    y1 = (*(*fAngleTable)(iTkin))(iTransfer-1);
    y2 = (*(*fAngleTable)(iTkin))(iTransfer);
 
    x1 = (*fAngleTable)(iTkin)->Energy(iTransfer-1);
    x2 = (*fAngleTable)(iTkin)->Energy(iTransfer);
 
    delta = y2 - y1;
    mean  = y2 + y1;
 
    if ( x1 == x2 ) randTransfer = x2;
    else
    {
      if ( delta < epsilon*mean ) 
      {
        randTransfer = x1 + ( x2 - x1 )*G4UniformRand();
      }
      else 
      {
        randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 );
      }
    }
  }
  return randTransfer;
}

Referenced by SampleSin2HalfTheta().

◆ Initialise()

void G4NeutronElectronElModel::Initialise ( )

Definition at line 120 of file G4NeutronElectronElModel.cc.

{
  G4double result = 0., sum, Tkin, dt, t1, t2;
  G4int iTkin, jTransfer;
  G4Integrator<G4NeutronElectronElModel, G4double(G4NeutronElectronElModel::*)(G4double)> integral;
 
  fAngleTable = new G4PhysicsTable(fEnergyBin);
 
  for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
  {
    Tkin  = fEnergyVector->GetLowEdgeEnergy(iTkin);
    fAm      = CalculateAm(Tkin);
    dt = 1./fAngleBin;
 
    G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fAngleBin);
 
    sum = 0.;
 
    for( jTransfer = 0; jTransfer < fAngleBin; jTransfer++)
    {
      t1 = dt*jTransfer;
      t2 = t1 + dt;
 
      result = integral.Legendre96( this, &G4NeutronElectronElModel::XscIntegrand, t1, t2 );
 
      sum += result;
      // G4cout<<sum<<", ";
      vectorT->PutValue(jTransfer, t1, sum);
    }
    // G4cout<<G4endl;   
    fAngleTable->insertAt(iTkin,vectorT);
  }
  return;
}

Referenced by G4NeutronElectronElModel().

◆ IsApplicable()

G4bool G4NeutronElectronElModel::IsApplicable	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	targetNucleus )

virtual

Reimplemented from G4HadronicInteraction.

Definition at line 110 of file G4NeutronElectronElModel.cc.

{
  G4String pName = aTrack.GetDefinition()->GetParticleName();
  G4double energy = aTrack.GetTotalEnergy();
 
  return (pName == "neutron" && energy >= fMinEnergy && energy <= fMaxEnergy);                            
}

◆ ModelDescription()

void G4NeutronElectronElModel::ModelDescription ( std::ostream & outFile ) const

virtual

Reimplemented from G4HadronElastic.

Definition at line 101 of file G4NeutronElectronElModel.cc.

{
  outFile << "G4NeutronElectronElModel is a neutrino-electron (neutral current) elastic scattering\n"
      << "model which uses the standard model \n"
      << "transfer parameterization.  The model is fully relativistic\n";
}

◆ SampleSin2HalfTheta()

G4double G4NeutronElectronElModel::SampleSin2HalfTheta ( G4double Tkin )

Definition at line 159 of file G4NeutronElectronElModel.cc.

{
  G4double result = 0., position; 
  G4int iTkin, iTransfer;
 
  for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
  {
    if( Tkin < fEnergyVector->Energy(iTkin) ) break;
  }  
  if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1;   // Tkin is more then theMaxEnergy
  if ( iTkin < 0 )           iTkin = 0; // against negative index, Tkin < theMinEnergy
 
    position = (*(*fAngleTable)(iTkin))(fAngleBin-1)*G4UniformRand();
 
    // G4cout<<"position = "<<position<<G4endl;
 
    for( iTransfer = 0; iTransfer < fAngleBin; iTransfer++)
    {
      if( position <= (*(*fAngleTable)(iTkin))(iTransfer) ) break;
    }
    if (iTransfer >= fAngleBin-1) iTransfer = fAngleBin-1;
 
    // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
 
    result = GetTransfer(iTkin, iTransfer, position);
 
    // G4cout<<"t = "<<t<<G4endl;
  
 
  return result;
}

Referenced by ApplyYourself().

◆ SetCutEnergy()

void G4NeutronElectronElModel::SetCutEnergy ( G4double ec )

inline

Definition at line 79 of file G4NeutronElectronElModel.hh.

79{fCutEnergy=ec;};

◆ XscIntegrand()

G4double G4NeutronElectronElModel::XscIntegrand ( G4double x )

Definition at line 239 of file G4NeutronElectronElModel.cc.

{
  G4double result = 1., q2, znq2, znf, znf2, znf4;
 
  znq2 = 1. + 2.*fee*x/fM;
 
  q2 = 4.*fee2*x/znq2;
 
  znf  = 1 + q2/fMv2;
  znf2 = znf*znf;
  znf4 = znf2*znf2;
 
  result /= ( x + fAm )*znq2*znq2*znf4; 
 
  result *= ( 1 - x )/( 1 + q2/4./fM2 ) + 2.*x;
 
  return result;
}

Referenced by Initialise().

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

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ G4NeutronElectronElModel()

◆ ~G4NeutronElectronElModel()

Member Function Documentation

◆ ApplyYourself()

◆ CalculateAm()

◆ GetAm()

◆ GetCutEnergy()

◆ GetTransfer()

◆ Initialise()

◆ IsApplicable()

◆ ModelDescription()

◆ SampleSin2HalfTheta()

◆ SetCutEnergy()

◆ XscIntegrand()