Geant4  10.02.p01
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 92921 2015-09-21 15:06:51Z 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;
91  for (G4int i=0; i<nModelsForRegion; ++i) {
92  theListOfModelIndexes[i] = indx[i];
93  lowKineticEnergy[i] = lowE[i];
94  }
96  theRegion = reg;
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
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] = 0;
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 >= 0 && 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 
226 const G4DataVector*
228  const G4ParticleDefinition* secondaryParticle,
229  G4double minSubRange,
230  G4int val)
231 {
232  verboseLevel = val;
233  G4String partname = p->GetParticleName();
234  // if(partname == "proton") verboseLevel = 2;
235  //else verboseLevel = 0;
236  if(1 < verboseLevel) {
237  G4cout << "G4EmModelManager::Initialise() for "
238  << partname << G4endl;
239  }
240  // Are models defined?
241  if(nEmModels < 1) {
243  ed << "No models found out for " << p->GetParticleName()
244  << " !";
245  G4Exception("G4EmModelManager::Initialise","em0002",
246  FatalException, ed);
247  }
248 
249  particle = p;
250  Clear(); // needed if run is not first
251 
252  G4RegionStore* regionStore = G4RegionStore::GetInstance();
253  const G4Region* world =
254  regionStore->GetRegion("DefaultRegionForTheWorld", false);
255 
256  // Identify the list of regions with different set of models
257  nRegions = 1;
258  std::vector<const G4Region*> setr;
259  setr.push_back(world);
260  G4bool isWorld = false;
261 
262  for (G4int ii=0; ii<nEmModels; ++ii) {
263  const G4Region* r = regions[ii];
264  if ( r == 0 || r == world) {
265  isWorld = true;
266  regions[ii] = world;
267  } else {
268  G4bool newRegion = true;
269  if (nRegions>1) {
270  for (G4int j=1; j<nRegions; ++j) {
271  if ( r == setr[j] ) { newRegion = false; }
272  }
273  }
274  if (newRegion) {
275  setr.push_back(r);
276  nRegions++;
277  }
278  }
279  }
280  // Are models defined?
281  if(!isWorld) {
283  ed << "No models defined for the World volume for "
284  << p->GetParticleName() << " !";
285  G4Exception("G4EmModelManager::Initialise","em0002",
286  FatalException, ed);
287  }
288 
289  G4ProductionCutsTable* theCoupleTable=
291  size_t numOfCouples = theCoupleTable->GetTableSize();
292 
293  // prepare vectors, shortcut for the case of only 1 model
294  // or only one region
295  if(nRegions > 1 && nEmModels > 1) {
296  idxOfRegionModels.resize(numOfCouples,0);
297  setOfRegionModels.resize((size_t)nRegions,0);
298  } else {
299  idxOfRegionModels.resize(1,0);
300  setOfRegionModels.resize(1,0);
301  }
302 
303  std::vector<G4int> modelAtRegion(nEmModels);
304  std::vector<G4int> modelOrd(nEmModels);
305  G4DataVector eLow(nEmModels+1);
306  G4DataVector eHigh(nEmModels);
307 
308  if(1 < verboseLevel) {
309  G4cout << " Nregions= " << nRegions
310  << " Nmodels= " << nEmModels << G4endl;
311  }
312 
313  // Order models for regions
314  for (G4int reg=0; reg<nRegions; ++reg) {
315  const G4Region* region = setr[reg];
316  G4int n = 0;
317 
318  for (G4int ii=0; ii<nEmModels; ++ii) {
319 
320  G4VEmModel* model = models[ii];
321  if ( region == regions[ii] ) {
322 
323  G4double tmin = model->LowEnergyLimit();
324  G4double tmax = model->HighEnergyLimit();
325  G4int ord = orderOfModels[ii];
326  G4bool push = true;
327  G4bool insert = false;
328  G4int idx = n;
329 
330  if(1 < verboseLevel) {
331  G4cout << "Model #" << ii
332  << " <" << model->GetName() << "> for region <";
333  if (region) G4cout << region->GetName();
334  G4cout << "> "
335  << " tmin(MeV)= " << tmin/MeV
336  << "; tmax(MeV)= " << tmax/MeV
337  << "; order= " << ord
338  << G4endl;
339  }
340 
341  if(n > 0) {
342 
343  // extend energy range to previous models
344  tmin = std::min(tmin, eHigh[n-1]);
345  tmax = std::max(tmax, eLow[0]);
346  //G4cout << "tmin= " << tmin << " tmax= "
347  // << tmax << " ord= " << ord <<G4endl;
348  // empty energy range
349  if( tmax - tmin <= eV) push = false;
350  // low-energy model
351  else if (tmax == eLow[0]) {
352  push = false;
353  insert = true;
354  idx = 0;
355  // resolve intersections
356  } else if(tmin < eHigh[n-1]) {
357  // compare order
358  for(G4int k=0; k<n; ++k) {
359  // new model has lower application
360  if(ord >= modelOrd[k]) {
361  if(tmin < eHigh[k] && tmin >= eLow[k]) tmin = eHigh[k];
362  if(tmax <= eHigh[k] && tmax > eLow[k]) tmax = eLow[k];
363  if(tmax > eHigh[k] && tmin < eLow[k]) {
364  if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k];
365  else tmax = eLow[k];
366  }
367  if( tmax - tmin <= eV) {
368  push = false;
369  break;
370  }
371  }
372  }
373  //G4cout << "tmin= " << tmin << " tmax= "
374  // << tmax << " push= " << push << " idx= " << idx <<G4endl;
375  if(push) {
376  if (tmax == eLow[0]) {
377  push = false;
378  insert = true;
379  idx = 0;
380  // continue resolve intersections
381  } else if(tmin < eHigh[n-1]) {
382  // last energy interval
383  if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) {
384  eHigh[n-1] = tmin;
385  // first energy interval
386  } else if(tmin <= eLow[0] && tmax < eHigh[0]) {
387  eLow[0] = tmax;
388  push = false;
389  insert = true;
390  idx = 0;
391  } else {
392  // find energy interval to replace
393  for(G4int k=0; k<n; ++k) {
394  if(tmin <= eLow[k] && tmax >= eHigh[k]) {
395  push = false;
396  modelAtRegion[k] = ii;
397  modelOrd[k] = ord;
398  isUsed[ii] = 1;
399  }
400  }
401  }
402  }
403  }
404  }
405  }
406  if(insert) {
407  for(G4int k=n-1; k>=idx; --k) {
408  modelAtRegion[k+1] = modelAtRegion[k];
409  modelOrd[k+1] = modelOrd[k];
410  eLow[k+1] = eLow[k];
411  eHigh[k+1] = eHigh[k];
412  }
413  }
414  //G4cout << "push= " << push << " insert= " << insert
415  //<< " idx= " << idx <<G4endl;
416  if (push || insert) {
417  ++n;
418  modelAtRegion[idx] = ii;
419  modelOrd[idx] = ord;
420  eLow[idx] = tmin;
421  eHigh[idx] = tmax;
422  isUsed[ii] = 1;
423  }
424  }
425  }
426  eLow[0] = 0.0;
427  eLow[n] = eHigh[n-1];
428 
429  if(1 < verboseLevel) {
430  G4cout << "New G4RegionModels set with " << n << " models for region <";
431  if (region) { G4cout << region->GetName(); }
432  G4cout << "> Elow(MeV)= ";
433  for(G4int ii=0; ii<=n; ++ii) {G4cout << eLow[ii]/MeV << " ";}
434  G4cout << G4endl;
435  }
436  G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region);
437  setOfRegionModels[reg] = rm;
438  // shortcut
439  if(1 == nEmModels) { break; }
440  }
441 
443  currModel = models[0];
444 
445  // Access to materials and build cuts
446  size_t idx = 1;
447  if(secondaryParticle) {
448  if( secondaryParticle == G4Gamma::Gamma() ) { idx = 0; }
449  else if( secondaryParticle == G4Electron::Electron()) { idx = 1; }
450  else if( secondaryParticle == G4Positron::Positron()) { idx = 2; }
451  else { idx = 3; }
452  }
453 
454  theCuts =
455  static_cast<const G4DataVector*>(theCoupleTable->GetEnergyCutsVector(idx));
456 
457  // for the second run the check on cuts should be repeated
458  if(theCutsNew) { *theCutsNew = *theCuts; }
459 
460  if(minSubRange < 1.0) {
461  if( !theSubCuts ) { theSubCuts = new G4DataVector(); }
462  theSubCuts->resize(numOfCouples,DBL_MAX);
463  }
464 
465  // G4cout << "========Start define cuts" << G4endl;
466  // define cut values
467  for(size_t i=0; i<numOfCouples; ++i) {
468 
469  const G4MaterialCutsCouple* couple =
470  theCoupleTable->GetMaterialCutsCouple(i);
471  const G4Material* material = couple->GetMaterial();
472  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
473 
474  G4int reg = 0;
475  if(nRegions > 1 && nEmModels > 1) {
476  reg = nRegions;
477  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
478  do {--reg;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
479  idxOfRegionModels[i] = reg;
480  }
481  if(1 < verboseLevel) {
482  G4cout << "G4EmModelManager::Initialise() for "
483  << material->GetName()
484  << " indexOfCouple= " << i
485  << " indexOfRegion= " << reg
486  << G4endl;
487  }
488 
489  G4double cut = (*theCuts)[i];
490  if(secondaryParticle) {
491 
492  // compute subcut
493  if( cut < DBL_MAX && minSubRange < 1.0) {
494  G4double subcut = minSubRange*cut;
495  G4double rcut = std::min(minSubRange*pcuts->GetProductionCut(idx),
497  G4double tcutmax =
498  theCoupleTable->ConvertRangeToEnergy(secondaryParticle,
499  material,rcut);
500  if(tcutmax < subcut) { subcut = tcutmax; }
501  (*theSubCuts)[i] = subcut;
502  }
503 
504  // note that idxOfRegionModels[] not always filled
505  G4int inn = 0;
506  G4int nnm = 1;
507  if(nRegions > 1 && nEmModels > 1) {
508  inn = idxOfRegionModels[i];
509  }
510  // check cuts and introduce upper limits
511  //G4cout << "idx= " << i << " cut(keV)= " << cut/keV << G4endl;
514 
515  //G4cout << "idx= " << i << " Nmod= " << nnm << G4endl;
516 
517  for(G4int jj=0; jj<nnm; ++jj) {
518  //G4cout << "jj= " << jj << " modidx= "
519  // << currRegionModel->ModelIndex(jj) << G4endl;
521  G4double cutlim = currModel->MinEnergyCut(particle,couple);
522  if(cutlim > cut) {
523  if(!theCutsNew) { theCutsNew = new G4DataVector(*theCuts); }
524  (*theCutsNew)[i] = cutlim;
525  /*
526  G4cout << "### " << partname << " energy loss model in "
527  << material->GetName()
528  << " Cut was changed from " << cut/keV << " keV to "
529  << cutlim/keV << " keV " << " due to "
530  << currModel->GetName() << G4endl;
531  */
532  }
533  }
534  }
535  }
536  if(theCutsNew) { theCuts = theCutsNew; }
537 
538  // initialize models
539  G4int nn = 0;
540  severalModels = true;
541  for(G4int jj=0; jj<nEmModels; ++jj) {
542  if(1 == isUsed[jj]) {
543  ++nn;
544  currModel = models[jj];
546  if(flucModels[jj]) { flucModels[jj]->InitialiseMe(particle); }
547  }
548  }
549  if(1 == nn) { severalModels = false; }
550 
551  if(1 < verboseLevel) {
552  G4cout << "G4EmModelManager for " << partname
553  << " is initialised; nRegions= " << nRegions
554  << " severalModels: " << severalModels
555  << G4endl;
556  }
557 
558  return theCuts;
559 }
560 
561 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
562 
564  const G4MaterialCutsCouple* couple,
565  G4EmTableType tType)
566 {
567  size_t i = couple->GetIndex();
568  G4double cut = (*theCuts)[i];
569  G4double emin = 0.0;
570 
571  if(fTotal == tType) { cut = DBL_MAX; }
572  else if(fSubRestricted == tType) {
573  emin = cut;
574  if(theSubCuts) { emin = (*theSubCuts)[i]; }
575  }
576 
577  if(1 < verboseLevel) {
578  G4cout << "G4EmModelManager::FillDEDXVector() for "
579  << couple->GetMaterial()->GetName()
580  << " cut(MeV)= " << cut
581  << " emin(MeV)= " << emin
582  << " Type " << tType
583  << " for " << particle->GetParticleName()
584  << G4endl;
585  }
586 
587  G4int reg = 0;
588  if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
589  const G4RegionModels* regModels = setOfRegionModels[reg];
590  G4int nmod = regModels->NumberOfModels();
591 
592  // Calculate energy losses vector
593 
594  //G4cout << "nmod= " << nmod << G4endl;
595  size_t totBinsLoss = aVector->GetVectorLength();
596  G4double del = 0.0;
597  G4int k0 = 0;
598 
599  for(size_t j=0; j<totBinsLoss; ++j) {
600 
601  G4double e = aVector->Energy(j);
602 
603  // Choose a model of energy losses
604  G4int k = 0;
605  if (nmod > 1) {
606  k = nmod;
607  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
608  do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
609  //G4cout << "k= " << k << G4endl;
610  if(k > 0 && k != k0) {
611  k0 = k;
612  G4double elow = regModels->LowEdgeEnergy(k);
613  G4double dedx1 = ComputeDEDX(models[regModels->ModelIndex(k-1)],
614  couple,elow,cut,emin);
615  G4double dedx2 = ComputeDEDX(models[regModels->ModelIndex(k)],
616  couple,elow,cut,emin);
617  del = 0.0;
618  if(dedx2 > 0.0) { del = (dedx1/dedx2 - 1.0)*elow; }
619  //G4cout << "elow= " << elow
620  // << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl;
621  }
622  }
623  G4double dedx =
624  ComputeDEDX(models[regModels->ModelIndex(k)],couple,e,cut,emin);
625  dedx *= (1.0 + del/e);
626 
627  if(2 < verboseLevel) {
628  G4cout << "Material= " << couple->GetMaterial()->GetName()
629  << " E(MeV)= " << e/MeV
630  << " dEdx(MeV/mm)= " << dedx*mm/MeV
631  << " del= " << del*mm/MeV<< " k= " << k
632  << " modelIdx= " << regModels->ModelIndex(k)
633  << G4endl;
634  }
635  if(dedx < 0.0) { dedx = 0.0; }
636  aVector->PutValue(j, dedx);
637  }
638 }
639 
640 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
641 
643  const G4MaterialCutsCouple* couple,
644  G4bool startFromNull,
645  G4EmTableType tType)
646 {
647  size_t i = couple->GetIndex();
648  G4double cut = (*theCuts)[i];
649  G4double tmax = DBL_MAX;
650  if (fSubRestricted == tType) {
651  tmax = cut;
652  if(theSubCuts) { cut = (*theSubCuts)[i]; }
653  }
654 
655  G4int reg = 0;
656  if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
657  const G4RegionModels* regModels = setOfRegionModels[reg];
658  G4int nmod = regModels->NumberOfModels();
659  if(1 < verboseLevel) {
660  G4cout << "G4EmModelManager::FillLambdaVector() for "
662  << " in " << couple->GetMaterial()->GetName()
663  << " Emin(MeV)= " << aVector->Energy(0)
664  << " Emax(MeV)= " << aVector->GetMaxEnergy()
665  << " cut= " << cut
666  << " Type " << tType
667  << " nmod= " << nmod
668  << " theSubCuts " << theSubCuts
669  << G4endl;
670  }
671 
672  // Calculate lambda vector
673  size_t totBinsLambda = aVector->GetVectorLength();
674  G4double del = 0.0;
675  G4int k0 = 0;
676  G4int k = 0;
677  G4VEmModel* mod = models[regModels->ModelIndex(0)];
678  for(size_t j=0; j<totBinsLambda; ++j) {
679 
680  G4double e = aVector->Energy(j);
681 
682  // Choose a model
683  if (nmod > 1) {
684  k = nmod;
685  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
686  do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
687  if(k > 0 && k != k0) {
688  k0 = k;
689  G4double elow = regModels->LowEdgeEnergy(k);
690  G4VEmModel* mod1 = models[regModels->ModelIndex(k-1)];
691  G4double xs1 = mod1->CrossSection(couple,particle,elow,cut,tmax);
692  mod = models[regModels->ModelIndex(k)];
693  G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax);
694  del = 0.0;
695  if(xs2 > 0.0) { del = (xs1/xs2 - 1.0)*elow; }
696  //G4cout << "New model k=" << k << " E(MeV)= " << e/MeV
697  // << " Elow(MeV)= " << elow/MeV << " del= " << del << G4endl;
698  }
699  }
700  G4double cross = mod->CrossSection(couple,particle,e,cut,tmax);
701  cross *= (1.0 + del/e);
702  if(fIsCrossSectionPrim == tType) { cross *= e; }
703 
704  if(j==0 && startFromNull) { cross = 0.0; }
705 
706  if(2 < verboseLevel) {
707  G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV
708  << " cross(1/mm)= " << cross*mm
709  << " del= " << del*mm << " k= " << k
710  << " modelIdx= " << regModels->ModelIndex(k)
711  << G4endl;
712  }
713  if(cross < 0.0) { cross = 0.0; }
714  aVector->PutValue(j, cross);
715  }
716 }
717 
718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
719 
721 {
722  if(verb == 0) { return; }
723  for(G4int i=0; i<nRegions; ++i) {
725  const G4Region* reg = r->Region();
726  G4int n = r->NumberOfModels();
727  if(n > 0) {
728  G4cout << " ===== EM models for the G4Region " << reg->GetName()
729  << " ======" << G4endl;;
730  for(G4int j=0; j<n; ++j) {
731  G4VEmModel* model = models[r->ModelIndex(j)];
732  G4double emin =
734  G4double emax =
736  G4cout << std::setw(20);
737  G4cout << model->GetName() << " : Emin= "
738  << std::setw(8) << G4BestUnit(emin,"Energy")
739  << " Emax= "
740  << std::setw(8) << G4BestUnit(emax,"Energy");
741  G4PhysicsTable* table = model->GetCrossSectionTable();
742  if(table) {
743  size_t kk = table->size();
744  for(size_t k=0; k<kk; ++k) {
745  G4PhysicsVector* v = (*table)[k];
746  if(v) {
747  G4int nn = v->GetVectorLength() - 1;
748  G4cout << " Table with " << nn << " bins Emin= "
749  << std::setw(6) << G4BestUnit(v->Energy(0),"Energy")
750  << " Emax= "
751  << std::setw(6) << G4BestUnit(v->Energy(nn),"Energy");
752  break;
753  }
754  }
755  }
757  if(an) { G4cout << " " << an->GetName(); }
758  if(fluoFlag && model->DeexcitationFlag()) {
759  G4cout << " FluoActive";
760  }
761  G4cout << G4endl;
762  }
763  }
764  if(1 == nEmModels) { break; }
765  }
766  if(theCutsNew) {
767  G4cout << " ===== Limit on energy threshold has been applied "
768  << G4endl;
769  }
770 }
771 
772 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double LowEnergyActivationLimit() const
Definition: G4VEmModel.hh:655
G4double HighEnergyActivationLimit() const
Definition: G4VEmModel.hh:648
G4double * lowKineticEnergy
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
static const double MeV
Definition: G4SIunits.hh:211
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
const G4String & GetName() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
std::vector< G4int > idxOfRegionModels
G4double GetMaxEnergy() const
void UpdateEmModel(const G4String &, G4double, G4double)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
G4VEmModel * GetModel(G4int, G4bool ver=false)
G4double GetProductionCut(G4int index) const
virtual void DefineForRegion(const G4Region *)
Definition: G4VEmModel.cc:332
const G4String & GetName() const
Definition: G4Material.hh:178
G4double LowEdgeEnergy(G4int n) const
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:617
const G4DataVector * theCuts
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
G4PhysicsTable * GetCrossSectionTable()
Definition: G4VEmModel.hh:826
#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:387
void FillLambdaVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4bool startFromNull=true, G4EmTableType t=fRestricted)
int G4int
Definition: G4Types.hh:78
std::vector< G4int > orderOfModels
const G4String & GetParticleName() const
static const G4double reg
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)=0
G4int ModelIndex(G4int n) const
static G4RegionStore * GetInstance()
G4int NumberOfModels() const
G4GLOB_DLL std::ostream G4cout
void FillDEDXVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4EmTableType t=fRestricted)
std::vector< G4int > isUsed
bool G4bool
Definition: G4Types.hh:79
G4RegionModels(G4int nMod, std::vector< G4int > &indx, G4DataVector &lowE, const G4Region *reg)
G4double CrossSection(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:501
const G4ParticleDefinition * particle
void PutValue(size_t index, G4double theValue)
const G4int n
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double Energy(size_t index) const
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
G4DataVector * theSubCuts
std::vector< G4VEmModel * > models
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:212
G4bool DeexcitationFlag() const
Definition: G4VEmModel.hh:683
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
G4double ComputeDEDX(G4VEmModel *model, const G4MaterialCutsCouple *, G4double kinEnergy, G4double cutEnergy, G4double minEnergy)
std::vector< const G4Region * > regions
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
G4RegionModels * currRegionModel
const G4String & GetName() const
Definition: G4VEmModel.hh:795
G4double ConvertRangeToEnergy(const G4ParticleDefinition *particle, const G4Material *material, G4double range)
G4VEmModel * currModel
double G4double
Definition: G4Types.hh:76
const G4Region * theRegion
const G4Region * Region() const
G4int * theListOfModelIndexes
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
std::vector< G4VEmFluctuationModel * > flucModels
static const double mm
Definition: G4SIunits.hh:114
const G4String & GetName() const
const G4Material * GetMaterial() const
G4DataVector * theCutsNew
std::vector< G4RegionModels * > setOfRegionModels