G4AdjointhIonisationModel Class Reference

#include <G4AdjointhIonisationModel.hh>

Inheritance diagram for G4AdjointhIonisationModel:

Public Member Functions
	G4AdjointhIonisationModel (G4ParticleDefinition *pDef)
	~G4AdjointhIonisationModel () override
void	SampleSecondaries (const G4Track &aTrack, G4bool isScatProjToProj, G4ParticleChange *fParticleChange) override
void	RapidSampleSecondaries (const G4Track &aTrack, G4bool isScatProjToProj, G4ParticleChange *fParticleChange)
G4double	DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.) override
G4double	AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool isScatProjToProj) override
G4double	GetSecondAdjEnergyMaxForScatProjToProj (G4double primAdjEnergy) override
G4double	GetSecondAdjEnergyMinForScatProjToProj (G4double primAdjEnergy, G4double tcut=0.) override
G4double	GetSecondAdjEnergyMaxForProdToProj (G4double primAdjEnergy) override
G4double	GetSecondAdjEnergyMinForProdToProj (G4double primAdjEnergy) override
	G4AdjointhIonisationModel (G4AdjointhIonisationModel &)=delete
G4AdjointhIonisationModel &	operator= (const G4AdjointhIonisationModel &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 47 of file G4AdjointhIonisationModel.hh.

Constructor & Destructor Documentation

G4AdjointhIonisationModel() [1/2]

G4AdjointhIonisationModel::G4AdjointhIonisationModel ( G4ParticleDefinition *  pDef )

explicit

Definition at line 42 of file G4AdjointhIonisationModel.cc.

  : G4VEmAdjointModel("Adjoint_hIonisation")
{
  fUseMatrix               = true;
  fUseMatrixPerElement     = true;
  fApplyCutInRange         = true;
  fOneMatrixForAllElements = true;
  fSecondPartSameType      = false;
 
  // The direct EM Model is taken as BetheBloch. It is only used for the
  // computation of the differential cross section.
  // The Bragg model could be used as an alternative as it offers the same
  // differential cross section
 
  fDirectModel              = new G4BetheBlochModel(pDef);
  fBraggDirectEMModel       = new G4BraggModel(pDef);
  fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron();
  fDirectPrimaryPart        = pDef;
 
  if(pDef == G4Proton::Proton())
  {
    fAdjEquivDirectPrimPart = G4AdjointProton::AdjointProton();
  }
 
  DefineProjectileProperty();
}

~G4AdjointhIonisationModel()

G4AdjointhIonisationModel::~G4AdjointhIonisationModel ( )

override

Definition at line 70 of file G4AdjointhIonisationModel.cc.

{}

G4AdjointhIonisationModel() [2/2]

G4AdjointhIonisationModel::G4AdjointhIonisationModel ( G4AdjointhIonisationModel &  )

delete

Member Function Documentation

AdjointCrossSection()

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

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 397 of file G4AdjointhIonisationModel.cc.

{
  if(fUseMatrix)
    return G4VEmAdjointModel::AdjointCrossSection(aCouple, primEnergy,
                                                  isScatProjToProj);
  DefineCurrentMaterial(aCouple);
 
  G4double Cross =
    fCurrentMaterial->GetElectronDensity() * twopi_mc2_rcl2 * fMass;
 
  if(!isScatProjToProj)
  {
    G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(primEnergy);
    G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(primEnergy);
    if(Emax_proj > Emin_proj && primEnergy > fTcutSecond)
    {
      Cross *= (1. / Emin_proj - 1. / Emax_proj) / primEnergy;
    }
    else
      Cross = 0.;
  }
  else
  {
    G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProj(primEnergy);
    G4double Emin_proj =
      GetSecondAdjEnergyMinForScatProjToProj(primEnergy, fTcutSecond);
    G4double diff1 = Emin_proj - primEnergy;
    G4double diff2 = Emax_proj - primEnergy;
    G4double t1 =
      (1. / diff1 + 1. / Emin_proj - 1. / diff2 - 1. / Emax_proj) / primEnergy;
    G4double t2 =
      2. * std::log(Emax_proj / Emin_proj) / primEnergy / primEnergy;
    Cross *= (t1 + t2);
  }
  fLastCS = Cross;
  return Cross;
}

DiffCrossSectionPerAtomPrimToSecond()

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

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 280 of file G4AdjointhIonisationModel.cc.

{  // Probably here the Bragg Model should be also used for
   // kinEnergyProj/nuc < 2 MeV
 
  G4double dSigmadEprod = 0.;
  G4double Emax_proj    = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd);
  G4double Emin_proj    = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd);
 
  // the produced particle should have a kinetic energy smaller than the
  // projectile
  if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj)
  {
    G4double Tmax = kinEnergyProj;
    G4double E1   = kinEnergyProd;
    //1.0006 factor seems to give the best diff CS, important impact on proton correction factor
    G4double E2   = kinEnergyProd *1.0006;
    G4double sigma1, sigma2;
    if(kinEnergyProj > 2. * MeV)
    {
      sigma1 = fDirectModel->ComputeCrossSectionPerAtom(
        fDirectPrimaryPart, kinEnergyProj, Z, A, E1, 1.e20);
      sigma2 = fDirectModel->ComputeCrossSectionPerAtom(
        fDirectPrimaryPart, kinEnergyProj, Z, A, E2, 1.e20);
    }
    else
    {
      sigma1 = fBraggDirectEMModel->ComputeCrossSectionPerAtom(
        fDirectPrimaryPart, kinEnergyProj, Z, A, E1, 1.e20);
      sigma2 = fBraggDirectEMModel->ComputeCrossSectionPerAtom(
        fDirectPrimaryPart, kinEnergyProj, Z, A, E2, 1.e20);
    }
 
    dSigmadEprod = (sigma1 - sigma2) / (E2 - E1);
    if(dSigmadEprod > 1.)
    {
      G4cout << "sigma1 " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
             << '\t' << sigma1 << G4endl;
      G4cout << "sigma2 " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
             << '\t' << sigma2 << G4endl;
      G4cout << "dsigma " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
             << '\t' << dSigmadEprod << G4endl;
    }
 
    // correction of differential cross section at high energy to correct for
    // the suppression of particle at secondary at high energy used in the Bethe
    // Bloch Model. This correction consists of multiplying by g the probability
    // function used to test the rejection of a secondary. Source code taken
    // from G4BetheBlochModel::SampleSecondaries
    G4double deltaKinEnergy = kinEnergyProd;
 
    // projectile formfactor - suppression of high energy
    // delta-electron production at high energy
    G4double x = fFormFact * deltaKinEnergy;
    if(x > 1.e-6)
    {
      G4double totEnergy = kinEnergyProj + fMass;
      G4double etot2     = totEnergy * totEnergy;
      G4double beta2 = kinEnergyProj * (kinEnergyProj + 2.0 * fMass) / etot2;
      G4double f     = 1.0 - beta2 * deltaKinEnergy / Tmax;
      G4double f1    = 0.0;
      if(0.5 == fSpin)
      {
        f1 = 0.5 * deltaKinEnergy * deltaKinEnergy / etot2;
        f += f1;
      }
      G4double x1 = 1.0 + x;
      G4double gg = 1.0 / (x1 * x1);
      if(0.5 == fSpin)
      {
        G4double x2 = 0.5 * electron_mass_c2 * deltaKinEnergy / (fMass * fMass);
        gg *= (1.0 + fMagMoment2 * (x2 - f1 / f) / (1.0 + x2));
      }
      if(gg > 1.0)
      {
        G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
               << G4endl;
        gg = 1.;
      }
      dSigmadEprod *= gg;
    }
  }
 
  return dSigmadEprod;
}

