Concrete monopole ionisation model. More...

#include <G4mplIonisationWithDeltaModel.h>

Inheritance diagram for G4mplIonisationWithDeltaModel:

Public Member Functions
	G4mplIonisationWithDeltaModel (G4double mCharge, const G4String &nam="mplIonisationWithDelta")
	Constructor.

virtual	~G4mplIonisationWithDeltaModel ()
	Destructor.

virtual void	Initialise (const G4ParticleDefinition *p, const G4DataVector &) override
	G4VEmModel::Initialise implementation.

virtual G4double	ComputeDEDXPerVolume (const G4Material material, const G4ParticleDefinition p, G4double kineticEnergy, G4double maxEnergy) override
	G4VEmModel::ComputeDEDXPerVolume implementation.

virtual G4double	ComputeCrossSectionPerElectron (const G4ParticleDefinition *p, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
	Compute cross section per electron for delta electrons emission.

virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *p, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
	Compute cross section per atom for delta electrons emission.

virtual void	SampleSecondaries (std::vector< G4DynamicParticle * > vdp, const G4MaterialCutsCouple materialCutsCouple, const G4DynamicParticle *dp, G4double tmin, G4double maxEnergy) override
	Create the sample of secondary delta electrons.

virtual G4double	SampleFluctuations (const G4MaterialCutsCouple couple, const G4DynamicParticle dp, G4double tcut, G4double tmax, G4double length, G4double meanLoss) override
	Create fluctuations in the energies lost to a secondary delta electron.

virtual G4double	Dispersion (const G4Material material, const G4DynamicParticle , G4double tcut, G4double tmax, G4double length) override
	Calculate dispersion.

virtual G4double	MinEnergyCut (const G4ParticleDefinition , const G4MaterialCutsCouple couple) override
	Threshold for zero value.

void	SetParticle (const G4ParticleDefinition *p)
	Read definition of the monopole.

Protected Member Functions
virtual G4double	MaxSecondaryEnergy (const G4ParticleDefinition *p, G4double kinEnergy) override
	Calculate maximum energy available for secondary particle emission.

Private Member Functions
G4double	ComputeDEDXAhlen (const G4Material *material, G4double bg2, G4double cut)
	Calculate dedx based on extrapolated Ahlen formula.

G4mplIonisationWithDeltaModel &	operator= (const G4mplIonisationWithDeltaModel &right)=delete
	Assignment operator should be hidden.

	G4mplIonisationWithDeltaModel (const G4mplIonisationWithDeltaModel &copy)=delete
	Copy constructor should be hidden.

Private Attributes
const G4ParticleDefinition *	monopole
	Monopole definition.

G4ParticleDefinition *	theElectron
	Electron definition.

G4ParticleChangeForLoss *	fParticleChange
	Pointer to ionising particle.

G4double	mass
	Mass of the monopole.

G4double	magCharge
	Monopole magnetic charge in e+ units.

G4double	twoln10
	log(100.0)

G4double	betalow
	Beta threshold for low asymptotic.

G4double	betalim
	Beta threshold for high energy (only Ahlen formula)

G4double	beta2lim
	Square of betalim.

G4double	bg2lim
	(beta*gamma)^2 for betalim

G4double	chargeSquare
	Square of magnetic charge in units of Dirac charge.

G4double	dedxlim
	dedx limit in asymptotic formula, not used

G4double	pi_hbarc2_over_mc2
	Convenient constants combination with mass.

G4double	nmpl
	Magn.

Static Private Attributes
static std::vector< G4double > *	dedx0 = nullptr
	Base dedx for each couple in current material.

Detailed Description

Concrete monopole ionisation model.

Monopole magnetic ionisation dedx is similar to electric one with the exception of beta*beta factor. This model is also capable of producing secondary delta electrons.

Definition at line 34 of file G4mplIonisationWithDeltaModel.h.

Constructor & Destructor Documentation

◆ G4mplIonisationWithDeltaModel() [1/2]

G4mplIonisationWithDeltaModel	(	G4double	mCharge,
		const G4String &	nam = `"mplIonisationWithDelta"`
	)

explicit

Constructor.

Parameters

mCharge	Magnetic charge of the monopole, in e+ units
nam	Name of the ionisation process in GEANT4

Definition at line 41 of file G4mplIonisationWithDeltaModel.cc.

  : G4VEmModel(nam), G4VEmFluctuationModel(nam),
    magCharge(mCharge),
    twoln10(log(100.0)),
    betalow(0.01),
    betalim(0.1),
    beta2lim(betalim * betalim),
    bg2lim(beta2lim * (1.0 + beta2lim)),
    pi_hbarc2_over_mc2(pi * hbarc * hbarc / electron_mass_c2)
{
  nmpl = magCharge * 2 * fine_structure_const;
  chargeSquare = magCharge * magCharge * 4 * fine_structure_const *
                 fine_structure_const;
  //NOTE Formulas below assume Dirac charge units for magnetic charge, g_D = 68.5e
  dedxlim = 45. * chargeSquare * GeV * cm2 / g;
  fParticleChange = nullptr;
  theElectron = G4Electron::Electron();
  B2INFO("Monopole ionisation model with d-electron production, Gmag= "  << magCharge / eplus);
  if (nmpl >= 6)
    B2WARNING("Monopole charge Gmag= " << magCharge / eplus << "e not reasonable. Please choose a value smaller than 411e.");
  monopole = nullptr;
  mass = 0.0;
}

◆ ~G4mplIonisationWithDeltaModel()

~G4mplIonisationWithDeltaModel ( )

virtual

Destructor.

Definition at line 66 of file G4mplIonisationWithDeltaModel.cc.

{
  if (IsMaster()) { delete dedx0; }
}

◆ G4mplIonisationWithDeltaModel() [2/2]

G4mplIonisationWithDeltaModel ( const G4mplIonisationWithDeltaModel & copy )

privatedelete

Copy constructor should be hidden.

Parameters

copy	Copy reference

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

G4double ComputeCrossSectionPerAtom	(	const G4ParticleDefinition *	p,
		G4double	kineticEnergy,
		G4double	Z,
		G4double	A,
		G4double	cutEnergy,
		G4double	maxEnergy
	)

overridevirtual

Compute cross section per atom for delta electrons emission.

Parameters

p	Pointer to monopole definition
kineticEnergy	Monopole kinetic energy value
Z	Atomic charge
A	This argument of base function is ignored
cutEnergy	Cut energy for ionisation process
maxEnergy	Maximum energy of secondary electron

Returns: Cross section

Definition at line 214 of file G4mplIonisationWithDeltaModel.cc.

{
  G4double cross =
    Z * ComputeCrossSectionPerElectron(p, kineticEnergy, cutEnergy, maxEnergy);
  return cross;
}

◆ ComputeCrossSectionPerElectron()

G4double ComputeCrossSectionPerElectron	(	const G4ParticleDefinition *	p,
		G4double	kineticEnergy,
		G4double	cutEnergy,
		G4double	maxEnergy
	)

virtual

Compute cross section per electron for delta electrons emission.

Parameters

p	Pointer to monopole definition
kineticEnergy	Monopole kinetic energy value
cutEnergy	Cut energy for ionisation process
maxEnergy	Maximum energy of secondary electron

Returns: Cross section

Definition at line 196 of file G4mplIonisationWithDeltaModel.cc.

{
  if (!monopole) { SetParticle(p); }
  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
  G4double maxEnergy = std::min(tmax, maxKinEnergy);
  G4double cutEnergy = std::max(LowEnergyLimit(), cut);
  G4double cross = (cutEnergy < maxEnergy)
                   ? (1.0 / cutEnergy - 1.0 / maxEnergy) * twopi_mc2_rcl2 * chargeSquare : 0.0; //cross section used in ATLAS
//   G4double cross = (cutEnergy < maxEnergy)
//                    ? (0.5 / cutEnergy - 0.5 / maxEnergy) * pi_hbarc2_over_mc2 * chargeSquare : 0.0; //one coming with GEANT4; causes 4 orders higher values
  return cross;
}

◆ ComputeDEDXAhlen()

G4double ComputeDEDXAhlen	(	const G4Material *	material,
		G4double	bg2,
		G4double	cut
	)

private

Calculate dedx based on extrapolated Ahlen formula.

References [1] Steven P. Ahlen: Energy loss of relativistic heavy ionizing particles, S.P. Ahlen, Rev. Mod. Phys 52(1980), p121 [2] K.A. Milton arXiv:hep-ex/0602040 [3] S.P. Ahlen and K. Kinoshita, Phys. Rev. D26 (1982) 2347

Parameters

material	Pointer to the material of current volume
bg2	(beta*gamma)^2
cut	Cut energy for ionisation process

Returns: dedx magnitude

Definition at line 161 of file G4mplIonisationWithDeltaModel.cc.

{
  G4double eDensity = material->GetElectronDensity();
  G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
 
  // Ahlen's formula for nonconductors, [1]p157, f(5.7)
  G4double dedx =
    0.5 * (G4Log(2.0 * electron_mass_c2 * bg2 * cutEnergy / (eexc * eexc)) - 1.0);//"Conventional" ionisation
//   G4double dedx =
//     1.0 * (G4Log(2.0 * electron_mass_c2 * bg2 * cutEnergy / (eexc * eexc)));//Fryberger magneticon double ionisation
 
 
  G4double k = 0;   // Cross-section correction
  if (nmpl >= 0.5) { k = 0.406; }
  if (nmpl >= 1) { k = 0.346; }
  if (nmpl >= 1.5) { k = 0.3; }
  const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};   // Bloch correction
  if (nmpl < 6)
    dedx += 0.5 * k - B[int(floor(nmpl + 0.5))];
 
 
  // density effect correction
  G4double x = G4Log(bg2) / twoln10;
  dedx -= material->GetIonisation()->DensityCorrection(x);
 
  // now compute the total ionization loss
  dedx *=  pi_hbarc2_over_mc2 * eDensity * chargeSquare;
 
  dedx = std::max(dedx, 0.0);
  return dedx;
}

