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