Referenced by RapidSampleSecondaries().

GetSecondAdjEnergyMaxForProdToProj()

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProj ( G4double  primAdjEnergy )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 454 of file G4AdjointhIonisationModel.cc.

{
  return GetHighEnergyLimit();
}

Referenced by AdjointCrossSection(), DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().

GetSecondAdjEnergyMaxForScatProjToProj()

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProj ( G4double  primAdjEnergy )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 438 of file G4AdjointhIonisationModel.cc.

{
  G4double Tmax = primAdjEnergy * fOnePlusRatio2 /
                  (fOneMinusRatio2 - 2. * fMassRatio * primAdjEnergy / fMass);
  return Tmax;
}

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

GetSecondAdjEnergyMinForProdToProj()

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProj ( G4double  primAdjEnergy )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 460 of file G4AdjointhIonisationModel.cc.

{
  G4double Tmin =
    (2. * primAdjEnergy - 4. * fMass +
     std::sqrt(4. * primAdjEnergy * primAdjEnergy + 16. * fMass * fMass +
               8. * primAdjEnergy * fMass * (1. / fMassRatio + fMassRatio))) /
    4.;
  return Tmin;
}

Referenced by AdjointCrossSection(), DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().

GetSecondAdjEnergyMinForScatProjToProj()

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProj ( G4double  primAdjEnergy, G4double  tcut = 0.  )