◆ ComputeDEDXPerVolume()

G4double ComputeDEDXPerVolume	(	const G4Material *	material,
		const G4ParticleDefinition *	p,
		G4double	kineticEnergy,
		G4double	maxEnergy
	)

overridevirtual

G4VEmModel::ComputeDEDXPerVolume implementation.

Combines dedx estimation from extrapolated Ahlen formula and formula for low beta.

Parameters

material	Pointer to the material of current volume
p	Pointer to monopole definition
kineticEnergy	Monopole kinetic energy value
maxEnergy	Maximum energy of secondary electron

Returns: dedx magnitude

Definition at line 119 of file G4mplIonisationWithDeltaModel.cc.

{
  if (!monopole) { SetParticle(p); }
  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
  G4double cutEnergy = std::min(tmax, maxEnergy);
  cutEnergy = std::max(LowEnergyLimit(), cutEnergy);
  G4double tau   = kineticEnergy / mass;
  G4double gam   = tau + 1.0;
  G4double bg2   = tau * (tau + 2.0);
  G4double beta2 = bg2 / (gam * gam);
  G4double beta  = sqrt(beta2);
 
  // low-energy asymptotic formula
  //G4double dedx  = dedxlim*beta*material->GetDensity();
  G4double dedx = (*dedx0)[CurrentCouple()->GetIndex()] * beta;
 
  // above asymptotic
  if (beta > betalow) {
 
    // high energy
    if (beta >= betalim) {
      dedx = ComputeDEDXAhlen(material, bg2, cutEnergy);
 
    } else {
 
      //G4double dedx1 = dedxlim*betalow*material->GetDensity();
      G4double dedx1 = (*dedx0)[CurrentCouple()->GetIndex()] * betalow;
      G4double dedx2 = ComputeDEDXAhlen(material, bg2lim, cutEnergy);
 
      // extrapolation between two formula
      G4double kapa2 = beta - betalow;
      G4double kapa1 = betalim - beta;
      dedx = (kapa1 * dedx1 + kapa2 * dedx2) / (kapa1 + kapa2);
    }
  }
  return dedx;
}

