Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EmModelManager.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: G4EmModelManager.cc 97742 2016-06-08 09:24:54Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4EmModelManager
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 07.05.2002
38 //
39 // Modifications:
40 //
41 // 23-12-02 V.Ivanchenko change interface in order to move
42 // to cut per region
43 // 20-01-03 Migrade to cut per region (V.Ivanchenko)
44 // 24-01-03 Make models region aware (V.Ivanchenko)
45 // 13-02-03 The set of models is defined for region (V.Ivanchenko)
46 // 06-03-03 Fix in energy intervals for models (V.Ivanchenko)
47 // 13-04-03 Add startFromNull (V.Ivanchenko)
48 // 13-05-03 Add calculation of precise range (V.Ivanchenko)
49 // 16-07-03 Replace G4Material by G4MaterialCutCouple in dE/dx and CrossSection
50 // calculation (V.Ivanchenko)
51 // 21-07-03 Add UpdateEmModel method (V.Ivanchenko)
52 // 03-11-03 Substitute STL vector for G4RegionModels (V.Ivanchenko)
53 // 26-01-04 Fix in energy range conditions (V.Ivanchenko)
54 // 24-03-05 Remove check or IsInCharge (V.Ivanchenko)
55 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
56 // 18-08-05 Fix cut for e+e- pair production (V.Ivanchenko)
57 // 29-11-05 Add protection for arithmetic operations with cut=DBL_MAX (V.Ivanchenko)
58 // 20-01-06 Introduce G4EmTableType and reducing number of methods (VI)
59 // 13-05-06 Add GetModel by index method (VI)
60 // 15-03-07 Add maxCutInRange (V.Ivanchenko)
61 // 12-04-07 Add verbosity at destruction (V.Ivanchenko)
62 // 08-04-08 Fixed and simplified initialisation of G4RegionModel (VI)
63 // 03-08-09 Create internal vectors only it is needed (VI)
64 // 14-07-11 Use pointer to the vector of cuts and not local copy (VI)
65 //
66 // Class Description:
67 //
68 // It is the unified energy loss process it calculates the continuous
69 // energy loss for charged particles using a set of Energy Loss
70 // models valid for different energy regions. There are a possibility
71 // to create and access to dE/dx and range tables, or to calculate
72 // that information on fly.
73 // -------------------------------------------------------------------
74 //
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
78 #include "G4EmModelManager.hh"
79 #include "G4SystemOfUnits.hh"
80 #include "G4PhysicsTable.hh"
81 #include "G4PhysicsVector.hh"
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
85 G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& indx,
86  G4DataVector& lowE, const G4Region* reg)
87 {
88  nModelsForRegion = nMod;
89  theListOfModelIndexes = new G4int [nModelsForRegion];
90  lowKineticEnergy = new G4double [nModelsForRegion+1];
91  for (G4int i=0; i<nModelsForRegion; ++i) {
92  theListOfModelIndexes[i] = indx[i];
93  lowKineticEnergy[i] = lowE[i];
94  }
95  lowKineticEnergy[nModelsForRegion] = lowE[nModelsForRegion];
96  theRegion = reg;
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
101 G4RegionModels::~G4RegionModels()
102 {
103  delete [] theListOfModelIndexes;
104  delete [] lowKineticEnergy;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 
109 #include "G4Step.hh"
110 #include "G4ParticleDefinition.hh"
111 #include "G4PhysicsVector.hh"
112 #include "G4MaterialCutsCouple.hh"
113 #include "G4ProductionCutsTable.hh"
114 #include "G4RegionStore.hh"
115 #include "G4Gamma.hh"
116 #include "G4Electron.hh"
117 #include "G4Positron.hh"
118 #include "G4UnitsTable.hh"
119 #include "G4DataVector.hh"
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122 
124  nEmModels(0),
125  nRegions(0),
126  particle(0),
127  verboseLevel(0)
128 {
129  maxSubCutInRange = 0.7*mm;
130  models.reserve(4);
131  flucModels.reserve(4);
132  regions.reserve(4);
133  orderOfModels.reserve(4);
134  isUsed.reserve(4);
135  severalModels = true;
136  fluoFlag = false;
137  currRegionModel = nullptr;
138  currModel = nullptr;
139  theCuts = nullptr;
140  theCutsNew = nullptr;
141  theSubCuts = nullptr;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 
147 {
148  verboseLevel = 0; // no verbosity at destruction
149  Clear();
150  delete theCutsNew;
151  delete theSubCuts;
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 
157 {
158  if(1 < verboseLevel) {
159  G4cout << "G4EmModelManager::Clear()" << G4endl;
160  }
161  size_t n = setOfRegionModels.size();
162  if(n > 0) {
163  for(size_t i=0; i<n; ++i) {
164  delete setOfRegionModels[i];
165  setOfRegionModels[i] = nullptr;
166  }
167  }
168 }
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 
173  G4VEmFluctuationModel* fm, const G4Region* r)
174 {
175  if(!p) {
176  G4cout << "G4EmModelManager::AddEmModel WARNING: no model defined."
177  << G4endl;
178  return;
179  }
180  models.push_back(p);
181  flucModels.push_back(fm);
182  regions.push_back(r);
183  orderOfModels.push_back(num);
184  isUsed.push_back(0);
185  p->DefineForRegion(r);
186  ++nEmModels;
187 }
188 
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
190 
192  G4double emin, G4double emax)
193 {
194  if (nEmModels > 0) {
195  for(G4int i=0; i<nEmModels; ++i) {
196  if(nam == models[i]->GetName()) {
197  models[i]->SetLowEnergyLimit(emin);
198  models[i]->SetHighEnergyLimit(emax);
199  break;
200  }
201  }
202  }
203  G4cout << "G4EmModelManager::UpdateEmModel WARNING: no model <"
204  << nam << "> is found out"
205  << G4endl;
206 }
207 
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209 
211 {
212  G4VEmModel* model = nullptr;
213  if(i < nEmModels) { model = models[i]; }
214  else if(verboseLevel > 0 && ver) {
215  G4cout << "G4EmModelManager::GetModel WARNING: "
216  << "index " << i << " is wrong Nmodels= "
217  << nEmModels;
218  if(particle) { G4cout << " for " << particle->GetParticleName(); }
219  G4cout<< G4endl;
220  }
221  return model;
222 }
223 
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
225 
227 {
228  G4RegionModels* rm = setOfRegionModels[idxOfRegionModels[idx]];
229  G4VEmModel* mod = models[rm->ModelIndex(k)];
230  return mod;
231 }
232 
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
234 
236 {
237  G4RegionModels* rm = setOfRegionModels[idxOfRegionModels[idx]];
238  return rm->NumberOfModels();
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
242 
243 const G4DataVector*
245  const G4ParticleDefinition* secondaryParticle,
246  G4double minSubRange,
247  G4int val)
248 {
249  verboseLevel = val;
250  G4String partname = p->GetParticleName();
251  if(1 < verboseLevel) {
252  G4cout << "G4EmModelManager::Initialise() for "
253  << partname << " Nmodels= " << nEmModels << G4endl;
254  }
255  // Are models defined?
256  if(nEmModels < 1) {
258  ed << "No models found out for " << p->GetParticleName()
259  << " !";
260  G4Exception("G4EmModelManager::Initialise","em0002",
261  FatalException, ed);
262  }
263 
264  particle = p;
265  Clear(); // needed if run is not first
266 
267  G4RegionStore* regionStore = G4RegionStore::GetInstance();
268  const G4Region* world =
269  regionStore->GetRegion("DefaultRegionForTheWorld", false);
270 
271  // Identify the list of regions with different set of models
272  nRegions = 1;
273  std::vector<const G4Region*> setr;
274  setr.push_back(world);
275  G4bool isWorld = false;
276 
277  for (G4int ii=0; ii<nEmModels; ++ii) {
278  const G4Region* r = regions[ii];
279  if ( r == 0 || r == world) {
280  isWorld = true;
281  regions[ii] = world;
282  } else {
283  G4bool newRegion = true;
284  if (nRegions>1) {
285  for (G4int j=1; j<nRegions; ++j) {
286  if ( r == setr[j] ) { newRegion = false; }
287  }
288  }
289  if (newRegion) {
290  setr.push_back(r);
291  nRegions++;
292  }
293  }
294  }
295  // Are models defined?
296  if(!isWorld) {
298  ed << "No models defined for the World volume for "
299  << p->GetParticleName() << " !";
300  G4Exception("G4EmModelManager::Initialise","em0002",
301  FatalException, ed);
302  }
303 
304  G4ProductionCutsTable* theCoupleTable=
306  size_t numOfCouples = theCoupleTable->GetTableSize();
307 
308  // prepare vectors, shortcut for the case of only 1 model
309  // or only one region
310  if(nRegions > 1 && nEmModels > 1) {
311  idxOfRegionModels.resize(numOfCouples,0);
312  setOfRegionModels.resize((size_t)nRegions,0);
313  } else {
314  idxOfRegionModels.resize(1,0);
315  setOfRegionModels.resize(1,0);
316  }
317 
318  std::vector<G4int> modelAtRegion(nEmModels);
319  std::vector<G4int> modelOrd(nEmModels);
320  G4DataVector eLow(nEmModels+1);
321  G4DataVector eHigh(nEmModels);
322 
323  if(1 < verboseLevel) {
324  G4cout << " Nregions= " << nRegions
325  << " Nmodels= " << nEmModels << G4endl;
326  }
327 
328  // Order models for regions
329  for (G4int reg=0; reg<nRegions; ++reg) {
330  const G4Region* region = setr[reg];
331  G4int n = 0;
332 
333  for (G4int ii=0; ii<nEmModels; ++ii) {
334 
335  G4VEmModel* model = models[ii];
336  if ( region == regions[ii] ) {
337 
338  G4double tmin = model->LowEnergyLimit();
339  G4double tmax = model->HighEnergyLimit();
340  G4int ord = orderOfModels[ii];
341  G4bool push = true;
342  G4bool insert = false;
343  G4int idx = n;
344 
345  if(1 < verboseLevel) {
346  G4cout << "Model #" << ii
347  << " <" << model->GetName() << "> for region <";
348  if (region) G4cout << region->GetName();
349  G4cout << "> "
350  << " tmin(MeV)= " << tmin/MeV
351  << "; tmax(MeV)= " << tmax/MeV
352  << "; order= " << ord
353  << "; tminAct= " << model->LowEnergyActivationLimit()/MeV
354  << "; tmaxAct= " << model->HighEnergyActivationLimit()/MeV
355  << G4endl;
356  }
357 
358  static const G4double limitdelta = 0.01*eV;
359  if(n > 0) {
360 
361  // extend energy range to previous models
362  tmin = std::min(tmin, eHigh[n-1]);
363  tmax = std::max(tmax, eLow[0]);
364  //G4cout << "tmin= " << tmin << " tmax= "
365  // << tmax << " ord= " << ord <<G4endl;
366  // empty energy range
367  if( tmax - tmin <= limitdelta) { push = false; }
368  // low-energy model
369  else if (tmax == eLow[0]) {
370  push = false;
371  insert = true;
372  idx = 0;
373  // resolve intersections
374  } else if(tmin < eHigh[n-1]) {
375  // compare order
376  for(G4int k=0; k<n; ++k) {
377  // new model has higher order parameter,
378  // so, its application area may be reduced
379  // to avoid intersections
380  if(ord >= modelOrd[k]) {
381  if(tmin < eHigh[k] && tmin >= eLow[k]) { tmin = eHigh[k]; }
382  if(tmax <= eHigh[k] && tmax > eLow[k]) { tmax = eLow[k]; }
383  if(tmax > eHigh[k] && tmin < eLow[k]) {
384  if(tmax - eHigh[k] > eLow[k] - tmin) { tmin = eHigh[k]; }
385  else { tmax = eLow[k]; }
386  }
387  if( tmax - tmin <= limitdelta) {
388  push = false;
389  break;
390  }
391  }
392  }
393  // this model has lower order parameter than possible
394  // other models, with which there may be intersections
395  // so, appliction area of such models may be reduced
396  //G4cout << "tmin= " << tmin << " tmax= "
397  // << tmax << " push= " << push << " idx= " << idx <<G4endl;
398  //G4cout << "n= " << n << " eLow[0]= " << eLow[0] << " eLow[n-1]= " << eLow[n-1]
399  // << " eHigh[0]= " << eHigh[0] << " eHigh[n-1]= " << eHigh[n-1] << G4endl;
400 
401  // insert below the first model
402  if (tmax <= eLow[0]) {
403  push = false;
404  insert = true;
405  idx = 0;
406  // resolve intersections
407  } else if(tmin < eHigh[n-1]) {
408  // last energy interval
409  if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) {
410  eHigh[n-1] = tmin;
411  // first energy interval
412  } else if(tmin <= eLow[0] && tmax < eHigh[0]) {
413  eLow[0] = tmax;
414  push = false;
415  insert = true;
416  idx = 0;
417  // loop over all models
418  } else {
419  for(G4int k=n-1; k>=0; --k) {
420  if(tmin <= eLow[k] && tmax >= eHigh[k]) {
421  // full overlap exclude previous model
422  isUsed[modelAtRegion[k]] = 0;
423  idx = k;
424  if(k < n-1) {
425  // shift upper models and change index
426  for(G4int kk=k; kk<n-1; ++kk) {
427  modelAtRegion[kk] = modelAtRegion[kk+1];
428  modelOrd[kk] = modelOrd[kk+1];
429  eLow[kk] = eLow[kk+1];
430  eHigh[kk] = eHigh[kk+1];
431  }
432  ++k;
433  }
434  --n;
435  } else {
436  // partially reduce previous model area
437  if(tmin <= eLow[k] && tmax > eLow[k]) {
438  eLow[k] = tmax;
439  idx = k;
440  insert = true;
441  push = false;
442  } else if(tmin < eHigh[k] && tmax >= eHigh[k]) {
443  eHigh[k] = tmin;
444  idx = k + 1;
445  if(idx < n) {
446  insert = true;
447  push = false;
448  }
449  } else if(tmin > eLow[k] && tmax < eHigh[k]) {
450  if(eHigh[k] - tmax > tmin - eLow[k]) {
451  eLow[k] = tmax;
452  idx = k;
453  insert = true;
454  push = false;
455  } else {
456  eHigh[k] = tmin;
457  idx = k + 1;
458  if(idx < n) {
459  insert = true;
460  push = false;
461  }
462  }
463  }
464  }
465  }
466  }
467  }
468  }
469  }
470  // provide space for the new model
471  if(insert) {
472  for(G4int k=n-1; k>=idx; --k) {
473  modelAtRegion[k+1] = modelAtRegion[k];
474  modelOrd[k+1] = modelOrd[k];
475  eLow[k+1] = eLow[k];
476  eHigh[k+1] = eHigh[k];
477  }
478  }
479  //G4cout << "push= " << push << " insert= " << insert
480  // << " idx= " << idx <<G4endl;
481  // the model is added
482  if (push || insert) {
483  ++n;
484  modelAtRegion[idx] = ii;
485  modelOrd[idx] = ord;
486  eLow[idx] = tmin;
487  eHigh[idx] = tmax;
488  isUsed[ii] = 1;
489  }
490  // exclude models with zero energy range
491  for(G4int k=n-1; k>=0; --k) {
492  if(eHigh[k] - eLow[k] <= limitdelta) {
493  isUsed[modelAtRegion[k]] = 0;
494  if(k < n-1) {
495  for(G4int kk=k; kk<n-1; ++kk) {
496  modelAtRegion[kk] = modelAtRegion[kk+1];
497  modelOrd[kk] = modelOrd[kk+1];
498  eLow[kk] = eLow[kk+1];
499  eHigh[kk] = eHigh[kk+1];
500  }
501  }
502  --n;
503  }
504  }
505  }
506  }
507  eLow[0] = 0.0;
508  eLow[n] = eHigh[n-1];
509 
510  if(1 < verboseLevel) {
511  G4cout << "### New G4RegionModels set with " << n << " models for region <";
512  if (region) { G4cout << region->GetName(); }
513  G4cout << "> Elow(MeV)= ";
514  for(G4int iii=0; iii<=n; ++iii) {G4cout << eLow[iii]/MeV << " ";}
515  G4cout << G4endl;
516  }
517  G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region);
518  setOfRegionModels[reg] = rm;
519  // shortcut
520  if(1 == nEmModels) { break; }
521  }
522 
523  currRegionModel = setOfRegionModels[0];
524  currModel = models[0];
525 
526  // Access to materials and build cuts
527  size_t idx = 1;
528  if(secondaryParticle) {
529  if( secondaryParticle == G4Gamma::Gamma() ) { idx = 0; }
530  else if( secondaryParticle == G4Electron::Electron()) { idx = 1; }
531  else if( secondaryParticle == G4Positron::Positron()) { idx = 2; }
532  else { idx = 3; }
533  }
534 
535  theCuts =
536  static_cast<const G4DataVector*>(theCoupleTable->GetEnergyCutsVector(idx));
537 
538  // for the second run the check on cuts should be repeated
539  if(theCutsNew) { *theCutsNew = *theCuts; }
540 
541  if(minSubRange < 1.0) {
542  if( !theSubCuts ) { theSubCuts = new G4DataVector(); }
543  theSubCuts->resize(numOfCouples,DBL_MAX);
544  }
545 
546  // G4cout << "========Start define cuts" << G4endl;
547  // define cut values
548  for(size_t i=0; i<numOfCouples; ++i) {
549 
550  const G4MaterialCutsCouple* couple =
551  theCoupleTable->GetMaterialCutsCouple(i);
552  const G4Material* material = couple->GetMaterial();
553  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
554 
555  G4int reg = 0;
556  if(nRegions > 1 && nEmModels > 1) {
557  reg = nRegions;
558  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
559  do {--reg;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
560  idxOfRegionModels[i] = reg;
561  }
562  if(1 < verboseLevel) {
563  G4cout << "G4EmModelManager::Initialise() for "
564  << material->GetName()
565  << " indexOfCouple= " << i
566  << " indexOfRegion= " << reg
567  << G4endl;
568  }
569 
570  G4double cut = (*theCuts)[i];
571  if(secondaryParticle) {
572 
573  // compute subcut
574  if( cut < DBL_MAX && minSubRange < 1.0) {
575  G4double subcut = minSubRange*cut;
576  G4double rcut = std::min(minSubRange*pcuts->GetProductionCut(idx),
577  maxSubCutInRange);
578  G4double tcutmax =
579  theCoupleTable->ConvertRangeToEnergy(secondaryParticle,
580  material,rcut);
581  if(tcutmax < subcut) { subcut = tcutmax; }
582  (*theSubCuts)[i] = subcut;
583  }
584 
585  // note that idxOfRegionModels[] not always filled
586  G4int inn = 0;
587  G4int nnm = 1;
588  if(nRegions > 1 && nEmModels > 1) {
589  inn = idxOfRegionModels[i];
590  }
591  // check cuts and introduce upper limits
592  //G4cout << "idx= " << i << " cut(keV)= " << cut/keV << G4endl;
593  currRegionModel = setOfRegionModels[inn];
594  nnm = currRegionModel->NumberOfModels();
595 
596  //G4cout << "idx= " << i << " Nmod= " << nnm << G4endl;
597 
598  for(G4int jj=0; jj<nnm; ++jj) {
599  //G4cout << "jj= " << jj << " modidx= "
600  // << currRegionModel->ModelIndex(jj) << G4endl;
601  currModel = models[currRegionModel->ModelIndex(jj)];
602  G4double cutlim = currModel->MinEnergyCut(particle,couple);
603  if(cutlim > cut) {
604  if(!theCutsNew) { theCutsNew = new G4DataVector(*theCuts); }
605  (*theCutsNew)[i] = cutlim;
606  /*
607  G4cout << "### " << partname << " energy loss model in "
608  << material->GetName()
609  << " Cut was changed from " << cut/keV << " keV to "
610  << cutlim/keV << " keV " << " due to "
611  << currModel->GetName() << G4endl;
612  */
613  }
614  }
615  }
616  }
617  if(theCutsNew) { theCuts = theCutsNew; }
618 
619  // initialize models
620  G4int nn = 0;
621  severalModels = true;
622  for(G4int jj=0; jj<nEmModels; ++jj) {
623  if(1 == isUsed[jj]) {
624  ++nn;
625  currModel = models[jj];
626  currModel->Initialise(particle, *theCuts);
627  if(flucModels[jj]) { flucModels[jj]->InitialiseMe(particle); }
628  }
629  }
630  if(1 == nn) { severalModels = false; }
631 
632  if(1 < verboseLevel) {
633  G4cout << "G4EmModelManager for " << partname
634  << " is initialised; nRegions= " << nRegions
635  << " severalModels: " << severalModels
636  << G4endl;
637  }
638 
639  return theCuts;
640 }
641 
642 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
643 
645  const G4MaterialCutsCouple* couple,
646  G4EmTableType tType)
647 {
648  size_t i = couple->GetIndex();
649  G4double cut = (*theCuts)[i];
650  G4double emin = 0.0;
651 
652  if(fTotal == tType) { cut = DBL_MAX; }
653  else if(fSubRestricted == tType) {
654  emin = cut;
655  if(theSubCuts) { emin = (*theSubCuts)[i]; }
656  }
657 
658  if(1 < verboseLevel) {
659  G4cout << "G4EmModelManager::FillDEDXVector() for "
660  << couple->GetMaterial()->GetName()
661  << " cut(MeV)= " << cut
662  << " emin(MeV)= " << emin
663  << " Type " << tType
664  << " for " << particle->GetParticleName()
665  << G4endl;
666  }
667 
668  G4int reg = 0;
669  if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
670  const G4RegionModels* regModels = setOfRegionModels[reg];
671  G4int nmod = regModels->NumberOfModels();
672 
673  // Calculate energy losses vector
674 
675  //G4cout << "nmod= " << nmod << G4endl;
676  size_t totBinsLoss = aVector->GetVectorLength();
677  G4double del = 0.0;
678  G4int k0 = 0;
679 
680  for(size_t j=0; j<totBinsLoss; ++j) {
681 
682  G4double e = aVector->Energy(j);
683 
684  // Choose a model of energy losses
685  G4int k = 0;
686  if (nmod > 1) {
687  k = nmod;
688  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
689  do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
690  //G4cout << "k= " << k << G4endl;
691  if(k > 0 && k != k0) {
692  k0 = k;
693  G4double elow = regModels->LowEdgeEnergy(k);
694  G4double dedx1 = ComputeDEDX(models[regModels->ModelIndex(k-1)],
695  couple,elow,cut,emin);
696  G4double dedx2 = ComputeDEDX(models[regModels->ModelIndex(k)],
697  couple,elow,cut,emin);
698  del = 0.0;
699  if(dedx2 > 0.0) { del = (dedx1/dedx2 - 1.0)*elow; }
700  //G4cout << "elow= " << elow
701  // << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl;
702  }
703  }
704  G4double dedx =
705  ComputeDEDX(models[regModels->ModelIndex(k)],couple,e,cut,emin);
706  dedx *= (1.0 + del/e);
707 
708  if(2 < verboseLevel) {
709  G4cout << "Material= " << couple->GetMaterial()->GetName()
710  << " E(MeV)= " << e/MeV
711  << " dEdx(MeV/mm)= " << dedx*mm/MeV
712  << " del= " << del*mm/MeV<< " k= " << k
713  << " modelIdx= " << regModels->ModelIndex(k)
714  << G4endl;
715  }
716  if(dedx < 0.0) { dedx = 0.0; }
717  aVector->PutValue(j, dedx);
718  }
719 }
720 
721 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
722 
724  const G4MaterialCutsCouple* couple,
725  G4bool startFromNull,
726  G4EmTableType tType)
727 {
728  size_t i = couple->GetIndex();
729  G4double cut = (*theCuts)[i];
730  G4double tmax = DBL_MAX;
731  if (fSubRestricted == tType) {
732  tmax = cut;
733  if(theSubCuts) { cut = (*theSubCuts)[i]; }
734  }
735 
736  G4int reg = 0;
737  if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
738  const G4RegionModels* regModels = setOfRegionModels[reg];
739  G4int nmod = regModels->NumberOfModels();
740  if(1 < verboseLevel) {
741  G4cout << "G4EmModelManager::FillLambdaVector() for "
742  << particle->GetParticleName()
743  << " in " << couple->GetMaterial()->GetName()
744  << " Emin(MeV)= " << aVector->Energy(0)
745  << " Emax(MeV)= " << aVector->GetMaxEnergy()
746  << " cut= " << cut
747  << " Type " << tType
748  << " nmod= " << nmod
749  << " theSubCuts " << theSubCuts
750  << G4endl;
751  }
752 
753  // Calculate lambda vector
754  size_t totBinsLambda = aVector->GetVectorLength();
755  G4double del = 0.0;
756  G4int k0 = 0;
757  G4int k = 0;
758  G4VEmModel* mod = models[regModels->ModelIndex(0)];
759  for(size_t j=0; j<totBinsLambda; ++j) {
760 
761  G4double e = aVector->Energy(j);
762 
763  // Choose a model
764  if (nmod > 1) {
765  k = nmod;
766  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
767  do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
768  if(k > 0 && k != k0) {
769  k0 = k;
770  G4double elow = regModels->LowEdgeEnergy(k);
771  G4VEmModel* mod1 = models[regModels->ModelIndex(k-1)];
772  G4double xs1 = mod1->CrossSection(couple,particle,elow,cut,tmax);
773  mod = models[regModels->ModelIndex(k)];
774  G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax);
775  del = 0.0;
776  if(xs2 > 0.0) { del = (xs1/xs2 - 1.0)*elow; }
777  //G4cout << "New model k=" << k << " E(MeV)= " << e/MeV
778  // << " Elow(MeV)= " << elow/MeV << " del= " << del << G4endl;
779  }
780  }
781  G4double cross = mod->CrossSection(couple,particle,e,cut,tmax);
782  cross *= (1.0 + del/e);
783  if(fIsCrossSectionPrim == tType) { cross *= e; }
784 
785  if(j==0 && startFromNull) { cross = 0.0; }
786 
787  if(2 < verboseLevel) {
788  G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV
789  << " cross(1/mm)= " << cross*mm
790  << " del= " << del*mm << " k= " << k
791  << " modelIdx= " << regModels->ModelIndex(k)
792  << G4endl;
793  }
794  if(cross < 0.0) { cross = 0.0; }
795  aVector->PutValue(j, cross);
796  }
797 }
798 
799 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
800 
802 {
803  if(verb == 0) { return; }
804  for(G4int i=0; i<nRegions; ++i) {
805  G4RegionModels* r = setOfRegionModels[i];
806  const G4Region* reg = r->Region();
807  G4int n = r->NumberOfModels();
808  if(n > 0) {
809  G4cout << " ===== EM models for the G4Region " << reg->GetName()
810  << " ======" << G4endl;;
811  for(G4int j=0; j<n; ++j) {
812  G4VEmModel* model = models[r->ModelIndex(j)];
813  G4double emin =
814  std::max(r->LowEdgeEnergy(j),model->LowEnergyActivationLimit());
815  G4double emax =
816  std::min(r->LowEdgeEnergy(j+1),model->HighEnergyActivationLimit());
817  if(emax > emin) {
818  G4cout << std::setw(20);
819  G4cout << model->GetName() << " : Emin= "
820  << std::setw(8) << G4BestUnit(emin,"Energy")
821  << " Emax= "
822  << std::setw(8) << G4BestUnit(emax,"Energy");
823  G4PhysicsTable* table = model->GetCrossSectionTable();
824  if(table) {
825  size_t kk = table->size();
826  for(size_t k=0; k<kk; ++k) {
827  G4PhysicsVector* v = (*table)[k];
828  if(v) {
829  G4int nn = v->GetVectorLength() - 1;
830  G4cout << " Table with " << nn << " bins Emin= "
831  << std::setw(6) << G4BestUnit(v->Energy(0),"Energy")
832  << " Emax= "
833  << std::setw(6) << G4BestUnit(v->Energy(nn),"Energy");
834  break;
835  }
836  }
837  }
839  if(an) { G4cout << " " << an->GetName(); }
840  if(fluoFlag && model->DeexcitationFlag()) {
841  G4cout << " FluoActive";
842  }
843  G4cout << G4endl;
844  }
845  }
846  }
847  if(1 == nEmModels) { break; }
848  }
849  if(theCutsNew) {
850  G4cout << " ===== Limit on energy threshold has been applied "
851  << G4endl;
852  }
853 }
854 
855 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4int NumberOfRegionModels(size_t index_couple) const
G4double LowEnergyActivationLimit() const
Definition: G4VEmModel.hh:657
G4double HighEnergyActivationLimit() const
Definition: G4VEmModel.hh:650
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
const G4String & GetName() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4double GetMaxEnergy() const
static constexpr double mm
Definition: G4SIunits.hh:115
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4VEmModel * GetModel(G4int idx, G4bool ver=false)
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
G4double GetProductionCut(G4int index) const
virtual void DefineForRegion(const G4Region *)
Definition: G4VEmModel.cc:340
const char * p
Definition: xmltok.h:285
void UpdateEmModel(const G4String &model_name, G4double emin, G4double emax)
const G4String & GetName() const
Definition: G4Material.hh:178
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:619
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
G4PhysicsTable * GetCrossSectionTable()
Definition: G4VEmModel.hh:833
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
size_t GetVectorLength() const
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
Definition: G4VEmModel.cc:395
void FillLambdaVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4bool startFromNull=true, G4EmTableType t=fRestricted)
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
static const G4double reg
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)=0
static G4RegionStore * GetInstance()
string material
Definition: eplot.py:19
G4VEmModel * GetRegionModel(G4int idx, size_t index_couple)
G4GLOB_DLL std::ostream G4cout
void FillDEDXVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4EmTableType t=fRestricted)
bool G4bool
Definition: G4Types.hh:79
G4double CrossSection(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:503
void PutValue(size_t index, G4double theValue)
static constexpr double eV
Definition: G4SIunits.hh:215
const G4int n
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double Energy(size_t index) const
tuple v
Definition: test.py:18
void DumpModelList(G4int verb)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
static const G4double emax
static G4ProductionCutsTable * GetProductionCutsTable()
G4bool DeexcitationFlag() const
Definition: G4VEmModel.hh:685
static G4Positron * Positron()
Definition: G4Positron.cc:94
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4EmTableType
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4String & GetName() const
Definition: G4VEmModel.hh:802
G4double ConvertRangeToEnergy(const G4ParticleDefinition *particle, const G4Material *material, G4double range)
double G4double
Definition: G4Types.hh:76
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
const XML_Char XML_Content * model
Definition: expat.h:151
const G4String & GetName() const
const G4Material * GetMaterial() const