Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VEmModel.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: G4VEmModel.cc 105120 2017-07-13 13:24:43Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4VEmModel
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 25.07.2005
38 //
39 // Modifications:
40 // 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko)
41 // 06.02.2006 add method ComputeMeanFreePath() (mma)
42 // 16.02.2009 Move implementations of virtual methods to source (VI)
43 //
44 //
45 // Class Description:
46 //
47 // Abstract interface to energy loss models
48 
49 // -------------------------------------------------------------------
50 //
51 
52 #include "G4VEmModel.hh"
53 #include "G4ElementData.hh"
54 #include "G4LossTableManager.hh"
55 #include "G4ProductionCutsTable.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "G4Log.hh"
60 #include "Randomize.hh"
61 #include <iostream>
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 
67 
69  flucModel(nullptr),anglModel(nullptr), name(nam), lowLimit(0.1*CLHEP::keV),
70  highLimit(100.0*CLHEP::TeV),eMinActive(0.0),eMaxActive(DBL_MAX),
71  polarAngleLimit(CLHEP::pi),secondaryThreshold(DBL_MAX),
72  theLPMflag(false),flagDeexcitation(false),flagForceBuildTable(false),
73  isMaster(true),fElementData(nullptr),pParticleChange(nullptr),xSectionTable(nullptr),
74  theDensityFactor(nullptr),theDensityIdx(nullptr),lossFlucFlag(true),fCurrentCouple(nullptr),
75  fCurrentElement(nullptr),fCurrentIsotope(nullptr),nsec(5)
76 {
77  xsec.resize(nsec);
78  nSelectors = 0;
79  elmSelectors = nullptr;
80  localElmSelectors = true;
81  localTable = true;
82  useAngularGenerator = false;
83  isLocked = false;
84  idxTable = 0;
85 
86  fEmManager = G4LossTableManager::Instance();
87  fEmManager->Register(this);
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
93 {
94  if(localElmSelectors) {
95  if(nSelectors > 0) {
96  for(G4int i=0; i<nSelectors; ++i) {
97  delete (*elmSelectors)[i];
98  }
99  }
100  delete elmSelectors;
101  }
102  delete anglModel;
103 
104  if(localTable && xSectionTable) {
106  delete xSectionTable;
107  xSectionTable = nullptr;
108  }
109  if(isMaster && fElementData) {
110  delete fElementData;
111  fElementData = nullptr;
112  }
113  fEmManager->DeRegister(this);
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
119 {
120  G4ParticleChangeForLoss* p = nullptr;
121  if (pParticleChange) {
122  p = static_cast<G4ParticleChangeForLoss*>(pParticleChange);
123  } else {
124  p = new G4ParticleChangeForLoss();
125  pParticleChange = p;
126  }
127  return p;
128 }
129 
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 
133 {
134  G4ParticleChangeForGamma* p = nullptr;
135  if (pParticleChange) {
136  p = static_cast<G4ParticleChangeForGamma*>(pParticleChange);
137  } else {
138  p = new G4ParticleChangeForGamma();
139  pParticleChange = p;
140  }
141  return p;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145 
147  const G4DataVector& cuts)
148 {
149  // using spline for element selectors should be investigated in details
150  // because small number of points may provide biased results
151  // large number of points requires significant increase of memory
152  G4bool spline = false;
153 
154  //G4cout << "IES: for " << GetName() << " Emin(MeV)= " << lowLimit/MeV
155  // << " Emax(MeV)= " << highLimit/MeV << G4endl;
156 
157  // two times less bins because probability functon is normalized
158  // so correspondingly is more smooth
159  if(highLimit <= lowLimit) { return; }
160 
162 
163  G4ProductionCutsTable* theCoupleTable=
165  G4int numOfCouples = theCoupleTable->GetTableSize();
166 
167  // prepare vector
168  if(!elmSelectors) {
169  elmSelectors = new std::vector<G4EmElementSelector*>;
170  }
171  if(numOfCouples > nSelectors) {
172  for(G4int i=nSelectors; i<numOfCouples; ++i) {
173  elmSelectors->push_back(nullptr);
174  }
175  nSelectors = numOfCouples;
176  }
177 
178  // initialise vector
179  for(G4int i=0; i<numOfCouples; ++i) {
180 
181  // no need in element selectors for infionite cuts
182  if(cuts[i] == DBL_MAX) { continue; }
183 
184  fCurrentCouple = theCoupleTable->GetMaterialCutsCouple(i);
185  const G4Material* material = fCurrentCouple->GetMaterial();
186 
187  // selector already exist check if should be deleted
188  G4bool create = true;
189  if((*elmSelectors)[i]) {
190  if(material == ((*elmSelectors)[i])->GetMaterial()) { create = false; }
191  else { delete (*elmSelectors)[i]; }
192  }
193  if(create) {
194  G4double emin = std::max(lowLimit,
195  MinPrimaryEnergy(material, part, cuts[i]));
196  G4double emax = std::max(highLimit, 10*emin);
197  static const G4double invlog106 = 1.0/(6*G4Log(10.));
198  G4int nbins = (G4int)(nbinsPerDec*G4Log(emax/emin)*invlog106);
199  nbins = std::max(nbins, 3);
200 
201  (*elmSelectors)[i] = new G4EmElementSelector(this,material,nbins,
202  emin,emax,spline);
203  }
204  ((*elmSelectors)[i])->Initialise(part, cuts[i]);
205  /*
206  G4cout << "G4VEmModel::InitialiseElmSelectors i= " << i
207  << " idx= " << fCurrentCouple->GetIndex()
208  << " " << part->GetParticleName()
209  << " for " << GetName() << " cut= " << cuts[i]
210  << " " << (*elmSelectors)[i] << G4endl;
211  ((*elmSelectors)[i])->Dump(part);
212  */
213  }
214 }
215 
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217 
219  G4VEmModel*)
220 {}
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
223 
225  const G4Material* material)
226 {
227  if(material) {
228  const G4ElementVector* theElementVector = material->GetElementVector();
229  G4int n = material->GetNumberOfElements();
230  for(G4int i=0; i<n; ++i) {
231  G4int Z = ((*theElementVector)[i])->GetZasInt();
232  InitialiseForElement(part, Z);
233  }
234  } else {
235  //G4cout << "G4VEmModel::InitialiseForMaterial for " << GetName();
236  //if(part) { G4cout << " and " << part->GetParticleName(); }
237  //G4cout << " with no material" << G4endl;
238  }
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 
244 {}
245 
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247 
249  const G4ParticleDefinition*,
251 {
252  return 0.0;
253 }
254 
255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
256 
258  const G4ParticleDefinition* p,
259  G4double ekin,
260  G4double emin,
261  G4double emax)
262 {
263  SetupForMaterial(p, material, ekin);
264  G4double cross = 0.0;
265  const G4ElementVector* theElementVector = material->GetElementVector();
266  const G4double* theAtomNumDensityVector =
267  material->GetVecNbOfAtomsPerVolume();
268  G4int nelm = material->GetNumberOfElements();
269  if(nelm > nsec) {
270  xsec.resize(nelm);
271  nsec = nelm;
272  }
273  for (G4int i=0; i<nelm; ++i) {
274  cross += theAtomNumDensityVector[i]*
275  ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
276  xsec[i] = cross;
277  }
278  return cross;
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
282 
284  const G4ParticleDefinition*,
285  G4double)
286 {
287  return 0.0;
288 }
289 
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291 
293 {}
294 
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
296 
298  const G4ParticleDefinition* pd,
299  G4double kinEnergy,
300  G4double tcut,
301  G4double tmax)
302 {
303  const G4ElementVector* theElementVector = material->GetElementVector();
304  G4int n = material->GetNumberOfElements() - 1;
305  fCurrentElement = (*theElementVector)[n];
306  if (n > 0) {
308  G4VEmModel::CrossSectionPerVolume(material,pd,kinEnergy,tcut,tmax);
309  for(G4int i=0; i<n; ++i) {
310  if (x <= xsec[i]) {
311  fCurrentElement = (*theElementVector)[i];
312  break;
313  }
314  }
315  }
316  return fCurrentElement;
317 }
318 
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320 
324 {
325  return 0.0;
326 }
327 
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
329 
330 G4double
332  G4int, G4int,
334 {
335  return 0.0;
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
339 
341 {}
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
344 
346 {
348  track.GetMaterial(), track.GetKineticEnergy());
349 }
350 
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352 
354  const G4Material*, G4double)
355 {
356  G4double q = p->GetPDGCharge()*inveplus;
357  return q*q;
358 }
359 
360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
361 
363  const G4Material*, G4double)
364 {
365  return p->GetPDGCharge();
366 }
367 
368 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
369 
371  const G4DynamicParticle*,
373 {}
374 
375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
376 
378  const G4ParticleDefinition* p, G4double e)
379 {
380  SetCurrentCouple(couple);
381  return e*e*CrossSectionPerVolume(couple->GetMaterial(),p,e,0.0,DBL_MAX);
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385 
387  const G4ParticleDefinition*,
388  G4double)
389 {
390  return 0.0;
391 }
392 
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
394 
396  const G4MaterialCutsCouple*)
397 {
398  return 0.0;
399 }
400 
401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
402 
404  G4double kineticEnergy)
405 {
406  return kineticEnergy;
407 }
408 
409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
410 
412  const G4Material*, G4double)
413 {}
414 
415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
416 
417 void
419 {
420  if(p && pParticleChange != p) { pParticleChange = p; }
421  flucModel = f;
422 }
423 
424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
425 
427 {
428  if(p != xSectionTable) {
429  if(xSectionTable && localTable) {
431  delete xSectionTable;
432  }
433  xSectionTable = p;
434  }
435  localTable = isLocal;
436 }
437 
438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
439 
440 void G4VEmModel::ModelDescription(std::ostream& outFile) const
441 {
442  outFile << "The description for this model has not been written yet.\n";
443 }
444 
445 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4int NumberOfBinsPerDecade() const
size_t idxTable
Definition: G4VEmModel.hh:428
const XML_Char * name
Definition: expat.h:151
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:257
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
Definition: G4VEmModel.cc:243
virtual void StartTracking(G4Track *)
Definition: G4VEmModel.cc:292
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
std::vector< G4Element * > G4ElementVector
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
Definition: G4VEmModel.cc:370
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
Definition: G4VEmModel.cc:418
void DeRegister(G4VEnergyLossProcess *p)
virtual void DefineForRegion(const G4Region *)
Definition: G4VEmModel.cc:340
const char * p
Definition: xmltok.h:285
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
Definition: G4VEmModel.cc:386
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
Definition: G4VEmModel.cc:218
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:411
tuple x
Definition: test.py:50
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
Definition: G4VEmModel.cc:395
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)=0
static constexpr double TeV
Definition: G4SIunits.hh:218
string material
Definition: eplot.py:19
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
G4double GetKineticEnergy() const
virtual void InitialiseForMaterial(const G4ParticleDefinition *, const G4Material *)
Definition: G4VEmModel.cc:224
#define G4UniformRand()
Definition: Randomize.hh:97
virtual void ModelDescription(std::ostream &outFile) const
Definition: G4VEmModel.cc:440
bool G4bool
Definition: G4Types.hh:79
static const G4double inveplus
Definition: G4VEmModel.hh:430
const G4ParticleDefinition * GetParticleDefinition() const
virtual ~G4VEmModel()
Definition: G4VEmModel.cc:92
static constexpr double eplus
void SetCrossSectionTable(G4PhysicsTable *, G4bool isLocal)
Definition: G4VEmModel.cc:426
void Register(G4VEnergyLossProcess *p)
virtual G4double GetPartialCrossSection(const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:283
const G4int n
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:321
virtual G4double ChargeSquareRatio(const G4Track &)
Definition: G4VEmModel.cc:345
G4Material * GetMaterial() const
static const G4double emax
G4double G4Log(G4double x)
Definition: G4Log.hh:230
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:353
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:362
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
Definition: G4VEmModel.cc:248
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:447
static G4EmParameters * Instance()
virtual G4double ComputeCrossSectionPerShell(const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:331
G4VParticleChange * pParticleChange
Definition: G4VEmModel.hh:424
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:403
static constexpr double pi
Definition: G4SIunits.hh:75
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4PhysicsTable * xSectionTable
Definition: G4VEmModel.hh:425
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
static constexpr double keV
Definition: G4SIunits.hh:216
G4ElementData * fElementData
Definition: G4VEmModel.hh:423
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
void clearAndDestroy()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:377
const G4Material * GetMaterial() const