ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VEmModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VEmModel.cc
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 // GEANT4 Class file
29 //
30 //
31 // File name: G4VEmModel
32 //
33 // Author: Vladimir Ivanchenko
34 //
35 // Creation date: 25.07.2005
36 //
37 // Modifications:
38 // 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko)
39 // 06.02.2006 add method ComputeMeanFreePath() (mma)
40 // 16.02.2009 Move implementations of virtual methods to source (VI)
41 //
42 //
43 // Class Description:
44 //
45 // Abstract interface to energy loss models
46 
47 // -------------------------------------------------------------------
48 //
49 
50 #include "G4VEmModel.hh"
51 #include "G4ElementData.hh"
52 #include "G4LossTableManager.hh"
53 #include "G4LossTableBuilder.hh"
54 #include "G4ProductionCutsTable.hh"
55 #include "G4ParticleChangeForLoss.hh"
56 #include "G4ParticleChangeForGamma.hh"
57 #include "G4SystemOfUnits.hh"
58 #include "G4Log.hh"
59 #include "Randomize.hh"
60 #include <iostream>
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 
65 G4VEmModel::G4VEmModel(const G4String& nam):
66  flucModel(nullptr),anglModel(nullptr), name(nam), lowLimit(0.1*CLHEP::keV),
67  highLimit(100.0*CLHEP::TeV),eMinActive(0.0),eMaxActive(DBL_MAX),
68  polarAngleLimit(CLHEP::pi),secondaryThreshold(DBL_MAX),
69  theLPMflag(false),flagDeexcitation(false),flagForceBuildTable(false),
70  isMaster(true),fElementData(nullptr),pParticleChange(nullptr),
71  xSectionTable(nullptr),pBaseMaterial(nullptr),idxTable(0),
72  lossFlucFlag(true),inveplus(1.0/CLHEP::eplus),pFactor(1.0),
73  fCurrentCouple(nullptr),fCurrentElement(nullptr),fCurrentIsotope(nullptr),
74  fTripletModel(nullptr),nsec(5)
75 {
76  xsec.resize(nsec);
77  nSelectors = 0;
78  elmSelectors = nullptr;
79  localElmSelectors = true;
80  localTable = true;
81  useAngularGenerator = false;
82  useBaseMaterials = true;
83  isLocked = false;
84 
85  fEmManager = G4LossTableManager::Instance();
86  fEmManager->Register(this);
87  G4LossTableBuilder* bld = fEmManager->GetTableBuilder();
88  theDensityFactor = bld->GetDensityFactors();
89  theDensityIdx = bld->GetCoupleIndexes();
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93 
94 G4VEmModel::~G4VEmModel()
95 {
96  if(localElmSelectors) {
97  for(G4int i=0; i<nSelectors; ++i) {
98  delete (*elmSelectors)[i];
99  }
100  delete elmSelectors;
101  }
102  delete anglModel;
103 
104  if(localTable && xSectionTable) {
105  xSectionTable->clearAndDestroy();
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 
118 G4ParticleChangeForLoss* G4VEmModel::GetParticleChangeForLoss()
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  if(fTripletModel) { fTripletModel->SetParticleChange(p); }
128  return p;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 G4ParticleChangeForGamma* G4VEmModel::GetParticleChangeForGamma()
134 {
135  G4ParticleChangeForGamma* p = nullptr;
136  if (pParticleChange) {
137  p = static_cast<G4ParticleChangeForGamma*>(pParticleChange);
138  } else {
139  p = new G4ParticleChangeForGamma();
140  pParticleChange = p;
141  }
142  if(fTripletModel) { fTripletModel->SetParticleChange(p); }
143  return p;
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 
148 void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* part,
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 = false;
155 
156  //G4cout << "IES: for " << GetName() << " Emin(MeV)= " << lowLimit/MeV
157  // << " Emax(MeV)= " << highLimit/MeV << G4endl;
158 
159  // two times less bins because probability functon is normalized
160  // so correspondingly is more smooth
161  if(highLimit <= lowLimit) { return; }
162 
163  G4int nbinsPerDec = G4EmParameters::Instance()->NumberOfBinsPerDecade();
164 
165  G4ProductionCutsTable* theCoupleTable=
166  G4ProductionCutsTable::GetProductionCutsTable();
167  G4int numOfCouples = theCoupleTable->GetTableSize();
168 
169  // prepare vector
170  if(!elmSelectors) {
171  elmSelectors = new std::vector<G4EmElementSelector*>;
172  }
173  if(numOfCouples > nSelectors) {
174  for(G4int i=nSelectors; i<numOfCouples; ++i) {
175  elmSelectors->push_back(nullptr);
176  }
177  nSelectors = numOfCouples;
178  }
179 
180  // initialise vector
181  for(G4int i=0; i<numOfCouples; ++i) {
182 
183  // no need in element selectors for infionite cuts
184  if(cuts[i] == DBL_MAX) { continue; }
185 
186  auto couple = theCoupleTable->GetMaterialCutsCouple(i);
187  auto material = couple->GetMaterial();
188  SetCurrentCouple(couple);
189 
190  // selector already exist check if should be deleted
191  G4bool create = true;
192  if((*elmSelectors)[i]) {
193  if(material == ((*elmSelectors)[i])->GetMaterial()) { create = false; }
194  else { delete (*elmSelectors)[i]; }
195  }
196  if(create) {
197  G4double emin = std::max(lowLimit,
198  MinPrimaryEnergy(material, part, cuts[i]));
199  G4double emax = std::max(highLimit, 10*emin);
200  static const G4double invlog106 = 1.0/(6*G4Log(10.));
201  G4int nbins = (G4int)(nbinsPerDec*G4Log(emax/emin)*invlog106);
202  nbins = std::max(nbins, 3);
203 
204  (*elmSelectors)[i] = new G4EmElementSelector(this,material,nbins,
205  emin,emax,spline);
206  }
207  ((*elmSelectors)[i])->Initialise(part, cuts[i]);
208  /*
209  G4cout << "G4VEmModel::InitialiseElmSelectors i= " << i
210  << " " << part->GetParticleName()
211  << " for " << GetName() << " cut= " << cuts[i]
212  << " " << (*elmSelectors)[i] << G4endl;
213  ((*elmSelectors)[i])->Dump(part);
214  */
215  }
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219 
220 void G4VEmModel::InitialiseLocal(const G4ParticleDefinition*,
221  G4VEmModel*)
222 {}
223 
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
225 
226 void G4VEmModel::InitialiseForMaterial(const G4ParticleDefinition* part,
227  const G4Material* material)
228 {
229  if(material) {
230  size_t n = material->GetNumberOfElements();
231  for(size_t i=0; i<n; ++i) {
232  G4int Z = material->GetElement(i)->GetZasInt();
233  InitialiseForElement(part, Z);
234  }
235  }
236 }
237 
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
239 
240 void G4VEmModel::InitialiseForElement(const G4ParticleDefinition*, G4int)
241 {}
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244 
245 G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
246  const G4ParticleDefinition*,
247  G4double,G4double)
248 {
249  return 0.0;
250 }
251 
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
253 
254 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
255  const G4ParticleDefinition* p,
256  G4double ekin,
257  G4double emin,
258  G4double emax)
259 {
260  SetupForMaterial(p, material, ekin);
261  G4double cross = 0.0;
262  const G4double* theAtomNumDensityVector =
263  material->GetVecNbOfAtomsPerVolume();
264  G4int nelm = material->GetNumberOfElements();
265  if(nelm > nsec) {
266  xsec.resize(nelm);
267  nsec = nelm;
268  }
269  for (G4int i=0; i<nelm; ++i) {
270  cross += theAtomNumDensityVector[i]*
271  ComputeCrossSectionPerAtom(p,material->GetElement(i),ekin,emin,emax);
272  xsec[i] = cross;
273  }
274  return cross;
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278 
279 G4double G4VEmModel::GetPartialCrossSection(const G4Material*, G4int,
280  const G4ParticleDefinition*,
281  G4double)
282 {
283  return 0.0;
284 }
285 
286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
287 
288 void G4VEmModel::StartTracking(G4Track*)
289 {}
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
292 
293 const G4Element* G4VEmModel::SelectRandomAtom(const G4Material* material,
294  const G4ParticleDefinition* pd,
295  G4double kinEnergy,
296  G4double tcut,
297  G4double tmax)
298 {
299  size_t n = material->GetNumberOfElements();
300  fCurrentElement = material->GetElement(0);
301  if (n > 1) {
302  G4double x = G4UniformRand()*
303  G4VEmModel::CrossSectionPerVolume(material,pd,kinEnergy,tcut,tmax);
304  for(size_t i=0; i<n; ++i) {
305  if (x <= xsec[i]) {
306  fCurrentElement = material->GetElement(i);
307  break;
308  }
309  }
310  }
311  return fCurrentElement;
312 }
313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
314 
315 G4int G4VEmModel::SelectRandomAtomNumber(const G4Material* mat)
316 {
317  // this algorith assumes that cross section is proportional to
318  // number electrons multiplied by number of atoms
319  size_t nn = mat->GetNumberOfElements();
320  fCurrentElement = mat->GetElement(0);
321  if(1 < nn) {
322  const G4double* at = mat->GetVecNbOfAtomsPerVolume();
323  G4double tot = mat->GetTotNbOfAtomsPerVolume()*G4UniformRand();
324  for(size_t i=0; i<nn; ++i) {
325  tot -= at[i];
326  if(tot <= 0.0) {
327  fCurrentElement = mat->GetElement(i);
328  break;
329  }
330  }
331  }
332  return fCurrentElement->GetZasInt();
333 }
334 
335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
336 
337 G4int G4VEmModel::SelectIsotopeNumber(const G4Element* elm)
338 {
339  SetCurrentElement(elm);
340  size_t ni = elm->GetNumberOfIsotopes();
341  fCurrentIsotope = elm->GetIsotope(0);
342  size_t idx = 0;
343  if(ni > 1) {
344  const G4double* ab = elm->GetRelativeAbundanceVector();
345  G4double x = G4UniformRand();
346  for(; idx<ni; ++idx) {
347  x -= ab[idx];
348  if (x <= 0.0) {
349  fCurrentIsotope = elm->GetIsotope(idx);
350  break;
351  }
352  }
353  }
354  return fCurrentIsotope->GetN();
355 }
356 
357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
358 
359 G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
360  G4double, G4double, G4double,
361  G4double, G4double)
362 {
363  return 0.0;
364 }
365 
366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
367 
368 G4double
369 G4VEmModel::ComputeCrossSectionPerShell(const G4ParticleDefinition*,
370  G4int, G4int,
371  G4double, G4double, G4double)
372 {
373  return 0.0;
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 
378 void G4VEmModel::DefineForRegion(const G4Region*)
379 {}
380 
381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
382 
383 G4double G4VEmModel::ChargeSquareRatio(const G4Track& track)
384 {
385  return GetChargeSquareRatio(track.GetParticleDefinition(),
386  track.GetMaterial(), track.GetKineticEnergy());
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390 
391 G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
392  const G4Material*, G4double)
393 {
394  G4double q = p->GetPDGCharge()*inveplus;
395  return q*q;
396 }
397 
398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
399 
400 G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
401  const G4Material*, G4double)
402 {
403  return p->GetPDGCharge();
404 }
405 
406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
407 
408 void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
409  const G4DynamicParticle*,
410  G4double&,G4double&,G4double)
411 {}
412 
413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
414 
415 G4double G4VEmModel::Value(const G4MaterialCutsCouple* couple,
416  const G4ParticleDefinition* p, G4double e)
417 {
418  SetCurrentCouple(couple);
419  return pFactor*e*e*CrossSectionPerVolume(pBaseMaterial,p,e,0.0,DBL_MAX);
420 }
421 
422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
423 
424 G4double G4VEmModel::MinPrimaryEnergy(const G4Material*,
425  const G4ParticleDefinition*,
426  G4double)
427 {
428  return 0.0;
429 }
430 
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
432 
433 G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,
434  const G4MaterialCutsCouple*)
435 {
436  return 0.0;
437 }
438 
439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
440 
441 G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
442  G4double kineticEnergy)
443 {
444  return kineticEnergy;
445 }
446 
447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
448 
449 void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
450  const G4Material*, G4double)
451 {}
452 
453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
454 
455 void
456 G4VEmModel::SetParticleChange(G4VParticleChange* p, G4VEmFluctuationModel* f)
457 {
458  if(p && pParticleChange != p) { pParticleChange = p; }
459  if(flucModel != f) { flucModel = f; }
460 }
461 
462 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
463 
464 void G4VEmModel::SetCrossSectionTable(G4PhysicsTable* p, G4bool isLocal)
465 {
466  if(p != xSectionTable) {
467  if(xSectionTable && localTable) {
468  xSectionTable->clearAndDestroy();
469  delete xSectionTable;
470  }
471  xSectionTable = p;
472  }
473  localTable = isLocal;
474 }
475 
476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
477 
478 void G4VEmModel::ModelDescription(std::ostream& outFile) const
479 {
480  outFile << "The description for this model has not been written yet.\n";
481 }
482 
483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......