Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // $Id: G4TablesForExtrapolator.cc 75145 2013-10-28 18:34:48Z vnivanch $
27 //
28 //---------------------------------------------------------------------------
29 //
30 // ClassName: G4TablesForExtrapolator
31 //
32 // Description: This class provide calculation of energy loss, fluctuation,
33 // and msc angle
34 //
35 // Author: 09.12.04 V.Ivanchenko
36 //
37 // Modification:
38 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
39 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
40 // 21-03-06 Add verbosity defined in the constructor and Initialisation
41 // start only when first public method is called (V.Ivanchenko)
42 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
43 // 12-05-06 SEt linLossLimit=0.001 (VI)
44 //
45 //----------------------------------------------------------------------------
46 //
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 
51 #include "G4PhysicalConstants.hh"
52 #include "G4SystemOfUnits.hh"
53 #include "G4PhysicsLogVector.hh"
54 #include "G4ParticleDefinition.hh"
55 #include "G4Material.hh"
56 #include "G4MaterialCutsCouple.hh"
57 #include "G4Electron.hh"
58 #include "G4Positron.hh"
59 #include "G4Proton.hh"
60 #include "G4MuonPlus.hh"
61 #include "G4MuonMinus.hh"
62 #include "G4EmParameters.hh"
63 #include "G4MollerBhabhaModel.hh"
64 #include "G4BetheBlochModel.hh"
68 #include "G4ProductionCuts.hh"
69 #include "G4LossTableBuilder.hh"
70 #include "G4WentzelVIModel.hh"
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
75  G4double e1, G4double e2)
76  : verbose(verb), nbins(bins), emin(e1), emax(e2)
77 {
78  Initialisation();
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82 
84 {
85  for(G4int i=0; i<nmat; i++) {delete couples[i];}
86 
87  dedxElectron->clearAndDestroy();
88  dedxPositron->clearAndDestroy();
89  dedxProton->clearAndDestroy();
90  dedxMuon->clearAndDestroy();
91  rangeElectron->clearAndDestroy();
92  rangePositron->clearAndDestroy();
93  rangeProton->clearAndDestroy();
94  rangeMuon->clearAndDestroy();
95  invRangeElectron->clearAndDestroy();
96  invRangePositron->clearAndDestroy();
97  invRangeProton->clearAndDestroy();
98  invRangeMuon->clearAndDestroy();
99  mscElectron->clearAndDestroy();
100 
101  delete dedxElectron;
102  delete dedxPositron;
103  delete dedxProton;
104  delete dedxMuon;
105  delete rangeElectron;
106  delete rangePositron;
107  delete rangeProton;
108  delete rangeMuon;
109  delete invRangeElectron;
110  delete invRangePositron;
111  delete invRangeProton;
112  delete invRangeMuon;
113  delete mscElectron;
114  delete pcuts;
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 
119 const G4PhysicsTable*
121 {
122  const G4PhysicsTable* table = nullptr;
123  switch (type)
124  {
125  case fDedxElectron:
126  table = dedxElectron;
127  break;
128  case fDedxPositron:
129  table = dedxPositron;
130  break;
131  case fDedxProton:
132  table = dedxProton;
133  break;
134  case fDedxMuon:
135  table = dedxMuon;
136  break;
137  case fRangeElectron:
138  table = rangeElectron;
139  break;
140  case fRangePositron:
141  table = rangePositron;
142  break;
143  case fRangeProton:
144  table = rangeProton;
145  break;
146  case fRangeMuon:
147  table = rangeMuon;
148  break;
149  case fInvRangeElectron:
150  table = invRangeElectron;
151  break;
152  case fInvRangePositron:
153  table = invRangePositron;
154  break;
155  case fInvRangeProton:
156  table = invRangeProton;
157  break;
158  case fInvRangeMuon:
159  table = invRangeMuon;
160  break;
161  case fMscElectron:
162  table = mscElectron;
163  }
164  return table;
165 }
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168 
169 void G4TablesForExtrapolator::Initialisation()
170 {
171  if(verbose>1) {
172  G4cout << "### G4TablesForExtrapolator::Initialisation" << G4endl;
173  }
174  currentParticle = nullptr;
175 
176  electron = G4Electron::Electron();
177  positron = G4Positron::Positron();
178  proton = G4Proton::Proton();
179  muonPlus = G4MuonPlus::MuonPlus();
180  muonMinus= G4MuonMinus::MuonMinus();
181 
182  mass = charge2 = 0.0;
183 
186  pcuts = new G4ProductionCuts();
187 
188  couples.resize(nmat,0);
189  cuts.resize(nmat,DBL_MAX);
190 
191  for(G4int i=0; i<nmat; ++i) {
192  couples[i] = new G4MaterialCutsCouple((*mtable)[i],pcuts);
193  }
194 
195  splineFlag = G4EmParameters::Instance()->Spline();
196 
197  dedxElectron = PrepareTable();
198  dedxPositron = PrepareTable();
199  dedxMuon = PrepareTable();
200  dedxProton = PrepareTable();
201  rangeElectron = PrepareTable();
202  rangePositron = PrepareTable();
203  rangeMuon = PrepareTable();
204  rangeProton = PrepareTable();
205  invRangeElectron = PrepareTable();
206  invRangePositron = PrepareTable();
207  invRangeMuon = PrepareTable();
208  invRangeProton = PrepareTable();
209  mscElectron = PrepareTable();
210 
211  G4LossTableBuilder builder;
212 
213  if(verbose>1) {
214  G4cout << "### G4TablesForExtrapolator Builds electron tables"
215  << G4endl;
216  }
217  ComputeElectronDEDX(electron, dedxElectron);
218  builder.BuildRangeTable(dedxElectron,rangeElectron);
219  builder.BuildInverseRangeTable(rangeElectron, invRangeElectron);
220 
221  if(verbose>1) {
222  G4cout << "### G4TablesForExtrapolator Builds positron tables"
223  << G4endl;
224  }
225  ComputeElectronDEDX(positron, dedxPositron);
226  builder.BuildRangeTable(dedxPositron, rangePositron);
227  builder.BuildInverseRangeTable(rangePositron, invRangePositron);
228 
229  if(verbose>1) {
230  G4cout << "### G4TablesForExtrapolator Builds muon tables" << G4endl;
231  }
232  ComputeMuonDEDX(muonPlus, dedxMuon);
233  builder.BuildRangeTable(dedxMuon, rangeMuon);
234  builder.BuildInverseRangeTable(rangeMuon, invRangeMuon);
235  /*
236  G4cout << "DEDX MUON" << G4endl
237  G4cout << *dedxMuon << G4endl;
238  G4cout << "RANGE MUON" << G4endl
239  G4cout << *rangeMuon << G4endl;
240  G4cout << "INVRANGE MUON" << G4endl
241  G4cout << *invRangeMuon << G4endl;
242  */
243  if(verbose>1) {
244  G4cout << "### G4TablesForExtrapolator Builds proton tables"
245  << G4endl;
246  }
247  ComputeProtonDEDX(proton, dedxProton);
248  builder.BuildRangeTable(dedxProton, rangeProton);
249  builder.BuildInverseRangeTable(rangeProton, invRangeProton);
250 
251  ComputeTrasportXS(electron, mscElectron);
252 }
253 
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255 
256 G4PhysicsTable* G4TablesForExtrapolator::PrepareTable()
257 {
258  G4PhysicsTable* table = new G4PhysicsTable();
259 
260  for(G4int i=0; i<nmat; ++i) {
261 
262  G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
263  v->SetSpline(splineFlag);
264  table->push_back(v);
265  }
266  return table;
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270 
271 void G4TablesForExtrapolator::ComputeElectronDEDX(
272  const G4ParticleDefinition* part,
273  G4PhysicsTable* table)
274 {
277  ioni->Initialise(part, cuts);
278  brem->Initialise(part, cuts);
279 
280  mass = electron_mass_c2;
281  charge2 = 1.0;
282  currentParticle = part;
283 
285  if(0<verbose) {
286  G4cout << "G4TablesForExtrapolator::ComputeElectronDEDX for "
287  << part->GetParticleName()
288  << G4endl;
289  }
290  for(G4int i=0; i<nmat; ++i) {
291 
292  const G4Material* mat = (*mtable)[i];
293  if(1<verbose) {
294  G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
295  }
296  const G4MaterialCutsCouple* couple = couples[i];
297  G4PhysicsVector* aVector = (*table)[i];
298 
299  for(G4int j=0; j<=nbins; ++j) {
300 
301  G4double e = aVector->Energy(j);
302  G4double dedx = ioni->ComputeDEDX(couple,part,e,e)
303  + brem->ComputeDEDX(couple,part,e,e);
304  if(1<verbose) {
305  G4cout << "j= " << j
306  << " e(MeV)= " << e/MeV
307  << " dedx(Mev/cm)= " << dedx*cm/MeV
308  << " dedx(Mev.cm2/g)= "
309  << dedx/((MeV*mat->GetDensity())/(g/cm2))
310  << G4endl;
311  }
312  aVector->PutValue(j,dedx);
313  }
314  if(splineFlag) { aVector->FillSecondDerivatives(); }
315  }
316  delete ioni;
317  delete brem;
318 }
319 
320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
321 
322 void
323 G4TablesForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part,
324  G4PhysicsTable* table)
325 {
326  G4BetheBlochModel* ioni = new G4BetheBlochModel();
329  ioni->Initialise(part, cuts);
330  pair->Initialise(part, cuts);
331  brem->Initialise(part, cuts);
332 
333  mass = part->GetPDGMass();
334  charge2 = 1.0;
335  currentParticle = part;
336 
338 
339  if(0<verbose) {
340  G4cout << "G4TablesForExtrapolator::ComputeMuonDEDX for "
341  << part->GetParticleName()
342  << G4endl;
343  }
344 
345  for(G4int i=0; i<nmat; ++i) {
346 
347  const G4Material* mat = (*mtable)[i];
348  if(1<verbose) {
349  G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
350  }
351  const G4MaterialCutsCouple* couple = couples[i];
352  G4PhysicsVector* aVector = (*table)[i];
353  for(G4int j=0; j<=nbins; j++) {
354 
355  G4double e = aVector->Energy(j);
356  G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
357  pair->ComputeDEDX(couple,part,e,e) +
358  brem->ComputeDEDX(couple,part,e,e);
359  aVector->PutValue(j,dedx);
360  if(1<verbose) {
361  G4cout << "j= " << j
362  << " e(MeV)= " << e/MeV
363  << " dedx(Mev/cm)= " << dedx*cm/MeV
364  << " dedx(Mev/(g/cm2)= "
365  << dedx/((MeV*mat->GetDensity())/(g/cm2))
366  << G4endl;
367  }
368  }
369  if(splineFlag) { aVector->FillSecondDerivatives(); }
370  }
371  delete ioni;
372 }
373 
374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
375 
376 void
377 G4TablesForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part,
378  G4PhysicsTable* table)
379 {
380  G4BetheBlochModel* ioni = new G4BetheBlochModel();
381  ioni->Initialise(part, cuts);
382 
383  mass = part->GetPDGMass();
384  charge2 = 1.0;
385  currentParticle = part;
386 
388 
389  if(0<verbose) {
390  G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for "
391  << part->GetParticleName()
392  << G4endl;
393  }
394 
395  for(G4int i=0; i<nmat; ++i) {
396 
397  const G4Material* mat = (*mtable)[i];
398  if(1<verbose)
399  G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
400  const G4MaterialCutsCouple* couple = couples[i];
401  G4PhysicsVector* aVector = (*table)[i];
402  for(G4int j=0; j<=nbins; j++) {
403 
404  G4double e = aVector->Energy(j);
405  G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
406  aVector->PutValue(j,dedx);
407  if(1<verbose) {
408  G4cout << "j= " << j
409  << " e(MeV)= " << e/MeV
410  << " dedx(Mev/cm)= " << dedx*cm/MeV
411  << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2))
412  << G4endl;
413  }
414  }
415  if(splineFlag) { aVector->FillSecondDerivatives(); }
416  }
417  delete ioni;
418 }
419 
420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
421 
422 void
423 G4TablesForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part,
424  G4PhysicsTable* table)
425 {
426  G4WentzelVIModel* msc = new G4WentzelVIModel();
428  msc->Initialise(part, cuts);
429 
430  mass = part->GetPDGMass();
431  charge2 = 1.0;
432  currentParticle = part;
433 
435 
436  if(0<verbose) {
437  G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for "
438  << part->GetParticleName()
439  << G4endl;
440  }
441 
442  for(G4int i=0; i<nmat; ++i) {
443 
444  const G4Material* mat = (*mtable)[i];
445  msc->SetCurrentCouple(couples[i]);
446  if(1<verbose)
447  G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
448  G4PhysicsVector* aVector = (*table)[i];
449  for(G4int j=0; j<=nbins; j++) {
450 
451  G4double e = aVector->Energy(j);
452  G4double xs = msc->CrossSectionPerVolume(mat,part,e);
453  aVector->PutValue(j,xs);
454  if(1<verbose) {
455  G4cout << "j= " << j << " e(MeV)= " << e/MeV
456  << " xs(1/mm)= " << xs*mm << G4endl;
457  }
458  }
459  if(splineFlag) { aVector->FillSecondDerivatives(); }
460  }
461  delete msc;
462 }
463 
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 
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)
Definition: G4VEmModel.cc:257
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:99
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4bool Spline() const
static constexpr double mm
Definition: G4SIunits.hh:115
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
static constexpr double cm2
Definition: G4SIunits.hh:120
const G4PhysicsTable * GetPhysicsTable(ExtTableType type) const
const G4String & GetName() const
Definition: G4Material.hh:178
void push_back(G4PhysicsVector *)
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
G4TablesForExtrapolator(G4int verb, G4int bins, G4double e1, G4double e2)
G4double GetDensity() const
Definition: G4Material.hh:180
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
int G4int
Definition: G4Types.hh:78
void FillSecondDerivatives()
const G4String & GetParticleName() const
void SetSpline(G4bool)
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
G4GLOB_DLL std::ostream G4cout
static constexpr double cm
Definition: G4SIunits.hh:119
void PutValue(size_t index, G4double theValue)
static G4Proton * Proton()
Definition: G4Proton.cc:93
float electron_mass_c2
Definition: hepunit.py:274
void BuildInverseRangeTable(const G4PhysicsTable *rangeTable, G4PhysicsTable *invRangeTable, G4bool isIonisation=false)
G4double Energy(size_t index) const
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
tuple v
Definition: test.py:18
static const G4double emax
static G4Positron * Positron()
Definition: G4Positron.cc:94
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
virtual G4double ComputeDEDX(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
Definition: G4VEmModel.hh:492
G4double GetPDGMass() const
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:447
static G4EmParameters * Instance()
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:100
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
static constexpr double pi
Definition: SystemOfUnits.h:54
void SetPolarAngleLimit(G4double)
Definition: G4VEmModel.hh:767
void clearAndDestroy()