◆ Dispersion()

G4double Dispersion	(	const G4Material *	material,
		const G4DynamicParticle *	dp,
		G4double	tcut,
		G4double	tmax,
		G4double	length
	)

overridevirtual

Calculate dispersion.

Parameters

material	Pointer to current material definition
dp	Secondary particle pointer
tcut	Tcut in dispersion formula
tmax	Tmax in dispersion formula?
length	Length in dispersion formula?

Returns: Dispersion

Definition at line 309 of file G4mplIonisationWithDeltaModel.cc.

{
  G4double siga = 0.0;
  G4double tau   = dp->GetKineticEnergy() / mass;
  if (tau > 0.0) {
    G4double electronDensity = material->GetElectronDensity();
    const G4double beta = dp->GetBeta();
    siga  = (tmax / (beta * beta) - 0.5 * tcut) * twopi_mc2_rcl2 * length
            * electronDensity * chargeSquare;
  }
  return siga;
}

◆ Initialise()

void Initialise	(	const G4ParticleDefinition *	p,
		const G4DataVector &
	)

overridevirtual

G4VEmModel::Initialise implementation.

Parameters

p	Monopole definition pointer.

Definition at line 84 of file G4mplIonisationWithDeltaModel.cc.

{
  if (!monopole) { SetParticle(p); }
  if (!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
  if (IsMaster()) {
    if (!dedx0) { dedx0 = new std::vector<G4double>; }
    G4ProductionCutsTable* theCoupleTable =
      G4ProductionCutsTable::GetProductionCutsTable();
    G4int numOfCouples = theCoupleTable->GetTableSize();
    G4int n = dedx0->size();
    if (n < numOfCouples) { dedx0->resize(numOfCouples); }
    G4Pow* g4calc = G4Pow::GetInstance();
 
    // initialise vector
    for (G4int i = 0; i < numOfCouples; ++i) {
 
      const G4Material* material =
        theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
      G4double eDensity = material->GetElectronDensity();
      G4double vF = electron_Compton_length * g4calc->A13(3.*pi * pi * eDensity);
      (*dedx0)[i] = pi_hbarc2_over_mc2 * eDensity * chargeSquare *
                    (G4Log(2 * vF / fine_structure_const) - 0.5) / vF;
    }
  }
}

◆ MaxSecondaryEnergy()

G4double MaxSecondaryEnergy	(	const G4ParticleDefinition *	p,
		G4double	kinEnergy
	)

overrideprotectedvirtual

Calculate maximum energy available for secondary particle emission.

Parameters

p	This argument of base function is ignored
kinEnergy	Kinetic energy of the ionising particle

Returns: Maximum energy available for secondary particle emission

Definition at line 327 of file G4mplIonisationWithDeltaModel.cc.

{
  G4double tau = kinEnergy / mass;
  return 2.0 * electron_mass_c2 * tau * (tau + 2.);
}

◆ MinEnergyCut()

G4double MinEnergyCut	(	const G4ParticleDefinition *	,
		const G4MaterialCutsCouple *	couple
	)

overridevirtual

Threshold for zero value.

Definition at line 112 of file G4mplIonisationWithDeltaModel.cc.

{
  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
}

◆ operator=()

G4mplIonisationWithDeltaModel & operator= ( const G4mplIonisationWithDeltaModel & right )

privatedelete

Assignment operator should be hidden.

Parameters

right Assign reference

◆ SampleFluctuations()

G4double SampleFluctuations	(	const G4MaterialCutsCouple *	couple,
		const G4DynamicParticle *	dp,
		G4double	tcut,
		G4double	tmax,
		G4double	length,
		G4double	meanLoss
	)

overridevirtual

Create fluctuations in the energies lost to a secondary delta electron.

Parameters

couple	Current volume couple to get its material
dp	Secondary particle pointer
tcut	Tcut in dispersion formula
tmax	Tmax in dispersion formula?
length	Length in dispersion formula?
meanLoss	Mean of energy loss

Returns: Fluctuated loss

Definition at line 279 of file G4mplIonisationWithDeltaModel.cc.

{
  G4double siga = Dispersion(couple->GetMaterial(), dp, tcut, tmax, length);
  G4double loss = meanLoss;
  siga = sqrt(siga);
  G4double twomeanLoss = meanLoss + meanLoss;
 
  if (twomeanLoss < siga) {
    G4double x;
    do {
      loss = twomeanLoss * G4UniformRand();
      x = (loss - meanLoss) / siga;
      // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
    } while (1.0 - 0.5 * x * x < G4UniformRand());
  } else {
    do {
      loss = G4RandGauss::shoot(meanLoss, siga);
      // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
    } while (0.0 > loss || loss > twomeanLoss);
  }
  return loss;
}

◆ SampleSecondaries()

void SampleSecondaries	(	std::vector< G4DynamicParticle * > *	vdp,
		const G4MaterialCutsCouple *	materialCutsCouple,
		const G4DynamicParticle *	dp,
		G4double	tmin,
		G4double	maxEnergy
	)

overridevirtual

Create the sample of secondary delta electrons.

Parameters

vdp	Pointer to the storage of sampled delta electrons
materialCutsCouple	This argument of base function is ignored
dp	Secondary particle pointer
tmin	Minimal energy of secondary electrons
maxEnergy	Maximum energy of secondary electrons

Definition at line 227 of file G4mplIonisationWithDeltaModel.cc.

{
  G4double kineticEnergy = dp->GetKineticEnergy();
  G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(), kineticEnergy);
 
  G4double maxKinEnergy = std::min(maxEnergy, tmax);
  if (minKinEnergy >= maxKinEnergy) { return; }
 
  G4double totEnergy     = kineticEnergy + mass;
  G4double etot2         = totEnergy * totEnergy;
  G4double beta2         = kineticEnergy * (kineticEnergy + 2.0 * mass) / etot2;
 
  // sampling without nuclear size effect
  G4double q = G4UniformRand();
  G4double deltaKinEnergy = minKinEnergy * maxKinEnergy
                            / (minKinEnergy * (1.0 - q) + maxKinEnergy * q);
 
  // delta-electron is produced
  G4double totMomentum = totEnergy * sqrt(beta2);
  G4double deltaMomentum =
    sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0 * electron_mass_c2));
  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
                  (deltaMomentum * totMomentum);
  cost = std::min(cost, 1.0);
 
  G4double sint = sqrt((1.0 - cost) * (1.0 + cost));
 
  G4double phi = twopi * G4UniformRand() ;
 
  G4ThreeVector deltaDirection(sint * cos(phi), sint * sin(phi), cost);
  G4ThreeVector direction = dp->GetMomentumDirection();
  deltaDirection.rotateUz(direction);
 
  // create G4DynamicParticle object for delta ray
  G4DynamicParticle* delta =
    new G4DynamicParticle(theElectron, deltaDirection, deltaKinEnergy);
 
  vdp->push_back(delta);
 
  // Change kinematics of primary particle
  kineticEnergy       -= deltaKinEnergy;
  G4ThreeVector finalP = direction * totMomentum - deltaDirection * deltaMomentum;
  finalP               = finalP.unit();
 
  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
  fParticleChange->SetProposedMomentumDirection(finalP);
}

