Geant4  10.02.p02
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 93264 2015-10-14 09:30:04Z 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 const G4double log106 = 6*G4Log(10.);
68 
70  flucModel(nullptr),anglModel(nullptr), name(nam), lowLimit(0.1*CLHEP::keV),
71  highLimit(100.0*CLHEP::TeV),eMinActive(0.0),eMaxActive(DBL_MAX),
72  polarAngleLimit(CLHEP::pi),secondaryThreshold(DBL_MAX),
73  theLPMflag(false),flagDeexcitation(false),flagForceBuildTable(false),
74  isMaster(true),fElementData(nullptr),pParticleChange(nullptr),xSectionTable(nullptr),
75  theDensityFactor(nullptr),theDensityIdx(nullptr),fCurrentCouple(nullptr),
76  fCurrentElement(nullptr),fCurrentIsotope(nullptr),nsec(5)
77 {
78  xsec.resize(nsec);
79  nSelectors = 0;
80  elmSelectors = nullptr;
81  localElmSelectors = true;
82  localTable = true;
83  useAngularGenerator = false;
84  isLocked = false;
85  idxTable = 0;
86 
88  fManager->Register(this);
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 
94 {
95  if(localElmSelectors) {
96  if(nSelectors > 0) {
97  for(G4int i=0; i<nSelectors; ++i) {
98  delete (*elmSelectors)[i];
99  }
100  }
101  delete elmSelectors;
102  }
103  delete anglModel;
104 
105  if(localTable && xSectionTable) {
107  delete xSectionTable;
108  xSectionTable = nullptr;
109  }
110  if(isMaster && fElementData) {
111  delete fElementData;
112  fElementData = nullptr;
113  }
114 
115  fManager->DeRegister(this);
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119 
121 {
122  G4ParticleChangeForLoss* p = nullptr;
123  if (pParticleChange) {
124  p = static_cast<G4ParticleChangeForLoss*>(pParticleChange);
125  } else {
126  p = new G4ParticleChangeForLoss();
127  pParticleChange = p;
128  }
129  return p;
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133 
135 {
136  G4ParticleChangeForGamma* p = nullptr;
137  if (pParticleChange) {
138  p = static_cast<G4ParticleChangeForGamma*>(pParticleChange);
139  } else {
140  p = new G4ParticleChangeForGamma();
141  pParticleChange = p;
142  }
143  return p;
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 
149  const G4DataVector& cuts)
150 {
151  // using spline for element selectors should be investigated in details
152  // because small number of points may provide biased results
153  // large number of points requires significant increase of memory
154  //G4bool spline = fManager->SplineFlag();
155  G4bool spline = false;
156 
157  //G4cout << "IES: for " << GetName() << " Emin(MeV)= " << lowLimit/MeV
158  // << " Emax(MeV)= " << highLimit/MeV << G4endl;
159 
160  // two times less bins because probability functon is normalized
161  // so correspondingly is more smooth
162  if(highLimit <= lowLimit) { return; }
163 
164  G4ProductionCutsTable* theCoupleTable=
166  G4int numOfCouples = theCoupleTable->GetTableSize();
167 
168  // prepare vector
169  if(!elmSelectors) {
170  elmSelectors = new std::vector<G4EmElementSelector*>;
171  }
172  if(numOfCouples > nSelectors) {
173  for(G4int i=nSelectors; i<numOfCouples; ++i) {
174  elmSelectors->push_back(nullptr);
175  }
176  nSelectors = numOfCouples;
177  }
178 
179  // initialise vector
180  for(G4int i=0; i<numOfCouples; ++i) {
181 
182  // no need in element selectors for infionite cuts
183  if(cuts[i] == DBL_MAX) { continue; }
184 
185  fCurrentCouple = theCoupleTable->GetMaterialCutsCouple(i);
186  const G4Material* material = fCurrentCouple->GetMaterial();
187 
188  // selector already exist check if should be deleted
189  G4bool create = true;
190  if((*elmSelectors)[i]) {
191  if(material == ((*elmSelectors)[i])->GetMaterial()) { create = false; }
192  else { delete (*elmSelectors)[i]; }
193  }
194  if(create) {
195  G4double emin = std::max(lowLimit,
196  MinPrimaryEnergy(material, part, cuts[i]));
197  G4double emax = std::max(highLimit, 10*emin);
199  *G4Log(emax/emin)/log106);
200  nbins = std::max(nbins, 3);
201 
202  (*elmSelectors)[i] = new G4EmElementSelector(this,material,nbins,
203  emin,emax,spline);
204  }
205  ((*elmSelectors)[i])->Initialise(part, cuts[i]);
206  /*
207  G4cout << "G4VEmModel::InitialiseElmSelectors i= " << i
208  << " idx= " << fCurrentCouple->GetIndex()
209  << " " << part->GetParticleName()
210  << " for " << GetName() << " cut= " << cuts[i]
211  << " " << (*elmSelectors)[i] << G4endl;
212  ((*elmSelectors)[i])->Dump(part);
213  */
214  }
215 }
216 
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218 
220  G4VEmModel*)
221 {}
222 
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
224 
226  const G4Material* material)
227 {
228  if(material) {
229  const G4ElementVector* theElementVector = material->GetElementVector();
230  G4int n = material->GetNumberOfElements();
231  for(G4int i=0; i<n; ++i) {
232  G4int Z = G4lrint(((*theElementVector)[i])->GetZ());
233  InitialiseForElement(part, Z);
234  }
235  } else {
236  //G4cout << "G4VEmModel::InitialiseForMaterial for " << GetName();
237  //if(part) { G4cout << " and " << part->GetParticleName(); }
238  //G4cout << " with no material" << G4endl;
239  }
240 }
241 
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
243 
245 {}
246 
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
248 
250  const G4ParticleDefinition*,
252 {
253  return 0.0;
254 }
255 
256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
257 
259  const G4ParticleDefinition* p,
260  G4double ekin,
261  G4double emin,
262  G4double emax)
263 {
264  SetupForMaterial(p, material, ekin);
265  G4double cross = 0.0;
266  const G4ElementVector* theElementVector = material->GetElementVector();
267  const G4double* theAtomNumDensityVector =
268  material->GetVecNbOfAtomsPerVolume();
269  G4int nelm = material->GetNumberOfElements();
270  if(nelm > nsec) {
271  xsec.resize(nelm);
272  nsec = nelm;
273  }
274  for (G4int i=0; i<nelm; ++i) {
275  cross += theAtomNumDensityVector[i]*
276  ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
277  xsec[i] = cross;
278  }
279  return cross;
280 }
281 
282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
283 
285 {}
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 
290  const G4ParticleDefinition* pd,
291  G4double kinEnergy,
292  G4double tcut,
293  G4double tmax)
294 {
295  const G4ElementVector* theElementVector = material->GetElementVector();
296  G4int n = material->GetNumberOfElements() - 1;
297  fCurrentElement = (*theElementVector)[n];
298  if (n > 0) {
300  G4VEmModel::CrossSectionPerVolume(material,pd,kinEnergy,tcut,tmax);
301  for(G4int i=0; i<n; ++i) {
302  if (x <= xsec[i]) {
303  fCurrentElement = (*theElementVector)[i];
304  break;
305  }
306  }
307  }
308  return fCurrentElement;
309 }
310 
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312 
316 {
317  return 0.0;
318 }
319 
320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
321 
322 G4double
324  G4int, G4int,
326 {
327  return 0.0;
328 }
329 
330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
331 
333 {}
334 
335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
336 
338 {
340  track.GetMaterial(), track.GetKineticEnergy());
341 }
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
344 
346  const G4Material*, G4double)
347 {
348  G4double q = p->GetPDGCharge()*inveplus;
349  return q*q;
350 }
351 
352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
353 
355  const G4Material*, G4double)
356 {
357  return p->GetPDGCharge();
358 }
359 
360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
361 
363  const G4DynamicParticle*,
365 {}
366 
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
368 
370  const G4ParticleDefinition* p, G4double e)
371 {
372  SetCurrentCouple(couple);
373  return e*e*CrossSectionPerVolume(couple->GetMaterial(),p,e,0.0,DBL_MAX);
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 
379  const G4ParticleDefinition*,
380  G4double)
381 {
382  return 0.0;
383 }
384 
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
386 
388  const G4MaterialCutsCouple*)
389 {
390  return 0.0;
391 }
392 
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
394 
396  G4double kineticEnergy)
397 {
398  return kineticEnergy;
399 }
400 
401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
402 
404  const G4Material*, G4double)
405 {}
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408 
409 void
411 {
412  if(p && pParticleChange != p) { pParticleChange = p; }
413  flucModel = f;
414 }
415 
416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
417 
419 {
420  if(p != xSectionTable) {
421  if(xSectionTable && localTable) {
423  delete xSectionTable;
424  }
425  xSectionTable = p;
426  }
427  localTable = isLocal;
428 }
429 
430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
431 
432 void G4VEmModel::ModelDescription(std::ostream& outFile) const
433 {
434  outFile << "The description for this model has not been written yet.\n";
435 }
436 
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool localTable
Definition: G4VEmModel.hh:412
size_t idxTable
Definition: G4VEmModel.hh:426
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:258
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
Definition: G4VEmModel.cc:244
virtual void StartTracking(G4Track *)
Definition: G4VEmModel.cc:284
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:120
std::vector< G4Element * > G4ElementVector
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
Definition: G4VEmModel.cc:362
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:148
void DeRegister(G4VEnergyLossProcess *p)
virtual void DefineForRegion(const G4Region *)
Definition: G4VEmModel.cc:332
G4String name
Definition: TRTMaterials.hh:40
G4bool isMaster
Definition: G4VEmModel.hh:410
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
Definition: G4VEmModel.cc:378
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
Definition: G4VEmModel.cc:219
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:403
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
Definition: G4VEmModel.cc:387
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
const G4Element * fCurrentElement
Definition: G4VEmModel.hh:435
G4VEmFluctuationModel * flucModel
Definition: G4VEmModel.hh:395
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:69
G4int nSelectors
Definition: G4VEmModel.hh:416
G4bool isLocked
Definition: G4VEmModel.hh:415
G4LossTableManager * fManager
Definition: G4VEmModel.hh:433
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)=0
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
G4double GetKineticEnergy() const
virtual void InitialiseForMaterial(const G4ParticleDefinition *, const G4Material *)
Definition: G4VEmModel.cc:225
#define G4UniformRand()
Definition: Randomize.hh:97
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=0)
Definition: G4VEmModel.cc:410
virtual void ModelDescription(std::ostream &outFile) const
Definition: G4VEmModel.cc:432
const G4MaterialCutsCouple * fCurrentCouple
Definition: G4VEmModel.hh:434
std::vector< G4EmElementSelector * > * elmSelectors
Definition: G4VEmModel.hh:417
bool G4bool
Definition: G4Types.hh:79
static const G4double inveplus
Definition: G4VEmModel.hh:427
const G4ParticleDefinition * GetParticleDefinition() const
virtual ~G4VEmModel()
Definition: G4VEmModel.cc:93
G4VEmAngularDistribution * anglModel
Definition: G4VEmModel.hh:396
void SetCrossSectionTable(G4PhysicsTable *, G4bool isLocal)
Definition: G4VEmModel.cc:418
void Register(G4VEnergyLossProcess *p)
const G4int n
G4int nsec
Definition: G4VEmModel.hh:438
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:313
virtual G4double ChargeSquareRatio(const G4Track &)
Definition: G4VEmModel.cc:337
G4Material * GetMaterial() const
static const G4double emax
G4double G4Log(G4double x)
Definition: G4Log.hh:230
std::vector< G4double > xsec
Definition: G4VEmModel.hh:439
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:345
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:354
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
Definition: G4VEmModel.cc:249
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static const double pi
Definition: G4SIunits.hh:74
int G4lrint(double ad)
Definition: templates.hh:163
G4bool localElmSelectors
Definition: G4VEmModel.hh:413
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4int GetNumberOfBinsPerDecade() const
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:445
const G4double x[NPOINTSGL]
G4double highLimit
Definition: G4VEmModel.hh:402
G4bool useAngularGenerator
Definition: G4VEmModel.hh:414
G4double lowLimit
Definition: G4VEmModel.hh:401
virtual G4double ComputeCrossSectionPerShell(const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:323
G4VParticleChange * pParticleChange
Definition: G4VEmModel.hh:422
const G4double log106
Definition: G4VEmModel.cc:67
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:395
static const double TeV
Definition: G4SIunits.hh:215
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const double keV
Definition: G4SIunits.hh:213
G4PhysicsTable * xSectionTable
Definition: G4VEmModel.hh:423
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:196
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
G4ElementData * fElementData
Definition: G4VEmModel.hh:421
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:134
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:369
const G4Material * GetMaterial() const