Geant4  10.01.p02
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 84398 2014-10-15 07:20:18Z 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 = 0;
138  currModel = 0;
139  theCuts = 0;
140  theCutsNew = 0;
141  theSubCuts = 0;
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 = 0;
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  do {--reg;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
478  idxOfRegionModels[i] = reg;
479  }
480  if(1 < verboseLevel) {
481  G4cout << "G4EmModelManager::Initialise() for "
482  << material->GetName()
483  << " indexOfCouple= " << i
484  << " indexOfRegion= " << reg
485  << G4endl;
486  }
487 
488  G4double cut = (*theCuts)[i];
489  if(secondaryParticle) {
490 
491  // compute subcut
492  if( cut < DBL_MAX && minSubRange < 1.0) {
493  G4double subcut = minSubRange*cut;
494  G4double rcut = std::min(minSubRange*pcuts->GetProductionCut(idx),
496  G4double tcutmax =
497  theCoupleTable->ConvertRangeToEnergy(secondaryParticle,
498  material,rcut);
499  if(tcutmax < subcut) { subcut = tcutmax; }
500  (*theSubCuts)[i] = subcut;
501  }
502 
503  // note that idxOfRegionModels[] not always filled
504  G4int inn = 0;
505  G4int nnm = 1;
506  if(nRegions > 1 && nEmModels > 1) {
507  inn = idxOfRegionModels[i];
508  }
509  // check cuts and introduce upper limits
510  //G4cout << "idx= " << i << " cut(keV)= " << cut/keV << G4endl;
513 
514  //G4cout << "idx= " << i << " Nmod= " << nnm << G4endl;
515 
516  for(G4int jj=0; jj<nnm; ++jj) {
517  //G4cout << "jj= " << jj << " modidx= "
518  // << currRegionModel->ModelIndex(jj) << G4endl;
520  G4double cutlim = currModel->MinEnergyCut(particle,couple);
521  if(cutlim > cut) {
522  if(!theCutsNew) { theCutsNew = new G4DataVector(*theCuts); }
523  (*theCutsNew)[i] = cutlim;
524  /*
525  G4cout << "### " << partname << " energy loss model in "
526  << material->GetName()
527  << " Cut was changed from " << cut/keV << " keV to "
528  << cutlim/keV << " keV " << " due to "
529  << currModel->GetName() << G4endl;
530  */
531  }
532  }
533  }
534  }
535  if(theCutsNew) { theCuts = theCutsNew; }
536 
537  // initialize models
538  G4int nn = 0;
539  severalModels = true;
540  for(G4int jj=0; jj<nEmModels; ++jj) {
541  if(1 == isUsed[jj]) {
542  ++nn;
543  currModel = models[jj];
545  if(flucModels[jj]) { flucModels[jj]->InitialiseMe(particle); }
546  }
547  }
548  if(1 == nn) { severalModels = false; }
549 
550  if(1 < verboseLevel) {
551  G4cout << "G4EmModelManager for " << partname
552  << " is initialised; nRegions= " << nRegions
553  << " severalModels: " << severalModels
554  << G4endl;
555  }
556 
557  return theCuts;
558 }
559 
560 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
561 
563  const G4MaterialCutsCouple* couple,
564  G4EmTableType tType)
565 {
566  size_t i = couple->GetIndex();
567  G4double cut = (*theCuts)[i];
568  G4double emin = 0.0;
569 
570  if(fTotal == tType) { cut = DBL_MAX; }
571  else if(fSubRestricted == tType) {
572  emin = cut;
573  if(theSubCuts) { emin = (*theSubCuts)[i]; }
574  }
575 
576  if(1 < verboseLevel) {
577  G4cout << "G4EmModelManager::FillDEDXVector() for "
578  << couple->GetMaterial()->GetName()
579  << " cut(MeV)= " << cut
580  << " emin(MeV)= " << emin
581  << " Type " << tType
582  << " for " << particle->GetParticleName()
583  << G4endl;
584  }
585 
586  G4int reg = 0;
587  if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
588  const G4RegionModels* regModels = setOfRegionModels[reg];
589  G4int nmod = regModels->NumberOfModels();
590 
591  // Calculate energy losses vector
592 
593  //G4cout << "nmod= " << nmod << G4endl;
594  size_t totBinsLoss = aVector->GetVectorLength();
595  G4double del = 0.0;
596  G4int k0 = 0;
597 
598  for(size_t j=0; j<totBinsLoss; ++j) {
599 
600  G4double e = aVector->Energy(j);
601 
602  // Choose a model of energy losses
603  G4int k = 0;
604  if (nmod > 1) {
605  k = nmod;
606  do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
607  //G4cout << "k= " << k << G4endl;
608  if(k > 0 && k != k0) {
609  k0 = k;
610  G4double elow = regModels->LowEdgeEnergy(k);
611  G4double dedx1 = ComputeDEDX(models[regModels->ModelIndex(k-1)],
612  couple,elow,cut,emin);
613  G4double dedx2 = ComputeDEDX(models[regModels->ModelIndex(k)],
614  couple,elow,cut,emin);
615  del = 0.0;
616  if(dedx2 > 0.0) { del = (dedx1/dedx2 - 1.0)*elow; }
617  //G4cout << "elow= " << elow
618  // << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl;
619  }
620  }
621  G4double dedx =
622  ComputeDEDX(models[regModels->ModelIndex(k)],couple,e,cut,emin);
623  dedx *= (1.0 + del/e);
624 
625  if(2 < verboseLevel) {
626  G4cout << "Material= " << couple->GetMaterial()->GetName()
627  << " E(MeV)= " << e/MeV
628  << " dEdx(MeV/mm)= " << dedx*mm/MeV
629  << " del= " << del*mm/MeV<< " k= " << k
630  << " modelIdx= " << regModels->ModelIndex(k)
631  << G4endl;
632  }
633  if(dedx < 0.0) { dedx = 0.0; }
634  aVector->PutValue(j, dedx);
635  }
636 }
637 
638 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
639 
641  const G4MaterialCutsCouple* couple,
642  G4bool startFromNull,
643  G4EmTableType tType)
644 {
645  size_t i = couple->GetIndex();
646  G4double cut = (*theCuts)[i];
647  G4double tmax = DBL_MAX;
648  if (fSubRestricted == tType) {
649  tmax = cut;
650  if(theSubCuts) { cut = (*theSubCuts)[i]; }
651  }
652 
653  G4int reg = 0;
654  if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
655  const G4RegionModels* regModels = setOfRegionModels[reg];
656  G4int nmod = regModels->NumberOfModels();
657  if(1 < verboseLevel) {
658  G4cout << "G4EmModelManager::FillLambdaVector() for "
660  << " in " << couple->GetMaterial()->GetName()
661  << " Emin(MeV)= " << aVector->Energy(0)
662  << " Emax(MeV)= " << aVector->GetMaxEnergy()
663  << " cut= " << cut
664  << " Type " << tType
665  << " nmod= " << nmod
666  << " theSubCuts " << theSubCuts
667  << G4endl;
668  }
669 
670  // Calculate lambda vector
671  size_t totBinsLambda = aVector->GetVectorLength();
672  G4double del = 0.0;
673  G4int k0 = 0;
674  G4int k = 0;
675  G4VEmModel* mod = models[regModels->ModelIndex(0)];
676  for(size_t j=0; j<totBinsLambda; ++j) {
677 
678  G4double e = aVector->Energy(j);
679 
680  // Choose a model
681  if (nmod > 1) {
682  k = nmod;
683  do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
684  if(k > 0 && k != k0) {
685  k0 = k;
686  G4double elow = regModels->LowEdgeEnergy(k);
687  G4VEmModel* mod1 = models[regModels->ModelIndex(k-1)];
688  G4double xs1 = mod1->CrossSection(couple,particle,elow,cut,tmax);
689  mod = models[regModels->ModelIndex(k)];
690  G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax);
691  del = 0.0;
692  if(xs2 > 0.0) { del = (xs1/xs2 - 1.0)*elow; }
693  //G4cout << "New model k=" << k << " E(MeV)= " << e/MeV
694  // << " Elow(MeV)= " << elow/MeV << " del= " << del << G4endl;
695  }
696  }
697  G4double cross = mod->CrossSection(couple,particle,e,cut,tmax);
698  cross *= (1.0 + del/e);
699  if(fIsCrossSectionPrim == tType) { cross *= e; }
700 
701  if(j==0 && startFromNull) { cross = 0.0; }
702 
703  if(2 < verboseLevel) {
704  G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV
705  << " cross(1/mm)= " << cross*mm
706  << " del= " << del*mm << " k= " << k
707  << " modelIdx= " << regModels->ModelIndex(k)
708  << G4endl;
709  }
710  if(cross < 0.0) { cross = 0.0; }
711  aVector->PutValue(j, cross);
712  }
713 }
714 
715 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716 
718 {
719  if(verb == 0) { return; }
720  for(G4int i=0; i<nRegions; ++i) {
722  const G4Region* reg = r->Region();
723  G4int n = r->NumberOfModels();
724  if(n > 0) {
725  G4cout << " ===== EM models for the G4Region " << reg->GetName()
726  << " ======" << G4endl;;
727  for(G4int j=0; j<n; ++j) {
728  G4VEmModel* model = models[r->ModelIndex(j)];
729  G4double emin =
731  G4double emax =
733  G4cout << std::setw(20);
734  G4cout << model->GetName() << " : Emin= "
735  << std::setw(8) << G4BestUnit(emin,"Energy")
736  << " Emax= "
737  << std::setw(8) << G4BestUnit(emax,"Energy");
738  G4PhysicsTable* table = model->GetCrossSectionTable();
739  if(table) {
740  size_t kk = table->size();
741  for(size_t k=0; k<kk; ++k) {
742  G4PhysicsVector* v = (*table)[k];
743  if(v) {
744  G4int nn = v->GetVectorLength() - 1;
745  G4cout << " Table with " << nn << " bins Emin= "
746  << std::setw(6) << G4BestUnit(v->Energy(0),"Energy")
747  << " Emax= "
748  << std::setw(6) << G4BestUnit(v->Energy(nn),"Energy");
749  break;
750  }
751  }
752  }
754  if(an) { G4cout << " " << an->GetName(); }
755  if(fluoFlag && model->DeexcitationFlag()) {
756  G4cout << " FluoActive";
757  }
758  G4cout << G4endl;
759  }
760  }
761  if(1 == nEmModels) { break; }
762  }
763  if(theCutsNew) {
764  G4cout << " ===== Limit on energy threshold has been applied "
765  << G4endl;
766  }
767 }
768 
769 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double LowEnergyActivationLimit() const
Definition: G4VEmModel.hh:637
G4double HighEnergyActivationLimit() const
Definition: G4VEmModel.hh:630
G4double * lowKineticEnergy
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:623
static const double MeV
Definition: G4SIunits.hh:193
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:616
G4VEmModel * GetModel(G4int, G4bool ver=false)
G4double GetProductionCut(G4int index) const
virtual void DefineForRegion(const G4Region *)
Definition: G4VEmModel.cc:324
const G4String & GetName() const
Definition: G4Material.hh:178
G4double LowEdgeEnergy(G4int n) const
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:599
const G4DataVector * theCuts
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
G4PhysicsTable * GetCrossSectionTable()
Definition: G4VEmModel.hh:808
#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:379
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:483
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)
G4DataVector * theSubCuts
std::vector< G4VEmModel * > models
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:194
G4bool DeexcitationFlag() const
Definition: G4VEmModel.hh:665
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:777
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:102
const G4String & GetName() const
const G4Material * GetMaterial() const
G4DataVector * theCutsNew
std::vector< G4RegionModels * > setOfRegionModels