overridevirtual

Reimplemented from G4VEmAdjointModel.

Definition at line 447 of file G4AdjointhIonisationModel.cc.

{
  return primAdjEnergy + tcut;
}

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

operator=()

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

delete

RapidSampleSecondaries()

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

Definition at line 143 of file G4AdjointhIonisationModel.cc.

{
  const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
 
  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
  G4double adjointPrimP         = theAdjointPrimary->GetTotalMomentum();
 
  if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
  {
    return;
  }
 
  G4double projectileKinEnergy = 0.;
  G4double eEnergy             = 0.;
  G4double newCS =
    fCurrentMaterial->GetElectronDensity() * twopi_mc2_rcl2 * fMass;
  if(!isScatProjToProj)
  {  // 1/E^2 distribution
 
    eEnergy       = adjointPrimKinEnergy;
    G4double Emax = GetSecondAdjEnergyMaxForProdToProj(adjointPrimKinEnergy);
    G4double Emin = GetSecondAdjEnergyMinForProdToProj(adjointPrimKinEnergy);
    if(Emin >= Emax)
      return;
    G4double a = 1. / Emax;
    G4double b = 1. / Emin;
    newCS      = newCS * (b - a) / eEnergy;
 
    projectileKinEnergy = 1. / (b - (b - a) * G4UniformRand());
  }
  else
  {
    G4double Emax =
      GetSecondAdjEnergyMaxForScatProjToProj(adjointPrimKinEnergy);
    G4double Emin =
      GetSecondAdjEnergyMinForScatProjToProj(adjointPrimKinEnergy, fTcutSecond);
    if(Emin >= Emax)
      return;
    G4double diff1 = Emin - adjointPrimKinEnergy;
    G4double diff2 = Emax - adjointPrimKinEnergy;
 
    G4double t1    = adjointPrimKinEnergy * (1. / diff1 - 1. / diff2);
    G4double t2    = adjointPrimKinEnergy * (1. / Emin - 1. / Emax);
    G4double t3    = 2. * std::log(Emax / Emin);
    G4double sum_t = t1 + t2 + t3;
    newCS      = newCS * sum_t / adjointPrimKinEnergy / adjointPrimKinEnergy;
    G4double t = G4UniformRand() * sum_t;
    if(t <= t1)
    {
      G4double q          = G4UniformRand() * t1 / adjointPrimKinEnergy;
      projectileKinEnergy = adjointPrimKinEnergy + 1. / (1. / diff1 - q);
    }
    else if(t <= t2)
    {
      G4double q          = G4UniformRand() * t2 / adjointPrimKinEnergy;
      projectileKinEnergy = 1. / (1. / Emin - q);
    }
    else
    {
      projectileKinEnergy = Emin * std::pow(Emax / Emin, G4UniformRand());
    }
    eEnergy = projectileKinEnergy - adjointPrimKinEnergy;
  }
 
  G4double diffCS_perAtom_Used = twopi_mc2_rcl2 * fMass * adjointPrimKinEnergy /
                                 projectileKinEnergy / projectileKinEnergy /
                                 eEnergy / eEnergy;
 
  // Weight correction
  // First w_corr is set to the ratio between adjoint total CS and fwd total CS
  G4double w_corr =
    G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
 
  w_corr *= newCS / fLastCS;
  // Then another correction is needed due to the fact that a biaised
  // differential CS has been used rather than the one consistent with the
  // direct model. Here we consider the true diffCS as the one obtained by the
  // numerical differentiation over Tcut of the direct CS
 
  G4double diffCS =
    DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy, 1, 1);
  w_corr *= diffCS / diffCS_perAtom_Used;
 
  if (isScatProjToProj && fTcutSecond>0.005) w_corr=1.;
 
  G4double new_weight = aTrack.GetWeight() * w_corr;
  fParticleChange->SetParentWeightByProcess(false);
  fParticleChange->SetSecondaryWeightByProcess(false);
  fParticleChange->ProposeParentWeight(new_weight);
 
  // Kinematic:
  // we consider a two body elastic scattering for the forward processes where
  // the projectile knocks on an e- at rest and gives it part of its energy
  G4double projectileM0          = fAdjEquivDirectPrimPart->GetPDGMass();
  G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy;
  G4double projectileP2 =
    projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0;
 
  // Companion
  G4double companionM0 = fAdjEquivDirectPrimPart->GetPDGMass();
  if(isScatProjToProj)
  {
    companionM0 = fAdjEquivDirectSecondPart->GetPDGMass();
  }
  G4double companionTotalEnergy =
    companionM0 + projectileKinEnergy - adjointPrimKinEnergy;
  G4double companionP2 =
    companionTotalEnergy * companionTotalEnergy - companionM0 * companionM0;
 
  // Projectile momentum
  G4double P_parallel =
    (adjointPrimP * adjointPrimP + projectileP2 - companionP2) /
    (2. * adjointPrimP);
  G4double P_perp = std::sqrt(projectileP2 - P_parallel * P_parallel);
  G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection();
  G4double phi               = G4UniformRand() * twopi;
  G4ThreeVector projectileMomentum =
    G4ThreeVector(P_perp * std::cos(phi), P_perp * std::sin(phi), P_parallel);
  projectileMomentum.rotateUz(dir_parallel);
 
  if(!isScatProjToProj)
  {  // kill the primary and add a secondary
    fParticleChange->ProposeTrackStatus(fStopAndKill);
    fParticleChange->AddSecondary(
      new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum));
  }
  else
  {
    fParticleChange->ProposeEnergy(projectileKinEnergy);
    fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
  }
}

