76 :maxEnergyTransfer(
DBL_MAX),verbose(verb),isInitialised(false)
89 mass = charge2 = electronDensity = radLength = bg2 = beta2
90 = kineticEnergy = tmax = 0;
93 dedxElectron = dedxPositron = dedxProton = rangeElectron
94 = rangePositron = rangeProton = invRangeElectron = invRangePositron
95 = invRangeProton = mscElectron = dedxMuon = rangeMuon = invRangeMuon = 0;
97 electron = positron = proton = muonPlus = muonMinus = 0;
104 for(
G4int i=0; i<nmat; i++) {
delete couples[i];}
109 delete rangeElectron;
110 delete rangePositron;
113 delete invRangeElectron;
114 delete invRangePositron;
115 delete invRangeProton;
129 if(!isInitialised) Initialisation();
130 G4double kinEnergyFinal = kinEnergy;
131 if(SetupKinematics(part, mat, kinEnergy)) {
135 kinEnergyFinal = 0.0;
136 }
else if(step < linLossLimit*r) {
137 kinEnergyFinal -= step*
ComputeDEDX(kinEnergy,part);
143 return kinEnergyFinal;
154 if(!isInitialised) { Initialisation(); }
155 G4double kinEnergyFinal = kinEnergy;
157 if(SetupKinematics(part, mat, kinEnergy)) {
161 if(step < linLossLimit*r) {
162 kinEnergyFinal += step*
ComputeDEDX(kinEnergy,part);
168 return kinEnergyFinal;
180 if(!isInitialised) Initialisation();
181 if(SetupKinematics(part, mat, kinEnergy)) {
182 if(part == electron || part == positron) {
183 G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
184 if(x < 0.2) { res *= (1.0 + 0.5*x + x*x/3.0); }
185 else if(x < 0.9999) { res = -
G4Log(1.0 - x)*stepLength/
x; }
202 if(!part || !mat || kinEnergy <
keV)
return false;
203 if(!isInitialised) Initialisation();
205 if(part != currentParticle) {
207 currentParticle =
part;
212 if(mat != currentMaterial) {
215 G4cout <<
"### G4EnergyLossForExtrapolator WARNING:index i= "
216 << i <<
" is out of table - NO extrapolation" <<
G4endl;
219 currentMaterial =
mat;
225 if(flag || kinEnergy != kineticEnergy) {
226 kineticEnergy = kinEnergy;
230 bg2 = tau * (tau + 2.0);
231 beta2 = bg2/(gam*gam);
233 if(part == electron) tmax *= 0.5;
234 else if(part != positron) {
238 if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
245 void G4EnergyLossForExtrapolator::Initialisation()
247 isInitialised =
true;
249 G4cout <<
"### G4EnergyLossForExtrapolator::Initialisation" <<
G4endl;
261 currentParticleName =
"";
267 couples.resize(nmat,0);
270 for(
G4int i=0; i<nmat; i++) {
274 dedxElectron = PrepareTable();
275 dedxPositron = PrepareTable();
276 dedxMuon = PrepareTable();
277 dedxProton = PrepareTable();
278 rangeElectron = PrepareTable();
279 rangePositron = PrepareTable();
280 rangeMuon = PrepareTable();
281 rangeProton = PrepareTable();
282 invRangeElectron = PrepareTable();
283 invRangePositron = PrepareTable();
284 invRangeMuon = PrepareTable();
285 invRangeProton = PrepareTable();
286 mscElectron = PrepareTable();
291 G4cout <<
"### G4EnergyLossForExtrapolator Builds electron tables"
294 ComputeElectronDEDX(electron, dedxElectron);
299 G4cout <<
"### G4EnergyLossForExtrapolator Builds positron tables"
302 ComputeElectronDEDX(positron, dedxPositron);
307 G4cout <<
"### G4EnergyLossForExtrapolator Builds muon tables" <<
G4endl;
309 ComputeMuonDEDX(muonPlus, dedxMuon);
314 G4cout <<
"### G4EnergyLossForExtrapolator Builds proton tables"
317 ComputeProtonDEDX(proton, dedxProton);
321 ComputeTrasportXS(electron, mscElectron);
330 for(
G4int i=0; i<nmat; i++) {
342 G4EnergyLossForExtrapolator::FindParticle(
const G4String&
name)
345 if(name != currentParticleName) {
348 G4cout <<
"### G4EnergyLossForExtrapolator WARNING:FindParticle fails to find "
364 if(part == electron) x = ComputeValue(kinEnergy, dedxElectron);
365 else if(part == positron) x = ComputeValue(kinEnergy, dedxPositron);
366 else if(part == muonPlus || part == muonMinus) {
367 x = ComputeValue(kinEnergy, dedxMuon);
370 x = ComputeValue(e, dedxProton)*charge2;
382 if(part == electron) x = ComputeValue(kinEnergy, rangeElectron);
383 else if(part == positron) x = ComputeValue(kinEnergy, rangePositron);
384 else if(part == muonPlus || part == muonMinus)
385 x = ComputeValue(kinEnergy, rangeMuon);
389 x = ComputeValue(e, rangeProton)/(charge2*massratio);
401 if(part == electron) x = ComputeValue(range, invRangeElectron);
402 else if(part == positron) x = ComputeValue(range, invRangePositron);
403 else if(part == muonPlus || part == muonMinus)
404 x = ComputeValue(range, invRangeMuon);
408 x = ComputeValue(r, invRangeProton)/massratio;
415 void G4EnergyLossForExtrapolator::ComputeElectronDEDX(
426 currentParticle =
part;
430 G4cout <<
"G4EnergyLossForExtrapolator::ComputeElectronDEDX for "
434 for(
G4int i=0; i<nmat; i++) {
442 for(
G4int j=0; j<=nbins; j++) {
449 <<
" e(MeV)= " << e/
MeV
450 <<
" dedx(Mev/cm)= " << dedx*
cm/
MeV
451 <<
" dedx(Mev.cm2/g)= "
477 currentParticle =
part;
482 G4cout <<
"G4EnergyLossForExtrapolator::ComputeMuonDEDX for "
487 for(
G4int i=0; i<nmat; i++) {
495 for(
G4int j=0; j<=nbins; j++) {
504 <<
" e(MeV)= " << e/
MeV
505 <<
" dedx(Mev/cm)= " << dedx*
cm/
MeV
525 currentParticle =
part;
530 G4cout <<
"G4EnergyLossForExtrapolator::ComputeProtonDEDX for "
535 for(
G4int i=0; i<nmat; i++) {
542 for(
G4int j=0; j<=nbins; j++) {
549 <<
" e(MeV)= " << e/
MeV
550 <<
" dedx(Mev/cm)= " << dedx*
cm/
MeV
571 currentParticle =
part;
576 G4cout <<
"G4EnergyLossForExtrapolator::ComputeProtonDEDX for "
581 for(
G4int i=0; i<nmat; i++) {
588 for(
G4int j=0; j<=nbins; j++) {
594 G4cout <<
"j= " << j <<
" e(MeV)= " << e/
MeV
595 <<
" xs(1/mm)= " << xs*
mm <<
G4endl;
void BuildRangeTable(const G4PhysicsTable *dedxTable, G4PhysicsTable *rangeTable, G4bool isIonisation=false)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
static G4MuonPlus * MuonPlus()
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4LossTableManager * Instance()
const G4String & GetName() const
void push_back(G4PhysicsVector *)
static G4MaterialTable * GetMaterialTable()
std::vector< G4Material * > G4MaterialTable
G4double GetDensity() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition const G4Material *G4double range
G4double GetElectronDensity() const
void PutValue(size_t index, G4double theValue)
static G4Proton * Proton()
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void BuildInverseRangeTable(const G4PhysicsTable *rangeTable, G4PhysicsTable *invRangeTable, G4bool isIonisation=false)
G4double Energy(size_t index) const
static size_t GetNumberOfMaterials()
G4double GetRadlen() const
G4double G4Log(G4double x)
static G4Positron * Positron()
G4double ComputeDEDX(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
G4double GetPDGMass() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4ParticleTable * GetParticleTable()
void SetCurrentCouple(const G4MaterialCutsCouple *)
static G4MuonMinus * MuonMinus()
static G4Electron * Electron()
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double GetPDGCharge() const
void SetPolarAngleLimit(G4double)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)