51 #include "G4mplAtlasIonisationModel.hh"
53 #include "Randomize.hh"
54 #include "G4LossTableManager.hh"
55 #include "G4ParticleChangeForLoss.hh"
56 #include "G4Version.hh"
58 #include "CLHEP/Units/SystemOfUnits.h"
59 #include "CLHEP/Units/PhysicalConstants.h"
65 G4mplAtlasIonisationModel::G4mplAtlasIonisationModel(G4double mCharge,
67 : G4VEmModel(nam),G4VEmFluctuationModel(nam),
72 twoln10(2.0*
log(10.0)),
75 bg2lim(beta2lim*(1.0 + beta2lim))
77 std::cout <<
"!!! G4mplAtlasIonisationModel constructor"<<std::endl;
80 nmpl = G4int(abs(magCharge)/68.0);
81 if(nmpl > 6) nmpl = 6;
82 else if(nmpl < 1) nmpl = 1;
85 chargeSquare = magCharge*magCharge;
90 G4mplAtlasIonisationModel::~G4mplAtlasIonisationModel()
92 std::cout <<
"!!! G4mplAtlasIonisationModel destructor"<<std::endl;
97 void G4mplAtlasIonisationModel::Initialise(
const G4ParticleDefinition*
p,
101 std::cout <<
"!!! G4mplAtlasIonisationModel::Initialise"<<std::endl;
104 mass = monopole->GetPDGMass();
107 fParticleChange =
reinterpret_cast<G4ParticleChangeForLoss*
>(pParticleChange);
109 fParticleChange =
new G4ParticleChangeForLoss();
114 G4double G4mplAtlasIonisationModel::ComputeDEDXPerVolume(
const G4Material* material,
115 const G4ParticleDefinition*,
116 G4double kineticEnergy,
119 G4double tau = kineticEnergy/
mass;
120 G4double gam = tau + 1.0;
121 G4double bg2 = tau * (tau+2.0);
122 G4double beta2 = bg2/(gam*gam);
124 G4double dedx0 = factlow*abs(beta2);
126 if(beta2 > beta2low) {
129 if(beta2 < beta2lim) {
134 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
135 G4double cden = material->GetIonisation()->GetCdensity();
136 G4double mden = material->GetIonisation()->GetMdensity();
137 G4double aden = material->GetIonisation()->GetAdensity();
138 G4double x0den = material->GetIonisation()->GetX0density();
139 G4double x1den = material->GetIonisation()->GetX1density();
143 G4double eDensity = material->GetElectronDensity();
148 if(nmpl > 1)
k = 0.346;
149 const G4double
B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};
154 G4double
x =
log(bg2)/twoln10;
156 dedx -= twoln10*
x - cden ;
157 if (
x < x1den ) dedx -= aden*
pow((x1den-
x),mden) ;
162 if (dedx < 0.0) dedx = 0.0 ;
167 if(beta2 < beta2lim) {
168 x =
log(dedx0) +
log(dedx/dedx0)*
log(b2/beta2low)/
log(beta2lim/beta2low);
179 G4double G4mplAtlasIonisationModel::SampleFluctuations(
180 #
if G4VERSION_NUMBER < 1100
181 const G4MaterialCutsCouple* material,
182 const G4DynamicParticle*
dp,
187 const G4MaterialCutsCouple* material,
188 const G4DynamicParticle*
dp,
192 const G4double meanLoss)
196 #if G4VERSION_NUMBER < 1100
197 G4double siga = Dispersion(material->GetMaterial(),
dp,tmax,
length);
199 G4double siga = Dispersion(material->GetMaterial(),
dp,tcut,tmax,
length);
201 G4double loss = meanLoss;
203 G4double twomeanLoss = meanLoss + meanLoss;
205 if(twomeanLoss < siga) {
208 loss = twomeanLoss*G4UniformRand();
209 x = (loss - meanLoss)/siga;
210 }
while (1.0 - 0.5*
x*
x < G4UniformRand());
213 loss = G4RandGauss::shoot(meanLoss,siga);
214 }
while (0.0 > loss || loss > twomeanLoss);