Referenced by SampleSecondaries().

SampleSecondaries()

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

overridevirtual

Implements G4VEmAdjointModel.

Definition at line 73 of file G4AdjointhIonisationModel.cc.

{
  if(!fUseMatrix)
    return RapidSampleSecondaries(aTrack, isScatProjToProj, fParticleChange);
 
  const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
 
  // Elastic inverse scattering
  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
  G4double adjointPrimP         = theAdjointPrimary->GetTotalMomentum();
 
  if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
  {
    return;
  }
 
  // Sample secondary energy
  G4double projectileKinEnergy =
    SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, isScatProjToProj);
  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(),
                        adjointPrimKinEnergy, projectileKinEnergy,
                        isScatProjToProj);
  // Caution!!! this weight correction should be always applied
 
  // Kinematic:
  // we consider a two body elastic scattering for the forward processes where
  // the projectile knock on an e- at rest and gives it part of its energy
  G4double projectileM0          = fAdjEquivDirectPrimPart->GetPDGMass();
  G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy;
  G4double projectileP2 =
    projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0;
 
  // Companion
  G4double companionM0 = fAdjEquivDirectPrimPart->GetPDGMass();
  if(isScatProjToProj)
  {
    companionM0 = fAdjEquivDirectSecondPart->GetPDGMass();
  }
  G4double companionTotalEnergy =
    companionM0 + projectileKinEnergy - adjointPrimKinEnergy;
  G4double companionP2 =
    companionTotalEnergy * companionTotalEnergy - companionM0 * companionM0;
 
  // Projectile momentum
  G4double P_parallel =
    (adjointPrimP * adjointPrimP + projectileP2 - companionP2) /
    (2. * adjointPrimP);
  G4double P_perp = std::sqrt(projectileP2 - P_parallel * P_parallel);
  G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection();
  G4double phi               = G4UniformRand() * twopi;
  G4ThreeVector projectileMomentum =
    G4ThreeVector(P_perp * std::cos(phi), P_perp * std::sin(phi), P_parallel);
  projectileMomentum.rotateUz(dir_parallel);
 
  if(!isScatProjToProj)
  {  // kill the primary and add a secondary
    fParticleChange->ProposeTrackStatus(fStopAndKill);
    fParticleChange->AddSecondary(
      new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum));
  }
  else
  {
    fParticleChange->ProposeEnergy(projectileKinEnergy);
    fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
  }
}

---

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

• G4AdjointhIonisationModel.hh
• G4AdjointhIonisationModel.cc