◆ SetParticle()

void SetParticle ( const G4ParticleDefinition * p )

Read definition of the monopole.

Definition at line 71 of file G4mplIonisationWithDeltaModel.cc.

{
  monopole = p;
  mass     = monopole->GetPDGMass();
  G4double emin =
    std::min(LowEnergyLimit(), 0.1 * mass * (1. / sqrt(1. - betalow * betalow) - 1.));
  G4double emax =
    std::max(HighEnergyLimit(), 10 * mass * (1. / sqrt(1. - beta2lim) - 1.));
  SetLowEnergyLimit(emin);
  SetHighEnergyLimit(emax);
}

Member Data Documentation

◆ beta2lim

G4double beta2lim

private

Square of betalim.

Definition at line 222 of file G4mplIonisationWithDeltaModel.h.

◆ betalim

G4double betalim

private

Beta threshold for high energy (only Ahlen formula)

Definition at line 221 of file G4mplIonisationWithDeltaModel.h.

◆ betalow

G4double betalow

private

Beta threshold for low asymptotic.

Definition at line 220 of file G4mplIonisationWithDeltaModel.h.

◆ bg2lim

G4double bg2lim

private

(beta*gamma)^2 for betalim

Definition at line 223 of file G4mplIonisationWithDeltaModel.h.

