Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4TablesForExtrapolator.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//---------------------------------------------------------------------------
27//
28// ClassName: G4TablesForExtrapolator
29//
30// Description: This class provide calculation of energy loss, fluctuation,
31// and msc angle
32//
33// Author: 09.12.04 V.Ivanchenko
34//
35// Modification:
36// 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
37// 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
38// 21-03-06 Add verbosity defined in the constructor and Initialisation
39// start only when first public method is called (V.Ivanchenko)
40// 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
41// 12-05-06 SEt linLossLimit=0.001 (VI)
42//
43//----------------------------------------------------------------------------
44//
45
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48#include "G4TablesForExtrapolator.hh"
49#include "G4PhysicalConstants.hh"
50#include "G4SystemOfUnits.hh"
51#include "G4PhysicsLogVector.hh"
52#include "G4ParticleDefinition.hh"
53#include "G4Material.hh"
54#include "G4MaterialCutsCouple.hh"
55#include "G4Electron.hh"
56#include "G4Positron.hh"
57#include "G4Proton.hh"
58#include "G4MuonPlus.hh"
59#include "G4MuonMinus.hh"
60#include "G4EmParameters.hh"
61#include "G4MollerBhabhaModel.hh"
62#include "G4BetheBlochModel.hh"
63#include "G4eBremsstrahlungRelModel.hh"
64#include "G4MuPairProductionModel.hh"
65#include "G4MuBremsstrahlungModel.hh"
66#include "G4ProductionCuts.hh"
67#include "G4LossTableBuilder.hh"
68#include "G4WentzelVIModel.hh"
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4TablesForExtrapolator::G4TablesForExtrapolator(G4int verb, G4int bins,
73 G4double e1, G4double e2)
74 : builder(nullptr), verbose(verb), nbins(bins), nmat(0), emin(e1), emax(e2)
75{
76 electron = G4Electron::Electron();
77 positron = G4Positron::Positron();
78 proton = G4Proton::Proton();
79 muonPlus = G4MuonPlus::MuonPlus();
80 muonMinus= G4MuonMinus::MuonMinus();
81 dedxElectron = dedxPositron = dedxProton = dedxMuon = rangeElectron =
82 rangePositron = rangeProton = rangeMuon = invRangeElectron =
83 invRangePositron = invRangeProton = invRangeMuon = mscElectron = nullptr;
84 pcuts = nullptr;
85 Initialisation();
86}
87
88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89
90G4TablesForExtrapolator::~G4TablesForExtrapolator()
91{
92 for(G4int i=0; i<nmat; i++) { delete couples[i]; }
93
94 dedxElectron->clearAndDestroy();
95 dedxPositron->clearAndDestroy();
96 dedxProton->clearAndDestroy();
97 dedxMuon->clearAndDestroy();
98 rangeElectron->clearAndDestroy();
99 rangePositron->clearAndDestroy();
100 rangeProton->clearAndDestroy();
101 rangeMuon->clearAndDestroy();
102 invRangeElectron->clearAndDestroy();
103 invRangePositron->clearAndDestroy();
104 invRangeProton->clearAndDestroy();
105 invRangeMuon->clearAndDestroy();
106 mscElectron->clearAndDestroy();
107
108 delete dedxElectron;
109 delete dedxPositron;
110 delete dedxProton;
111 delete dedxMuon;
112 delete rangeElectron;
113 delete rangePositron;
114 delete rangeProton;
115 delete rangeMuon;
116 delete invRangeElectron;
117 delete invRangePositron;
118 delete invRangeProton;
119 delete invRangeMuon;
120 delete mscElectron;
121 delete pcuts;
122 delete builder;
123}
124
125//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
126
127const G4PhysicsTable*
128G4TablesForExtrapolator::GetPhysicsTable(ExtTableType type) const
129{
130 const G4PhysicsTable* table = nullptr;
131 switch (type)
132 {
133 case fDedxElectron:
134 table = dedxElectron;
135 break;
136 case fDedxPositron:
137 table = dedxPositron;
138 break;
139 case fDedxProton:
140 table = dedxProton;
141 break;
142 case fDedxMuon:
143 table = dedxMuon;
144 break;
145 case fRangeElectron:
146 table = rangeElectron;
147 break;
148 case fRangePositron:
149 table = rangePositron;
150 break;
151 case fRangeProton:
152 table = rangeProton;
153 break;
154 case fRangeMuon:
155 table = rangeMuon;
156 break;
157 case fInvRangeElectron:
158 table = invRangeElectron;
159 break;
160 case fInvRangePositron:
161 table = invRangePositron;
162 break;
163 case fInvRangeProton:
164 table = invRangeProton;
165 break;
166 case fInvRangeMuon:
167 table = invRangeMuon;
168 break;
169 case fMscElectron:
170 table = mscElectron;
171 }
172 return table;
173}
174
175//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
176
177void G4TablesForExtrapolator::Initialisation()
178{
179 if(verbose>1) {
180 G4cout << "### G4TablesForExtrapolator::Initialisation" << G4endl;
181 }
182 currentParticle = nullptr;
183 mass = charge2 = 0.0;
184
185 nmat = G4Material::GetNumberOfMaterials();
186 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
187 if(!pcuts) { pcuts = new G4ProductionCuts(); }
188
189 G4int i0 = couples.size();
190 if(0 == i0) {
191 couples.reserve(nmat);
192 cuts.reserve(nmat);
193 }
194 for(G4int i=i0; i<nmat; ++i) {
195 couples.push_back(new G4MaterialCutsCouple((*mtable)[i],pcuts));
196 cuts.push_back(DBL_MAX);
197 }
198
199 splineFlag = G4EmParameters::Instance()->Spline();
200
201 dedxElectron = PrepareTable(dedxElectron);
202 dedxPositron = PrepareTable(dedxPositron);
203 dedxMuon = PrepareTable(dedxMuon);
204 dedxProton = PrepareTable(dedxProton);
205 rangeElectron = PrepareTable(rangeElectron);
206 rangePositron = PrepareTable(rangePositron);
207 rangeMuon = PrepareTable(rangeMuon);
208 rangeProton = PrepareTable(rangeProton);
209 invRangeElectron = PrepareTable(invRangeElectron);
210 invRangePositron = PrepareTable(invRangePositron);
211 invRangeMuon = PrepareTable(invRangeMuon);
212 invRangeProton = PrepareTable(invRangeProton);
213 mscElectron = PrepareTable(mscElectron);
214
215 if(!builder) { builder = new G4LossTableBuilder(); }
216 builder->InitialiseBaseMaterials();
217
218 if(verbose>1) {
219 G4cout << "### G4TablesForExtrapolator Builds electron tables"
220 << G4endl;
221 }
222 ComputeElectronDEDX(electron, dedxElectron);
223 builder->BuildRangeTable(dedxElectron,rangeElectron);
224 builder->BuildInverseRangeTable(rangeElectron, invRangeElectron);
225
226 if(verbose>1) {
227 G4cout << "### G4TablesForExtrapolator Builds positron tables"
228 << G4endl;
229 }
230 ComputeElectronDEDX(positron, dedxPositron);
231 builder->BuildRangeTable(dedxPositron, rangePositron);
232 builder->BuildInverseRangeTable(rangePositron, invRangePositron);
233
234 if(verbose>1) {
235 G4cout << "### G4TablesForExtrapolator Builds muon tables" << G4endl;
236 }
237 ComputeMuonDEDX(muonPlus, dedxMuon);
238 builder->BuildRangeTable(dedxMuon, rangeMuon);
239 builder->BuildInverseRangeTable(rangeMuon, invRangeMuon);
240 /*
241 G4cout << "DEDX MUON" << G4endl
242 G4cout << *dedxMuon << G4endl;
243 G4cout << "RANGE MUON" << G4endl
244 G4cout << *rangeMuon << G4endl;
245 G4cout << "INVRANGE MUON" << G4endl
246 G4cout << *invRangeMuon << G4endl;
247 */
248 if(verbose>1) {
249 G4cout << "### G4TablesForExtrapolator Builds proton tables"
250 << G4endl;
251 }
252 ComputeProtonDEDX(proton, dedxProton);
253 builder->BuildRangeTable(dedxProton, rangeProton);
254 builder->BuildInverseRangeTable(rangeProton, invRangeProton);
255
256 ComputeTrasportXS(electron, mscElectron);
257}
258
259//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260
261G4PhysicsTable* G4TablesForExtrapolator::PrepareTable(G4PhysicsTable* ptr)
262{
263 G4PhysicsTable* table;
264 G4int i0 = 0;
265 if(ptr) {
266 table = ptr;
267 i0 = ptr->size();
268 } else {
269 table = new G4PhysicsTable();
270 }
271 for(G4int i=i0; i<nmat; ++i) {
272 G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
273 v->SetSpline(splineFlag);
274 table->push_back(v);
275 }
276 return table;
277}
278
279//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
280
281void G4TablesForExtrapolator::ComputeElectronDEDX(
282 const G4ParticleDefinition* part,
283 G4PhysicsTable* table)
284{
285 G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel();
286 G4eBremsstrahlungRelModel* brem = new G4eBremsstrahlungRelModel();
287 ioni->Initialise(part, cuts);
288 brem->Initialise(part, cuts);
289 ioni->SetUseBaseMaterials(false);
290 brem->SetUseBaseMaterials(false);
291
292 mass = electron_mass_c2;
293 charge2 = 1.0;
294 currentParticle = part;
295
296 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
297 if(0<verbose) {
298 G4cout << "G4TablesForExtrapolator::ComputeElectronDEDX for "
299 << part->GetParticleName()
300 << G4endl;
301 }
302 for(G4int i=0; i<nmat; ++i) {
303
304 const G4Material* mat = (*mtable)[i];
305 if(1<verbose) {
306 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
307 }
308 const G4MaterialCutsCouple* couple = couples[i];
309 G4PhysicsVector* aVector = (*table)[i];
310
311 for(G4int j=0; j<=nbins; ++j) {
312
313 G4double e = aVector->Energy(j);
314 G4double dedx = ioni->ComputeDEDX(couple,part,e,e)
315 + brem->ComputeDEDX(couple,part,e,e);
316 if(1<verbose) {
317 G4cout << "j= " << j
318 << " e(MeV)= " << e/MeV
319 << " dedx(Mev/cm)= " << dedx*cm/MeV
320 << " dedx(Mev.cm2/g)= "
321 << dedx/((MeV*mat->GetDensity())/(g/cm2))
322 << G4endl;
323 }
324 aVector->PutValue(j,dedx);
325 }
326 if(splineFlag) { aVector->FillSecondDerivatives(); }
327 }
328 delete ioni;
329 delete brem;
330}
331
332//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333
334void
335G4TablesForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part,
336 G4PhysicsTable* table)
337{
338 G4BetheBlochModel* ioni = new G4BetheBlochModel();
339 G4MuPairProductionModel* pair = new G4MuPairProductionModel();
340 G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel();
341 ioni->Initialise(part, cuts);
342 pair->Initialise(part, cuts);
343 brem->Initialise(part, cuts);
344 ioni->SetUseBaseMaterials(false);
345 pair->SetUseBaseMaterials(false);
346 brem->SetUseBaseMaterials(false);
347
348 mass = part->GetPDGMass();
349 charge2 = 1.0;
350 currentParticle = part;
351
352 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
353
354 if(0<verbose) {
355 G4cout << "G4TablesForExtrapolator::ComputeMuonDEDX for "
356 << part->GetParticleName()
357 << G4endl;
358 }
359
360 for(G4int i=0; i<nmat; ++i) {
361
362 const G4Material* mat = (*mtable)[i];
363 if(1<verbose) {
364 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
365 }
366 const G4MaterialCutsCouple* couple = couples[i];
367 G4PhysicsVector* aVector = (*table)[i];
368 for(G4int j=0; j<=nbins; j++) {
369
370 G4double e = aVector->Energy(j);
371 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
372 pair->ComputeDEDX(couple,part,e,e) +
373 brem->ComputeDEDX(couple,part,e,e);
374 aVector->PutValue(j,dedx);
375 if(1<verbose) {
376 G4cout << "j= " << j
377 << " e(MeV)= " << e/MeV
378 << " dedx(Mev/cm)= " << dedx*cm/MeV
379 << " dedx(Mev/(g/cm2)= "
380 << dedx/((MeV*mat->GetDensity())/(g/cm2))
381 << G4endl;
382 }
383 }
384 if(splineFlag) { aVector->FillSecondDerivatives(); }
385 }
386 delete ioni;
387}
388
389//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390
391void
392G4TablesForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part,
393 G4PhysicsTable* table)
394{
395 G4BetheBlochModel* ioni = new G4BetheBlochModel();
396 ioni->Initialise(part, cuts);
397 ioni->SetUseBaseMaterials(false);
398
399 mass = part->GetPDGMass();
400 charge2 = 1.0;
401 currentParticle = part;
402
403 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
404
405 if(0<verbose) {
406 G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for "
407 << part->GetParticleName()
408 << G4endl;
409 }
410
411 for(G4int i=0; i<nmat; ++i) {
412
413 const G4Material* mat = (*mtable)[i];
414 if(1<verbose)
415 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
416 const G4MaterialCutsCouple* couple = couples[i];
417 G4PhysicsVector* aVector = (*table)[i];
418 for(G4int j=0; j<=nbins; j++) {
419
420 G4double e = aVector->Energy(j);
421 G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
422 aVector->PutValue(j,dedx);
423 if(1<verbose) {
424 G4cout << "j= " << j
425 << " e(MeV)= " << e/MeV
426 << " dedx(Mev/cm)= " << dedx*cm/MeV
427 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2))
428 << G4endl;
429 }
430 }
431 if(splineFlag) { aVector->FillSecondDerivatives(); }
432 }
433 delete ioni;
434}
435
436//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
437
438void
439G4TablesForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part,
440 G4PhysicsTable* table)
441{
442 G4WentzelVIModel* msc = new G4WentzelVIModel();
443 msc->SetPolarAngleLimit(CLHEP::pi);
444 msc->Initialise(part, cuts);
445 msc->SetUseBaseMaterials(false);
446
447 mass = part->GetPDGMass();
448 charge2 = 1.0;
449 currentParticle = part;
450
451 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
452
453 if(0<verbose) {
454 G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for "
455 << part->GetParticleName()
456 << G4endl;
457 }
458
459 for(G4int i=0; i<nmat; ++i) {
460
461 const G4Material* mat = (*mtable)[i];
462 msc->SetCurrentCouple(couples[i]);
463 if(1<verbose)
464 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
465 G4PhysicsVector* aVector = (*table)[i];
466 for(G4int j=0; j<=nbins; j++) {
467
468 G4double e = aVector->Energy(j);
469 G4double xs = msc->CrossSectionPerVolume(mat,part,e);
470 aVector->PutValue(j,xs);
471 if(1<verbose) {
472 G4cout << "j= " << j << " e(MeV)= " << e/MeV
473 << " xs(1/mm)= " << xs*mm << G4endl;
474 }
475 }
476 if(splineFlag) { aVector->FillSecondDerivatives(); }
477 }
478 delete msc;
479}
480
481//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
482
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:44
fInvRangeElectron
@ fInvRangeElectron
Definition: G4TablesForExtrapolator.hh:67
fInvRangePositron
@ fInvRangePositron
Definition: G4TablesForExtrapolator.hh:68
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4BetheBlochModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4BetheBlochModel.cc:103
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition: G4EmParameters.cc:67
G4EmParameters::Spline
G4bool Spline() const
Definition: G4EmParameters.cc:215
G4LossTableBuilder::BuildInverseRangeTable
void BuildInverseRangeTable(const G4PhysicsTable *rangeTable, G4PhysicsTable *invRangeTable, G4bool useBM=false)
Definition: G4LossTableBuilder.cc:259
G4LossTableBuilder::BuildRangeTable
void BuildRangeTable(const G4PhysicsTable *dedxTable, G4PhysicsTable *rangeTable, G4bool useBM=false)
Definition: G4LossTableBuilder.cc:174
G4LossTableBuilder::InitialiseBaseMaterials
void InitialiseBaseMaterials(const G4PhysicsTable *table=nullptr)
Definition: G4LossTableBuilder.cc:291
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:178
G4Material::GetNumberOfMaterials
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:644
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:637
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:175
G4MollerBhabhaModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4MollerBhabhaModel.cc:106
G4MuBremsstrahlungModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4MuBremsstrahlungModel.cc:144
G4MuPairProductionModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4MuPairProductionModel.cc:163
G4MuonMinus::MuonMinus
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:99
G4MuonPlus::MuonPlus
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:98
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:108
G4PhysicsTable::push_back
void push_back(G4PhysicsVector *)
G4PhysicsTable::clearAndDestroy
void clearAndDestroy()
G4PhysicsVector::Energy
G4double Energy(std::size_t index) const
G4PhysicsVector::PutValue
void PutValue(std::size_t index, G4double theValue)
G4PhysicsVector::FillSecondDerivatives
void FillSecondDerivatives()
Definition: G4PhysicsVector.cc:361
G4PhysicsVector::SetSpline
void SetSpline(G4bool)
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:93
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:92
G4TablesForExtrapolator::G4TablesForExtrapolator
G4TablesForExtrapolator(G4int verb, G4int bins, G4double e1, G4double e2)
Definition: G4TablesForExtrapolator.cc:72
G4TablesForExtrapolator::GetPhysicsTable
const G4PhysicsTable * GetPhysicsTable(ExtTableType type) const
Definition: G4TablesForExtrapolator.cc:128
G4VEmModel::SetPolarAngleLimit
void SetPolarAngleLimit(G4double)
Definition: G4VEmModel.hh:792
G4VEmModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:254
G4VEmModel::SetCurrentCouple
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:465
G4VEmModel::ComputeDEDX
virtual G4double ComputeDEDX(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
Definition: G4VEmModel.hh:518
G4VEmModel::SetUseBaseMaterials
void SetUseBaseMaterials(G4bool val)
Definition: G4VEmModel.hh:743
G4WentzelVIModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4WentzelVIModel.cc:122
G4eBremsstrahlungRelModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4eBremsstrahlungRelModel.cc:170
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62