G4AdjointComptonModel Class Reference

#include <G4AdjointComptonModel.hh>

Inheritance diagram for G4AdjointComptonModel:

Public Member Functions
	G4AdjointComptonModel ()
	~G4AdjointComptonModel () override
void	SampleSecondaries (const G4Track &aTrack, G4bool isScatProjToProj, G4ParticleChange *fParticleChange) override
void	RapidSampleSecondaries (const G4Track &aTrack, G4bool isScatProjToProj, G4ParticleChange *fParticleChange)
G4double	DiffCrossSectionPerAtomPrimToScatPrim (G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.) override
G4double	DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.) override
G4double	GetSecondAdjEnergyMaxForScatProjToProj (G4double primAdjEnergy) override
G4double	GetSecondAdjEnergyMinForProdToProj (G4double primAdjEnergy) override
G4double	AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool isScatProjToProj) override
void	SetDirectProcess (G4VEmProcess *aProcess)
	G4AdjointComptonModel (G4AdjointComptonModel &)=delete
G4AdjointComptonModel &	operator= (const G4AdjointComptonModel &right)=delete
Public Member Functions inherited from G4VEmAdjointModel
	G4VEmAdjointModel (const G4String &nam)
virtual	~G4VEmAdjointModel ()
virtual void	SampleSecondaries (const G4Track &aTrack, G4bool isScatProjToProj, G4ParticleChange *fParticleChange)=0
virtual G4double	AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool isScatProjToProj)
virtual G4double	DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
virtual G4double	DiffCrossSectionPerAtomPrimToScatPrim (G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
virtual G4double	DiffCrossSectionPerVolumePrimToSecond (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
virtual G4double	DiffCrossSectionPerVolumePrimToScatPrim (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
virtual G4double	GetSecondAdjEnergyMaxForScatProjToProj (G4double primAdjEnergy)
virtual G4double	GetSecondAdjEnergyMinForScatProjToProj (G4double primAdjEnergy, G4double tcut=0.)
virtual G4double	GetSecondAdjEnergyMaxForProdToProj (G4double primAdjEnergy)
virtual G4double	GetSecondAdjEnergyMinForProdToProj (G4double primAdjEnergy)
void	DefineCurrentMaterial (const G4MaterialCutsCouple *couple)
std::vector< std::vector< double > * >	ComputeAdjointCrossSectionVectorPerAtomForSecond (G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
std::vector< std::vector< double > * >	ComputeAdjointCrossSectionVectorPerAtomForScatProj (G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
std::vector< std::vector< double > * >	ComputeAdjointCrossSectionVectorPerVolumeForSecond (G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
std::vector< std::vector< double > * >	ComputeAdjointCrossSectionVectorPerVolumeForScatProj (G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
void	SetCSMatrices (std::vector< G4AdjointCSMatrix * > *Vec1CSMatrix, std::vector< G4AdjointCSMatrix * > *Vec2CSMatrix)
G4ParticleDefinition *	GetAdjointEquivalentOfDirectPrimaryParticleDefinition ()
G4ParticleDefinition *	GetAdjointEquivalentOfDirectSecondaryParticleDefinition ()
G4double	GetHighEnergyLimit ()
G4double	GetLowEnergyLimit ()
void	SetHighEnergyLimit (G4double aVal)
void	SetLowEnergyLimit (G4double aVal)
void	DefineDirectEMModel (G4VEmModel *aModel)
void	SetAdjointEquivalentOfDirectPrimaryParticleDefinition (G4ParticleDefinition *aPart)
void	SetAdjointEquivalentOfDirectSecondaryParticleDefinition (G4ParticleDefinition *aPart)
void	SetSecondPartOfSameType (G4bool aBool)
G4bool	GetSecondPartOfSameType ()
void	SetUseMatrix (G4bool aBool)
void	SetUseMatrixPerElement (G4bool aBool)
void	SetUseOnlyOneMatrixForAllElements (G4bool aBool)
void	SetApplyCutInRange (G4bool aBool)
G4bool	GetUseMatrix ()
G4bool	GetUseMatrixPerElement ()
G4bool	GetUseOnlyOneMatrixForAllElements ()
G4bool	GetApplyCutInRange ()
G4String	GetName ()
virtual void	SetCSBiasingFactor (G4double aVal)
void	SetCorrectWeightForPostStepInModel (G4bool aBool)
void	SetAdditionalWeightCorrectionFactorForPostStepOutsideModel (G4double factor)
	G4VEmAdjointModel (G4VEmAdjointModel &)=delete
G4VEmAdjointModel &	operator= (const G4VEmAdjointModel &right)=delete

Additional Inherited Members
Protected Member Functions inherited from G4VEmAdjointModel
G4double	DiffCrossSectionFunction1 (G4double kinEnergyProj)
G4double	DiffCrossSectionFunction2 (G4double kinEnergyProj)
G4double	SampleAdjSecEnergyFromCSMatrix (std::size_t MatrixIndex, G4double prim_energy, G4bool isScatProjToProj)
G4double	SampleAdjSecEnergyFromCSMatrix (G4double prim_energy, G4bool isScatProjToProj)
void	SelectCSMatrix (G4bool isScatProjToProj)
virtual G4double	SampleAdjSecEnergyFromDiffCrossSectionPerAtom (G4double prim_energy, G4bool isScatProjToProj)
virtual void	CorrectPostStepWeight (G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool isScatProjToProj)
Protected Attributes inherited from G4VEmAdjointModel
G4AdjointCSManager *	fCSManager
G4VEmModel *	fDirectModel = nullptr
const G4String	fName
G4Material *	fSelectedMaterial = nullptr
G4Material *	fCurrentMaterial = nullptr
G4MaterialCutsCouple *	fCurrentCouple = nullptr
G4ParticleDefinition *	fAdjEquivDirectPrimPart = nullptr
G4ParticleDefinition *	fAdjEquivDirectSecondPart = nullptr
G4ParticleDefinition *	fDirectPrimaryPart = nullptr
std::vector< G4AdjointCSMatrix * > *	fCSMatrixProdToProjBackScat = nullptr
std::vector< G4AdjointCSMatrix * > *	fCSMatrixProjToProjBackScat = nullptr
std::vector< G4double >	fElementCSScatProjToProj
std::vector< G4double >	fElementCSProdToProj
G4double	fKinEnergyProdForIntegration = 0.
G4double	fKinEnergyScatProjForIntegration = 0.
G4double	fLastCS = 0.
G4double	fLastAdjointCSForScatProjToProj = 0.
G4double	fLastAdjointCSForProdToProj = 0.
G4double	fPreStepEnergy = 0.
G4double	fTcutPrim = 0.
G4double	fTcutSecond = 0.
G4double	fHighEnergyLimit = 0.
G4double	fLowEnergyLimit = 0.
G4double	fCsBiasingFactor = 1.
G4double	fOutsideWeightFactor = 1.
G4int	fASelectedNucleus = 0
G4int	fZSelectedNucleus = 0
std::size_t	fCSMatrixUsed = 0
G4bool	fSecondPartSameType = false
G4bool	fInModelWeightCorr
G4bool	fApplyCutInRange = true
G4bool	fUseMatrix = false
G4bool	fUseMatrixPerElement = false
G4bool	fOneMatrixForAllElements = false

Detailed Description

Definition at line 42 of file G4AdjointComptonModel.hh.

Constructor & Destructor Documentation

G4AdjointComptonModel() [1/2]

G4AdjointComptonModel::G4AdjointComptonModel

(

)

Definition at line 42 of file G4AdjointComptonModel.cc.

  : G4VEmAdjointModel("AdjointCompton")
{
  SetApplyCutInRange(false);
  SetUseMatrix(false);
  SetUseMatrixPerElement(true);
  SetUseOnlyOneMatrixForAllElements(true);
  fAdjEquivDirectPrimPart   = G4AdjointGamma::AdjointGamma();
  fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron();
  fDirectPrimaryPart        = G4Gamma::Gamma();
  fSecondPartSameType       = false;
  fDirectModel =
    new G4KleinNishinaCompton(G4Gamma::Gamma(), "ComptonDirectModel");
}

~G4AdjointComptonModel()

G4AdjointComptonModel::~G4AdjointComptonModel ( )

override

Definition at line 58 of file G4AdjointComptonModel.cc.

{}

G4AdjointComptonModel() [2/2]

G4AdjointComptonModel::G4AdjointComptonModel ( G4AdjointComptonModel &  )

delete

Member Function Documentation

AdjointCrossSection()

G4double G4AdjointComptonModel::AdjointCrossSection ( const G4MaterialCutsCouple *  aCouple, G4double  primEnergy, G4bool  isScatProjToProj  )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 325 of file G4AdjointComptonModel.cc.

{
  if(fUseMatrix)
    return G4VEmAdjointModel::AdjointCrossSection(aCouple, primEnergy,
                                                  isScatProjToProj);
  DefineCurrentMaterial(aCouple);
 
  G4float Cross     = 0.;
  G4float Emax_proj = 0.;
  G4float Emin_proj = 0.;
  if(!isScatProjToProj)
  {
    Emax_proj = GetSecondAdjEnergyMaxForProdToProj(primEnergy);
    Emin_proj = GetSecondAdjEnergyMinForProdToProj(primEnergy);
    if(Emax_proj > Emin_proj)
    {
      Cross = 0.1 *
              std::log((Emax_proj - G4float(primEnergy)) * Emin_proj /
                       Emax_proj / (Emin_proj - primEnergy)) *
              (1. + 2. * std::log(G4float(1. + electron_mass_c2 / primEnergy)));
    }
  }
  else
  {
    Emax_proj = GetSecondAdjEnergyMaxForScatProjToProj(primEnergy);
    Emin_proj = GetSecondAdjEnergyMinForScatProjToProj(primEnergy, 0.);
    if(Emax_proj > Emin_proj)
    {
      Cross = 0.1 * std::log(Emax_proj / Emin_proj);
    }
  }
 
  Cross *= fCurrentMaterial->GetElectronDensity() * twopi_mc2_rcl2;
  fLastCS = Cross;
  return double(Cross);
}

DiffCrossSectionPerAtomPrimToScatPrim()

G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim ( G4double  kinEnergyProj, G4double  kinEnergyScatProj, G4double  Z, G4double  A = 0.  )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 263 of file G4AdjointComptonModel.cc.

{
  // Based on Klein Nishina formula
  // In the forward case (see G4KleinNishinaCompton) the cross section is
  // parametrised but the secondaries are sampled from the Klein Nishina
  // differential cross section. The differential cross section used here
  // is therefore the cross section multiplied by the normalised
  // differential Klein Nishina cross section
 
  // Klein Nishina Cross Section
  G4double epsilon           = gamEnergy0 / electron_mass_c2;
  G4double one_plus_two_epsi = 1. + 2. * epsilon;
  if(gamEnergy1 > gamEnergy0 || gamEnergy1 < gamEnergy0 / one_plus_two_epsi)
  {
    return 0.;
  }
 
  G4double CS = std::log(one_plus_two_epsi) *
                (1. - 2. * (1. + epsilon) / (epsilon * epsilon));
  CS +=
    4. / epsilon + 0.5 * (1. - 1. / (one_plus_two_epsi * one_plus_two_epsi));
  CS /= epsilon;
  // Note that the pi*re2*Z factor is neglected because it is suppressed when
  // computing dCS_dE1/CS in the differential cross section
 
  // Klein Nishina Differential Cross Section
  G4double epsilon1 = gamEnergy1 / electron_mass_c2;
  G4double v        = epsilon1 / epsilon;
  G4double term1    = 1. + 1. / epsilon - 1. / epsilon1;
  G4double dCS_dE1  = 1. / v + v + term1 * term1 - 1.;
  dCS_dE1 *= 1. / epsilon / gamEnergy0;
 
  // Normalised to the CS used in G4
  fDirectCS = fDirectModel->ComputeCrossSectionPerAtom(
    G4Gamma::Gamma(), gamEnergy0, Z, 0., 0., 0.);
 
  dCS_dE1 *= fDirectCS / CS;
 
  return dCS_dE1;
}

Referenced by DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().

DiffCrossSectionPerAtomPrimToSecond()

G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond ( G4double  kinEnergyProj, G4double  kinEnergyProd, G4double  Z, G4double  A = 0.  )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 251 of file G4AdjointComptonModel.cc.

{
  G4double gamEnergy1   = gamEnergy0 - kinEnergyElec;
  G4double dSigmadEprod = 0.;
  if(gamEnergy1 > 0.)
    dSigmadEprod =
      DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0, gamEnergy1, Z, A);
  return dSigmadEprod;
}

GetSecondAdjEnergyMaxForScatProjToProj()

G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProj ( G4double  primAdjEnergy )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 306 of file G4AdjointComptonModel.cc.

{
  G4double inv_e_max = 1. / primAdjEnergy - 2. / electron_mass_c2;
  G4double e_max     = GetHighEnergyLimit();
  if(inv_e_max > 0.)
    e_max = std::min(1. / inv_e_max, e_max);
  return e_max;
}

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

GetSecondAdjEnergyMinForProdToProj()

G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProj ( G4double  primAdjEnergy )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 317 of file G4AdjointComptonModel.cc.

{
  G4double half_e = primAdjEnergy / 2.;
  return half_e + std::sqrt(half_e * (electron_mass_c2 + half_e));
}

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

operator=()

G4AdjointComptonModel & G4AdjointComptonModel::operator= ( const G4AdjointComptonModel &  right )

delete

RapidSampleSecondaries()

void G4AdjointComptonModel::RapidSampleSecondaries	(	const G4Track &	aTrack,
		G4bool	isScatProjToProj,
		G4ParticleChange *	fParticleChange
	)

Definition at line 137 of file G4AdjointComptonModel.cc.

{
  const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
 
  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
 
  if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
  {
    return;
  }
 
  G4double diffCSUsed =
    0.1 * fCurrentMaterial->GetElectronDensity() * twopi_mc2_rcl2;
  G4double gammaE1 = 0.;
  G4double gammaE2 = 0.;
  if(!isScatProjToProj)
  {
    G4double Emax = GetSecondAdjEnergyMaxForProdToProj(adjointPrimKinEnergy);
    G4double Emin = GetSecondAdjEnergyMinForProdToProj(adjointPrimKinEnergy);
    if(Emin >= Emax)
      return;
    G4double f1 = (Emin - adjointPrimKinEnergy) / Emin;
    G4double f2 = (Emax - adjointPrimKinEnergy) / Emax / f1;
    gammaE1 = adjointPrimKinEnergy / (1. - f1 * std::pow(f2, G4UniformRand()));
    gammaE2 = gammaE1 - adjointPrimKinEnergy;
    diffCSUsed =
      diffCSUsed *
      (1. + 2. * std::log(1. + electron_mass_c2 / adjointPrimKinEnergy)) *
      adjointPrimKinEnergy / gammaE1 / gammaE2;
  }
  else
  {
    G4double Emax =
      GetSecondAdjEnergyMaxForScatProjToProj(adjointPrimKinEnergy);
    G4double Emin =
      GetSecondAdjEnergyMinForScatProjToProj(adjointPrimKinEnergy, fTcutSecond);
    if(Emin >= Emax)
      return;
    gammaE2    = adjointPrimKinEnergy;
    gammaE1    = Emin * std::pow(Emax / Emin, G4UniformRand());
    diffCSUsed = diffCSUsed / gammaE1;
  }
 
  // Weight correction
  // First w_corr is set to the ratio between adjoint total CS and fwd total CS
  G4double w_corr = fOutsideWeightFactor;
  if(fInModelWeightCorr)
  {
    w_corr =
      G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
  }
  // Then another correction is needed because a biased differential CS has
  // been used rather than the one consistent with the direct model
 
  G4double diffCS =
    DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2, 1, 0.);
  if(diffCS > 0.)
    diffCS /= fDirectCS;  // here we have the normalised diffCS
  // And we remultiply by the lambda of the forward process
  diffCS *= fDirectProcess->GetCrossSection(gammaE1, fCurrentCouple);
 
  w_corr *= diffCS / diffCSUsed;
 
  G4double new_weight = aTrack.GetWeight() * w_corr;
  fParticleChange->SetParentWeightByProcess(false);
  fParticleChange->SetSecondaryWeightByProcess(false);
  fParticleChange->ProposeParentWeight(new_weight);
 
  G4double cos_th = 1. + electron_mass_c2 * (1. / gammaE1 - 1. / gammaE2);
  if(!isScatProjToProj)
  {
    G4double p_elec = theAdjointPrimary->GetTotalMomentum();
    cos_th          = (gammaE1 - gammaE2 * cos_th) / p_elec;
  }
  G4double sin_th = 0.;
  if(std::abs(cos_th) > 1.)
  {
    if(cos_th > 0.)
    {
      cos_th = 1.;
    }
    else
      cos_th = -1.;
  }
  else
    sin_th = std::sqrt(1. - cos_th * cos_th);
 
  // gamma0 momentum
  G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection();
  G4double phi               = G4UniformRand() * twopi;
  G4ThreeVector gammaMomentum1 =
    gammaE1 *
    G4ThreeVector(std::cos(phi) * sin_th, std::sin(phi) * sin_th, cos_th);
  gammaMomentum1.rotateUz(dir_parallel);
 
  if(!isScatProjToProj)
  {  // kill the primary and add a secondary
    fParticleChange->ProposeTrackStatus(fStopAndKill);
    fParticleChange->AddSecondary(
      new G4DynamicParticle(fAdjEquivDirectPrimPart, gammaMomentum1));
  }
  else
  {
    fParticleChange->ProposeEnergy(gammaE1);
    fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
  }
}

Referenced by SampleSecondaries().

SampleSecondaries()

void G4AdjointComptonModel::SampleSecondaries ( const G4Track &  aTrack, G4bool  isScatProjToProj, G4ParticleChange *  fParticleChange  )

overridevirtual

Implements G4VEmAdjointModel.

Definition at line 61 of file G4AdjointComptonModel.cc.

{
  if(!fUseMatrix)
    return RapidSampleSecondaries(aTrack, isScatProjToProj, fParticleChange);
 
  // A recall of the compton scattering law:
  // Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th))
  // Therefore Egamma2_max= Egamma2(cos_th=1) = Egamma1
  // and       Egamma2_min= Egamma2(cos_th=-1) =
  //             Egamma1/(1+2.(Egamma1/E0_electron))
  const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
  if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
  {
    return;
  }
 
  // Sample secondary energy
  G4double gammaE1;
  gammaE1 =
    SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, isScatProjToProj);
 
  // gammaE2
  G4double gammaE2 = adjointPrimKinEnergy;
  if(!isScatProjToProj)
    gammaE2 = gammaE1 - adjointPrimKinEnergy;
 
  // Cos th
  G4double cos_th = 1. + electron_mass_c2 * (1. / gammaE1 - 1. / gammaE2);
  if(!isScatProjToProj)
  {
    cos_th =
      (gammaE1 - gammaE2 * cos_th) / theAdjointPrimary->GetTotalMomentum();
  }
  G4double sin_th = 0.;
  if(std::abs(cos_th) > 1.)
  {
    if(cos_th > 0.)
    {
      cos_th = 1.;
    }
    else
      cos_th = -1.;
    sin_th = 0.;
  }
  else
    sin_th = std::sqrt(1. - cos_th * cos_th);
 
  // gamma0 momentum
  G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection();
  G4double phi               = G4UniformRand() * twopi;
  G4ThreeVector gammaMomentum1 =
    gammaE1 *
    G4ThreeVector(std::cos(phi) * sin_th, std::sin(phi) * sin_th, cos_th);
  gammaMomentum1.rotateUz(dir_parallel);
 
  // correct the weight of particles before adding the secondary
  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(),
                        adjointPrimKinEnergy, gammaE1, isScatProjToProj);
 
  if(!isScatProjToProj)
  {  // kill the primary and add a secondary
    fParticleChange->ProposeTrackStatus(fStopAndKill);
    fParticleChange->AddSecondary(
      new G4DynamicParticle(fAdjEquivDirectPrimPart, gammaMomentum1));
  }
  else
  {
    fParticleChange->ProposeEnergy(gammaE1);
    fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
  }
}

SetDirectProcess()

void G4AdjointComptonModel::SetDirectProcess ( G4VEmProcess *  aProcess )

inline

Definition at line 72 of file G4AdjointComptonModel.hh.

  {
    fDirectProcess = aProcess;
  };

---

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

• G4AdjointComptonModel.hh
• G4AdjointComptonModel.cc