◆ chargeSquare

G4double chargeSquare

private

Square of magnetic charge in units of Dirac charge.

Definition at line 224 of file G4mplIonisationWithDeltaModel.h.

◆ dedx0

std::vector< G4double > * dedx0 = nullptr

staticprivate

Base dedx for each couple in current material.

Definition at line 229 of file G4mplIonisationWithDeltaModel.h.

◆ dedxlim

G4double dedxlim

private

dedx limit in asymptotic formula, not used

Definition at line 225 of file G4mplIonisationWithDeltaModel.h.

◆ fParticleChange

G4ParticleChangeForLoss* fParticleChange

private

Pointer to ionising particle.

Definition at line 214 of file G4mplIonisationWithDeltaModel.h.

◆ magCharge

G4double magCharge

private

Monopole magnetic charge in e+ units.

Definition at line 217 of file G4mplIonisationWithDeltaModel.h.

◆ mass

G4double mass

private

Mass of the monopole.

Definition at line 216 of file G4mplIonisationWithDeltaModel.h.

◆ monopole

const G4ParticleDefinition* monopole

private

Monopole definition.

Definition at line 212 of file G4mplIonisationWithDeltaModel.h.

◆ nmpl

G4double nmpl

private

Magn.

charge in units of Dirac charge g_D, nmpl=1 -> g = g_D

Definition at line 227 of file G4mplIonisationWithDeltaModel.h.

◆ pi_hbarc2_over_mc2

G4double pi_hbarc2_over_mc2

private

Convenient constants combination with mass.

Definition at line 226 of file G4mplIonisationWithDeltaModel.h.

◆ theElectron

G4ParticleDefinition* theElectron

private

Electron definition.

Definition at line 213 of file G4mplIonisationWithDeltaModel.h.

◆ twoln10

G4double twoln10

private

log(100.0)

Definition at line 219 of file G4mplIonisationWithDeltaModel.h.

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

simulation/monopoles/include/G4mplIonisationWithDeltaModel.h
simulation/monopoles/src/G4mplIonisationWithDeltaModel.cc

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ G4mplIonisationWithDeltaModel() [1/2]

◆ ~G4mplIonisationWithDeltaModel()

◆ G4mplIonisationWithDeltaModel() [2/2]

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

◆ ComputeCrossSectionPerElectron()

◆ ComputeDEDXAhlen()

◆ ComputeDEDXPerVolume()

◆ Dispersion()

◆ Initialise()

◆ MaxSecondaryEnergy()

◆ MinEnergyCut()

◆ operator=()

◆ SampleFluctuations()

◆ SampleSecondaries()

◆ SetParticle()

Member Data Documentation

◆ beta2lim

◆ betalim

◆ betalow

◆ bg2lim

◆ chargeSquare

◆ dedx0

◆ dedxlim

◆ fParticleChange

◆ magCharge

◆ mass

◆ monopole

◆ nmpl

◆ pi_hbarc2_over_mc2

◆ theElectron

◆ twoln10