Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4EmModelManager.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/utils/src/G4EmModelManager.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4EmModelManager.cc (Version 9.4.p2)


  1 //                                             <<   1 //
  2 // ******************************************* <<   2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * License and Disclaimer                                           *
  4 // *                                           <<   4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License <<   7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.      <<   9 // * include a list of copyright holders.                             *
 10 // *                                           <<  10 // *                                                                  *
 11 // * Neither the authors of this software syst <<  11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin <<  12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran <<  13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum <<  14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio <<  16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                           <<  17 // *                                                                  *
 18 // * This  code  implementation is the result  <<  18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio <<  19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri <<  20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag <<  21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati <<  22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof <<  23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ******************************************* <<  24 // ********************************************************************
 25 //                                             <<  25 //
 26 //                                             <<  26 // $Id: G4EmModelManager.cc,v 1.63 2010-10-15 10:22:13 vnivanch Exp $
 27 // ------------------------------------------- <<  27 // GEANT4 tag $Name: geant4-09-04-patch-02 $
 28 //                                             <<  28 //
 29 // GEANT4 Class file                           <<  29 // -------------------------------------------------------------------
 30 //                                             <<  30 //
 31 //                                             <<  31 // GEANT4 Class file
 32 // File name:     G4EmModelManager             <<  32 //
 33 //                                             <<  33 //
 34 // Author:        Vladimir Ivanchenko          <<  34 // File name:     G4EmModelManager
 35 //                                             <<  35 //
 36 // Creation date: 07.05.2002                   <<  36 // Author:        Vladimir Ivanchenko
 37 //                                             <<  37 //
 38 // Modifications: V.Ivanchenko                 <<  38 // Creation date: 07.05.2002
 39 //                                             <<  39 //
 40 // Class Description:                          <<  40 // Modifications:
 41 //                                             <<  41 //
 42 // It is the unified energy loss process it ca <<  42 // 23-12-02 V.Ivanchenko change interface in order to move
 43 // energy loss for charged particles using a s <<  43 //                           to cut per region
 44 // models valid for different energy regions.  <<  44 // 20-01-03 Migrade to cut per region (V.Ivanchenko)
 45 // to create and access to dE/dx and range tab <<  45 // 24-01-03 Make models region aware (V.Ivanchenko)
 46 // that information on fly.                    <<  46 // 13-02-03 The set of models is defined for region (V.Ivanchenko)
 47 // ------------------------------------------- <<  47 // 06-03-03 Fix in energy intervals for models (V.Ivanchenko)
 48 //                                             <<  48 // 13-04-03 Add startFromNull (V.Ivanchenko)
 49 //....oooOO0OOooo........oooOO0OOooo........oo <<  49 // 13-05-03 Add calculation of precise range (V.Ivanchenko)
 50 //....oooOO0OOooo........oooOO0OOooo........oo <<  50 // 16-07-03 Replace G4Material by G4MaterialCutCouple in dE/dx and CrossSection
 51                                                <<  51 //          calculation (V.Ivanchenko)
 52 #include "G4EmModelManager.hh"                 <<  52 // 21-07-03 Add UpdateEmModel method (V.Ivanchenko)
 53 #include "G4SystemOfUnits.hh"                  <<  53 // 03-11-03 Substitute STL vector for G4RegionModels (V.Ivanchenko)
 54 #include "G4PhysicsTable.hh"                   <<  54 // 26-01-04 Fix in energy range conditions (V.Ivanchenko)
 55 #include "G4PhysicsVector.hh"                  <<  55 // 24-03-05 Remove check or IsInCharge (V.Ivanchenko)
 56 #include "G4VMscModel.hh"                      <<  56 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
 57                                                <<  57 // 18-08-05 Fix cut for e+e- pair production (V.Ivanchenko)
 58 #include "G4Step.hh"                           <<  58 // 29-11-05 Add protection for arithmetic operations with cut=DBL_MAX (V.Ivanchenko)
 59 #include "G4ParticleDefinition.hh"             <<  59 // 20-01-06 Introduce G4EmTableType and reducing number of methods (VI)
 60 #include "G4PhysicsVector.hh"                  <<  60 // 13-05-06 Add GetModel by index method (VI)
 61 #include "G4MaterialCutsCouple.hh"             <<  61 // 15-03-07 Add maxCutInRange (V.Ivanchenko)
 62 #include "G4ProductionCutsTable.hh"            <<  62 // 12-04-07 Add verbosity at destruction (V.Ivanchenko)
 63 #include "G4RegionStore.hh"                    <<  63 // 08-04-08 Fixed and simplified initialisation of G4RegionModel (VI)
 64 #include "G4Gamma.hh"                          <<  64 // 03-08-09 Create internal vectors only it is needed (VI)
 65 #include "G4Electron.hh"                       <<  65 //
 66 #include "G4Positron.hh"                       <<  66 // Class Description:
 67 #include "G4UnitsTable.hh"                     <<  67 //
 68 #include "G4DataVector.hh"                     <<  68 // It is the unified energy loss process it calculates the continuous
 69                                                <<  69 // energy loss for charged particles using a set of Energy Loss
 70 //....oooOO0OOooo........oooOO0OOooo........oo <<  70 // models valid for different energy regions. There are a possibility
 71                                                <<  71 // to create and access to dE/dx and range tables, or to calculate
 72 G4RegionModels::G4RegionModels(G4int nMod, std <<  72 // that information on fly.
 73                                G4DataVector& l <<  73 // -------------------------------------------------------------------
 74 {                                              <<  74 //
 75   nModelsForRegion      = nMod;                <<  75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 76   theListOfModelIndexes = new G4int [nModelsFo <<  76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 77   lowKineticEnergy      = new G4double [nModel <<  77 
 78   for (G4int i=0; i<nModelsForRegion; ++i) {   <<  78 #include "G4EmModelManager.hh"
 79     theListOfModelIndexes[i] = indx[i];        <<  79 
 80     lowKineticEnergy[i] = lowE[i];             <<  80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 81   }                                            <<  81 
 82   lowKineticEnergy[nModelsForRegion] = lowE[nM <<  82 G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& indx, 
 83   theRegion = reg;                             <<  83              G4DataVector& lowE, const G4Region* reg)
 84 }                                              <<  84 {
 85                                                <<  85   nModelsForRegion      = nMod;
 86 //....oooOO0OOooo........oooOO0OOooo........oo <<  86   theListOfModelIndexes = new G4int [nModelsForRegion];
 87                                                <<  87   lowKineticEnergy      = new G4double [nModelsForRegion+1];
 88 G4RegionModels::~G4RegionModels()              <<  88   for (G4int i=0; i<nModelsForRegion; ++i) {
 89 {                                              <<  89     theListOfModelIndexes[i] = indx[i];
 90   delete [] theListOfModelIndexes;             <<  90     lowKineticEnergy[i] = lowE[i];
 91   delete [] lowKineticEnergy;                  <<  91   }
 92 }                                              <<  92   lowKineticEnergy[nModelsForRegion] = lowE[nModelsForRegion];
 93                                                <<  93   theRegion = reg;
 94 //....oooOO0OOooo........oooOO0OOooo........oo <<  94 }
 95 //....oooOO0OOooo........oooOO0OOooo........oo <<  95 
 96                                                <<  96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 97 G4EmModelManager::G4EmModelManager()           <<  97 
 98 {                                              <<  98 G4RegionModels::~G4RegionModels()
 99   models.reserve(4);                           <<  99 {
100   flucModels.reserve(4);                       << 100   delete [] theListOfModelIndexes;
101   regions.reserve(4);                          << 101   delete [] lowKineticEnergy;
102   orderOfModels.reserve(4);                    << 102 }
103   isUsed.reserve(4);                           << 103 
104 }                                              << 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105                                                << 105 
106 //....oooOO0OOooo........oooOO0OOooo........oo << 106 #include "G4Step.hh"
107                                                << 107 #include "G4ParticleDefinition.hh"
108 G4EmModelManager::~G4EmModelManager()          << 108 #include "G4PhysicsVector.hh"
109 {                                              << 109 #include "G4Gamma.hh"
110   verboseLevel = 0; // no verbosity at destruc << 110 #include "G4Positron.hh"
111   Clear();                                     << 111 #include "G4MaterialCutsCouple.hh"
112   delete theCutsNew;                           << 112 #include "G4ProductionCutsTable.hh"
113 }                                              << 113 #include "G4RegionStore.hh"
114                                                << 114 #include "G4Gamma.hh"
115 //....oooOO0OOooo........oooOO0OOooo........oo << 115 #include "G4Positron.hh"
116                                                << 116 #include "G4UnitsTable.hh"
117 void G4EmModelManager::Clear()                 << 117 #include "G4DataVector.hh"
118 {                                              << 118 
119   if(1 < verboseLevel) {                       << 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120     G4cout << "G4EmModelManager::Clear()" << G << 120 
121   }                                            << 121 G4EmModelManager::G4EmModelManager():
122   std::size_t n = setOfRegionModels.size();    << 122   nEmModels(0),
123   for(std::size_t i=0; i<n; ++i) {             << 123   nRegions(0),
124     delete setOfRegionModels[i];               << 124   particle(0),
125     setOfRegionModels[i] = nullptr;            << 125   verboseLevel(0)
126   }                                            << 126 {
127 }                                              << 127   maxSubCutInRange = 0.7*mm;
128                                                << 128   models.reserve(4);
129 //....oooOO0OOooo........oooOO0OOooo........oo << 129   flucModels.reserve(4);
130                                                << 130   regions.reserve(4);
131 void G4EmModelManager::AddEmModel(G4int num, G << 131   orderOfModels.reserve(4);
132                                   G4VEmFluctua << 132   isUsed.reserve(4);
133 {                                              << 133   severalModels = true;
134   if(nullptr == p) {                           << 134   currRegionModel = 0;
135     G4cout << "G4EmModelManager::AddEmModel WA << 135   currModel = 0;
136            << G4endl;                          << 136 }
137     return;                                    << 137 
138   }                                            << 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139   models.push_back(p);                         << 139 
140   flucModels.push_back(fm);                    << 140 G4EmModelManager::~G4EmModelManager()
141   regions.push_back(r);                        << 141 {
142   orderOfModels.push_back(num);                << 142   verboseLevel = 0; // no verbosity at destruction
143   isUsed.push_back(0);                         << 143   Clear();
144   p->DefineForRegion(r);                       << 144 }
145   ++nEmModels;                                 << 145 
146 }                                              << 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147                                                << 147 
148 //....oooOO0OOooo........oooOO0OOooo........oo << 148 void G4EmModelManager::Clear()
149                                                << 149 {
150 G4VEmModel* G4EmModelManager::GetModel(G4int i << 150   if(1 < verboseLevel) {
151 {                                              << 151     G4cout << "G4EmModelManager::Clear()" << G4endl;
152   G4VEmModel* model = nullptr;                 << 152   }
153   if(idx >= 0 && idx < nEmModels) { model = mo << 153   size_t n = setOfRegionModels.size();
154   else if(verboseLevel > 0 && ver) {           << 154   if(n > 0) {
155     G4cout << "G4EmModelManager::GetModel WARN << 155     for(size_t i=0; i<n; ++i) {
156            << "index " << idx << " is wrong Nm << 156       delete setOfRegionModels[i];
157            << nEmModels;                       << 157       setOfRegionModels[i] = 0;
158     if(nullptr != particle) {                  << 158     }
159       G4cout << " for " << particle->GetPartic << 159   }
160     }                                          << 160 }
161     G4cout<< G4endl;                           << 161 
162   }                                            << 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
163   return model;                                << 163 
164 }                                              << 164 void G4EmModelManager::AddEmModel(G4int num, G4VEmModel* p,
165                                                << 165                                   G4VEmFluctuationModel* fm, const G4Region* r)
166 //....oooOO0OOooo........oooOO0OOooo........oo << 166 {
167                                                << 167   if(!p) {
168 G4VEmModel* G4EmModelManager::GetRegionModel(G << 168     G4cout << "G4EmModelManager::AddEmModel WARNING: no model defined." 
169 {                                              << 169      << G4endl;
170   G4RegionModels* rm = setOfRegionModels[idxOf << 170     return;
171   return (k < rm->NumberOfModels()) ? models[r << 171   }
172 }                                              << 172   models.push_back(p);
173                                                << 173   flucModels.push_back(fm);
174 //....oooOO0OOooo........oooOO0OOooo........oo << 174   regions.push_back(r);
175                                                << 175   orderOfModels.push_back(num);
176 G4int G4EmModelManager::NumberOfRegionModels(s << 176   isUsed.push_back(0);
177 {                                              << 177   p->DefineForRegion(r);
178   G4RegionModels* rm = setOfRegionModels[idxOf << 178   nEmModels++;
179   return rm->NumberOfModels();                 << 179 }
180 }                                              << 180 
181                                                << 181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
182 //....oooOO0OOooo........oooOO0OOooo........oo << 182 
183                                                << 183 void G4EmModelManager::UpdateEmModel(const G4String& nam, 
184 const G4DataVector*                            << 184              G4double emin, G4double emax)
185 G4EmModelManager::Initialise(const G4ParticleD << 185 {
186                              const G4ParticleD << 186   if (nEmModels > 0) {
187                              G4int verb)       << 187     for(G4int i=0; i<nEmModels; ++i) {
188 {                                              << 188       if(nam == models[i]->GetName()) {
189   verboseLevel = verb;                         << 189         models[i]->SetLowEnergyLimit(emin);
190   if(1 < verboseLevel) {                       << 190         models[i]->SetHighEnergyLimit(emax);
191     G4cout << "G4EmModelManager::Initialise()  << 191   break;
192            << p->GetParticleName() << "  Nmode << 192       }
193   }                                            << 193     }
194   // Are models defined?                       << 194   }
195   if(nEmModels < 1) {                          << 195   G4cout << "G4EmModelManager::UpdateEmModel WARNING: no model <"
196     G4ExceptionDescription ed;                 << 196          << nam << "> is found out"
197     ed << "No models found out for " << p->Get << 197    << G4endl;
198        << " !";                                << 198 }
199     G4Exception("G4EmModelManager::Initialise" << 199 
200                 FatalException, ed);           << 200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
201   }                                            << 201 
202                                                << 202 G4VEmModel* G4EmModelManager::GetModel(G4int i, G4bool ver)
203   particle = p;                                << 203 {
204   Clear(); // needed if run is not first       << 204   G4VEmModel* m = 0;
205                                                << 205   if(i >= 0 && i < nEmModels) {m = models[i];}
206   G4RegionStore* regionStore = G4RegionStore:: << 206   else if(verboseLevel > 0 && ver) { 
207   const G4Region* world =                      << 207     G4cout << "G4EmModelManager::GetModel WARNING: "
208     regionStore->GetRegion("DefaultRegionForTh << 208      << "index " << i << " is wrong Nmodels= "
209                                                << 209      << nEmModels;
210   // Identify the list of regions with differe << 210     if(particle) G4cout << " for " << particle->GetParticleName(); 
211   nRegions = 1;                                << 211     G4cout<< G4endl;
212   std::vector<const G4Region*> setr;           << 212   }
213   setr.push_back(world);                       << 213   return m;
214   G4bool isWorld = false;                      << 214 }
215                                                << 215 
216   for (G4int ii=0; ii<nEmModels; ++ii) {       << 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
217     const G4Region* r = regions[ii];           << 217 
218     if ( r == nullptr || r == world) {         << 218 const G4DataVector* 
219       isWorld = true;                          << 219 G4EmModelManager::Initialise(const G4ParticleDefinition* p,
220       regions[ii] = world;                     << 220            const G4ParticleDefinition* secondaryParticle,
221     } else {                                   << 221            G4double minSubRange,
222       G4bool newRegion = true;                 << 222            G4int val)
223       if (nRegions>1) {                        << 223 {
224         for (G4int j=1; j<nRegions; ++j) {     << 224   verboseLevel = val;
225           if ( r == setr[j] ) { newRegion = fa << 225   G4String partname = p->GetParticleName();
226         }                                      << 226   if(1 < verboseLevel) {
227       }                                        << 227     G4cout << "G4EmModelManager::Initialise() for "
228       if (newRegion) {                         << 228            << partname << G4endl;
229         setr.push_back(r);                     << 229   }
230         ++nRegions;                            << 230   // Are models defined?
231       }                                        << 231   if(nEmModels < 1) {
232     }                                          << 232     G4Exception("G4EmModelManager::Initialise: no model defined for " + partname);
233   }                                            << 233   }
234   // Are models defined?                       << 234 
235   if(!isWorld) {                               << 235   particle = p;
236     G4ExceptionDescription ed;                 << 236   Clear(); // needed if run is not first
237     ed << "No models defined for the World vol << 237 
238        << p->GetParticleName() << " !";        << 238   G4RegionStore* regionStore = G4RegionStore::GetInstance();
239     G4Exception("G4EmModelManager::Initialise" << 239   const G4Region* world = 
240                 FatalException, ed);           << 240     regionStore->GetRegion("DefaultRegionForTheWorld", false);
241   }                                            << 241 
242                                                << 242   // Identify the list of regions with different set of models
243   G4ProductionCutsTable* theCoupleTable=       << 243   nRegions = 1;
244     G4ProductionCutsTable::GetProductionCutsTa << 244   std::vector<const G4Region*> setr;
245   std::size_t numOfCouples = theCoupleTable->G << 245   setr.push_back(world);
246                                                << 246   G4bool isWorld = false;
247   // prepare vectors, shortcut for the case of << 247 
248   // or only one region                        << 248   for (G4int ii=0; ii<nEmModels; ++ii) {
249   if(nRegions > 1 && nEmModels > 1) {          << 249     const G4Region* r = regions[ii];
250     idxOfRegionModels.resize(numOfCouples,0);  << 250     if ( r == 0 || r == world) {
251     setOfRegionModels.resize((std::size_t)nReg << 251       isWorld = true;
252   } else {                                     << 252       regions[ii] = world;
253     idxOfRegionModels.resize(1,0);             << 253     } else {
254     setOfRegionModels.resize(1,nullptr);       << 254       G4bool newRegion = true;
255   }                                            << 255       if (nRegions>1) {
256                                                << 256         for (G4int j=1; j<nRegions; ++j) {
257   std::vector<G4int>    modelAtRegion(nEmModel << 257     if ( r == setr[j] ) newRegion = false;
258   std::vector<G4int>    modelOrd(nEmModels);   << 258         }
259   G4DataVector          eLow(nEmModels+1);     << 259       }
260   G4DataVector          eHigh(nEmModels);      << 260       if (newRegion) {
261                                                << 261         setr.push_back(r);
262   if(1 < verboseLevel) {                       << 262   nRegions++;
263     G4cout << "    Nregions= " << nRegions     << 263       }
264            << "  Nmodels= " << nEmModels << G4 << 264     }
265   }                                            << 265   }
266                                                << 266   // Are models defined?
267   // Order models for regions                  << 267   if(!isWorld) {
268   for (G4int reg=0; reg<nRegions; ++reg) {     << 268     G4Exception("G4EmModelManager::Initialise: no models defined for " + 
269     const G4Region* region = setr[reg];        << 269     partname + " in the World volume");
270     G4int n = 0;                               << 270   }
271                                                << 271 
272     for (G4int ii=0; ii<nEmModels; ++ii) {     << 272   G4ProductionCutsTable* theCoupleTable=
273                                                << 273     G4ProductionCutsTable::GetProductionCutsTable();
274       G4VEmModel* model = models[ii];          << 274   size_t numOfCouples = theCoupleTable->GetTableSize();
275       if ( region == regions[ii] ) {           << 275 
276                                                << 276   // prepare vectors, shortcut for the case of only 1 model
277         G4double tmin = model->LowEnergyLimit( << 277   if(nRegions > 1 && nEmModels > 1) {
278         G4double tmax = model->HighEnergyLimit << 278     if(numOfCouples > idxOfRegionModels.size()) {
279         G4int  ord    = orderOfModels[ii];     << 279       idxOfRegionModels.resize(numOfCouples);
280         G4bool push   = true;                  << 280     }
281         G4bool insert = false;                 << 281   }
282         G4int  idx    = n;                     << 282   size_t nr = 1;
283                                                << 283   if(nEmModels > 1) { nr = nRegions; }
284         if(1 < verboseLevel) {                 << 284   if(nr > setOfRegionModels.size()) { setOfRegionModels.resize(nr); }
285           G4cout << "Model #" << ii            << 285 
286                  << " <" << model->GetName() < << 286   std::vector<G4int>    modelAtRegion(nEmModels);
287           if (region) G4cout << region->GetNam << 287   std::vector<G4int>    modelOrd(nEmModels);
288           G4cout << ">  "                      << 288   G4DataVector          eLow(nEmModels+1);
289                  << " tmin(MeV)= " << tmin/MeV << 289   G4DataVector          eHigh(nEmModels);
290                  << "; tmax(MeV)= " << tmax/Me << 290 
291                  << "; order= " << ord         << 291   // Order models for regions
292      << "; tminAct= " << model->LowEnergyActiv << 292   for (G4int reg=0; reg<nRegions; ++reg) {
293      << "; tmaxAct= " << model->HighEnergyActi << 293     const G4Region* region = setr[reg];
294                  << G4endl;                    << 294     G4int n = 0;
295         }                                      << 295 
296                                                << 296     for (G4int ii=0; ii<nEmModels; ++ii) {
297   static const G4double limitdelta = 0.01*eV;  << 297 
298         if(n > 0) {                            << 298       G4VEmModel* model = models[ii];
299                                                << 299       if ( region == regions[ii] ) {
300           // extend energy range to previous m << 300 
301           tmin = std::min(tmin, eHigh[n-1]);   << 301   G4double tmin = model->LowEnergyLimit();
302           tmax = std::max(tmax, eLow[0]);      << 302   G4double tmax = model->HighEnergyLimit();
303           //G4cout << "tmin= " << tmin << "  t << 303   G4int  ord    = orderOfModels[ii];
304           //           << tmax << "  ord= " << << 304   G4bool push   = true;
305           // empty energy range                << 305   G4bool insert = false;
306           if( tmax - tmin <= limitdelta) { pus << 306   G4int  idx    = n;
307           // low-energy model                  << 307 
308           else if (tmax == eLow[0]) {          << 308   if(1 < verboseLevel) {
309             push = false;                      << 309     G4cout << "Model #" << ii 
310             insert = true;                     << 310      << " <" << model->GetName() << "> for region <";
311             idx = 0;                           << 311     if (region) G4cout << region->GetName();
312             // resolve intersections           << 312     G4cout << ">  "
313           } else if(tmin < eHigh[n-1]) {       << 313      << " tmin(MeV)= " << tmin/MeV
314             // compare order                   << 314      << "; tmax(MeV)= " << tmax/MeV
315             for(G4int k=0; k<n; ++k) {         << 315      << "; order= " << ord
316               // new model has higher order pa << 316      << G4endl;
317         // so, its application area may be red << 317   }
318         // to avoid intersections              << 318         
319               if(ord >= modelOrd[k]) {         << 319   if(n > 0) {
320                 if(tmin < eHigh[k]  && tmin >= << 320 
321                 if(tmax <= eHigh[k] && tmax >  << 321     // extend energy range to previous models
322                 if(tmax > eHigh[k] && tmin < e << 322     tmin = std::min(tmin, eHigh[n-1]);
323                   if(tmax - eHigh[k] > eLow[k] << 323     tmax = std::max(tmax, eLow[0]);
324                   else { tmax = eLow[k]; }     << 324     //G4cout << "tmin= " << tmin << "  tmax= " 
325                 }                              << 325     //     << tmax << "  ord= " << ord <<G4endl;
326                 if( tmax - tmin <= limitdelta) << 326     // empty energy range
327                   push = false;                << 327     if( tmax - tmin <= eV) push = false;
328                   break;                       << 328     // low-energy model
329                 }                              << 329     else if (tmax == eLow[0]) {
330               }                                << 330       push = false;
331             }                                  << 331       insert = true;
332       // this model has lower order parameter  << 332       idx = 0;
333       // other models, with which there may be << 333       // resolve intersections
334       // so, appliction area of such models ma << 334     } else if(tmin < eHigh[n-1]) { 
335                                                << 335       // compare order
336       // insert below the first model          << 336       for(G4int k=0; k<n; ++k) {
337       if (tmax <= eLow[0]) {                   << 337         // new model has lower application 
338         push = false;                          << 338         if(ord >= modelOrd[k]) {
339         insert = true;                         << 339     if(tmin < eHigh[k]  && tmin >= eLow[k]) tmin = eHigh[k];
340         idx = 0;                               << 340     if(tmax <= eHigh[k] && tmax >  eLow[k]) tmax = eLow[k];
341         // resolve intersections               << 341     if(tmax > eHigh[k] && tmin < eLow[k]) {
342       } else if(tmin < eHigh[n-1]) {           << 342       if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k];
343         // last energy interval                << 343       else tmax = eLow[k];
344         if(tmin > eLow[n-1] && tmax >= eHigh[n << 344     }
345     eHigh[n-1] = tmin;                         << 345     if( tmax - tmin <= eV) {
346     // first energy interval                   << 346       push = false;
347         } else if(tmin <= eLow[0] && tmax < eH << 347       break;
348     eLow[0] = tmax;                            << 348     }
349     push = false;                              << 349         }
350     insert = true;                             << 350       }
351     idx = 0;                                   << 351       //G4cout << "tmin= " << tmin << "  tmax= " 
352     // loop over all models                    << 352       //     << tmax << "  push= " << push << " idx= " << idx <<G4endl;
353         } else {                               << 353       if(push) {
354     for(G4int k=n-1; k>=0; --k) {              << 354         if (tmax == eLow[0]) {
355       if(tmin <= eLow[k] && tmax >= eHigh[k])  << 355     push = false;
356         // full overlap exclude previous model << 356     insert = true;
357         isUsed[modelAtRegion[k]] = 0;          << 357     idx = 0;
358         idx = k;                               << 358     // continue resolve intersections
359         if(k < n-1) {                          << 359         } else if(tmin < eHigh[n-1]) { 
360           // shift upper models and change ind << 360     // last energy interval
361           for(G4int kk=k; kk<n-1; ++kk) {      << 361     if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) {
362       modelAtRegion[kk] = modelAtRegion[kk+1]; << 362       eHigh[n-1] = tmin;
363       modelOrd[kk] = modelOrd[kk+1];           << 363       // first energy interval
364       eLow[kk] = eLow[kk+1];                   << 364     } else if(tmin <= eLow[0] && tmax < eHigh[0]) {
365       eHigh[kk] = eHigh[kk+1];                 << 365       eLow[0] = tmax;
366           }                                    << 366       push = false;
367                       ++k;                     << 367       insert = true;
368         }                                      << 368       idx = 0;
369         --n;                                   << 369     } else {
370       } else {                                 << 370       // find energy interval to replace
371         // partially reduce previous model are << 371       for(G4int k=0; k<n; ++k) { 
372         if(tmin <= eLow[k] && tmax > eLow[k])  << 372         if(tmin <= eLow[k] && tmax >= eHigh[k]) {
373           eLow[k] = tmax;                      << 373           push = false;
374                       idx = k;                 << 374           modelAtRegion[k] = ii;
375                       insert = true;           << 375           modelOrd[k] = ord;
376                       push = false;            << 376           isUsed[ii] = 1;
377         } else if(tmin < eHigh[k] && tmax >= e << 377         } 
378           eHigh[k] = tmin;                     << 378       }
379                       idx = k + 1;             << 379     }
380                       if(idx < n) {            << 380         }
381       insert = true;                           << 381       }
382       push = false;                            << 382     }
383           }                                    << 383   }
384         } else if(tmin > eLow[k] && tmax < eHi << 384   if(insert) {
385           if(eHigh[k] - tmax > tmin - eLow[k]) << 385     for(G4int k=n-1; k>=idx; --k) {    
386       eLow[k] = tmax;                          << 386       modelAtRegion[k+1] = modelAtRegion[k];
387                         idx = k;               << 387       modelOrd[k+1] = modelOrd[k];
388                         insert = true;         << 388       eLow[k+1]  = eLow[k];
389       push = false;                            << 389       eHigh[k+1] = eHigh[k];
390           } else {                             << 390     }
391       eHigh[k] = tmin;                         << 391   }
392       idx = k + 1;                             << 392   //G4cout << "push= " << push << " insert= " << insert 
393       if(idx < n) {                            << 393   //<< " idx= " << idx <<G4endl;
394         insert = true;                         << 394   if (push || insert) {
395         push = false;                          << 395     ++n;
396       }                                        << 396     modelAtRegion[idx] = ii;
397           }                                    << 397     modelOrd[idx] = ord;
398         }                                      << 398     eLow[idx]  = tmin;
399                   }                            << 399     eHigh[idx] = tmax;
400                 }                              << 400     isUsed[ii] = 1;
401               }                                << 401   }
402             }                                  << 402       }
403           }                                    << 403     }
404         }                                      << 404     eLow[0] = 0.0;
405   // provide space for the new model           << 405     eLow[n] = eHigh[n-1];
406         if(insert) {                           << 406 
407           for(G4int k=n-1; k>=idx; --k) {      << 407     if(1 < verboseLevel) {
408             modelAtRegion[k+1] = modelAtRegion << 408       G4cout << "New G4RegionModels set with " << n << " models for region <";
409             modelOrd[k+1] = modelOrd[k];       << 409       if (region) G4cout << region->GetName();
410             eLow[k+1]  = eLow[k];              << 410       G4cout << ">  Elow(MeV)= ";
411             eHigh[k+1] = eHigh[k];             << 411       for(G4int ii=0; ii<=n; ++ii) {G4cout << eLow[ii]/MeV << " ";}
412           }                                    << 412       G4cout << G4endl;
413         }                                      << 413     }
414         //G4cout << "push= " << push << " inse << 414     G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region);
415   //       << " idx= " << idx <<G4endl;        << 415     setOfRegionModels[reg] = rm;
416   // the model is added                        << 416     if(1 == nEmModels) break;
417         if (push || insert) {                  << 417   }
418           ++n;                                 << 418 
419           modelAtRegion[idx] = ii;             << 419   currRegionModel = setOfRegionModels[0];
420           modelOrd[idx] = ord;                 << 420 
421           eLow[idx]  = tmin;                   << 421   // Access to materials and build cuts
422           eHigh[idx] = tmax;                   << 422   size_t idx = 1;
423           isUsed[ii] = 1;                      << 423   if(secondaryParticle) {
424         }                                      << 424     if( secondaryParticle == G4Gamma::Gamma() )           { idx = 0; }
425   // exclude models with zero energy range     << 425     else if( secondaryParticle == G4Positron::Positron()) { idx = 2; }
426   for(G4int k=n-1; k>=0; --k) {                << 426   }
427     if(eHigh[k] - eLow[k] <= limitdelta) {     << 427 
428       isUsed[modelAtRegion[k]] = 0;            << 428   if(numOfCouples > theCuts.size()) { theCuts.resize(numOfCouples); }
429       if(k < n-1) {                            << 429   if(minSubRange < 1.0 && numOfCouples > theSubCuts.size()) {
430         for(G4int kk=k; kk<n-1; ++kk) {        << 430     theSubCuts.resize(numOfCouples);
431     modelAtRegion[kk] = modelAtRegion[kk+1];   << 431   }
432     modelOrd[kk] = modelOrd[kk+1];             << 432   for(size_t i=0; i<numOfCouples; ++i) {
433     eLow[kk] = eLow[kk+1];                     << 433 
434     eHigh[kk] = eHigh[kk+1];                   << 434     const G4MaterialCutsCouple* couple = 
435         }                                      << 435       theCoupleTable->GetMaterialCutsCouple(i);
436       }                                        << 436     const G4Material* material = couple->GetMaterial();
437       --n;                                     << 437     const G4ProductionCuts* pcuts = couple->GetProductionCuts();
438     }                                          << 438  
439   }                                            << 439     G4int reg = 0;
440       }                                        << 440     if(nRegions > 1 && nEmModels > 1) {
441     }                                          << 441       reg = nRegions;
442     eLow[0] = 0.0;                             << 442       do {--reg;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
443     eLow[n] = eHigh[n-1];                      << 443       idxOfRegionModels[i] = reg;
444                                                << 444     }
445     if(1 < verboseLevel) {                     << 445     if(1 < verboseLevel) {
446       G4cout << "### New G4RegionModels set wi << 446       G4cout << "G4EmModelManager::Initialise() for "
447              << " models for region <";        << 447        << material->GetName()
448       if (region) { G4cout << region->GetName( << 448        << " indexOfCouple= " << i
449       G4cout << ">  Elow(MeV)= ";              << 449        << " indexOfRegion= " << reg
450       for(G4int iii=0; iii<=n; ++iii) {G4cout  << 450        << G4endl;
451       G4cout << G4endl;                        << 451     }
452     }                                          << 452 
453     auto rm = new G4RegionModels(n, modelAtReg << 453     G4double cut = (*theCoupleTable->GetEnergyCutsVector(idx))[i]; 
454     setOfRegionModels[reg] = rm;               << 454     G4double subcut = DBL_MAX;
455     // shortcut                                << 455     if(secondaryParticle) {
456     if(1 == nEmModels) { break; }              << 456 
457   }                                            << 457       // compute subcut 
458                                                << 458       if( cut < DBL_MAX && minSubRange < 1.0) {
459   currRegionModel = setOfRegionModels[0];      << 459   subcut = minSubRange*cut;
460   currModel = models[0];                       << 460         G4double rcut = std::min(minSubRange*pcuts->GetProductionCut(idx), 
461                                                << 461          maxSubCutInRange);
462   // Access to materials and build cuts        << 462   G4double tcutmax = 
463   std::size_t idx = 1;                         << 463     theCoupleTable->ConvertRangeToEnergy(secondaryParticle,material,rcut);
464   if(nullptr != secondaryParticle) {           << 464   if(tcutmax < subcut) { subcut = tcutmax; }
465     if( secondaryParticle == G4Gamma::Gamma()  << 465       }
466     else if( secondaryParticle == G4Electron:: << 466     }
467     else if( secondaryParticle == G4Positron:: << 467     theCuts[i] = cut;
468     else { idx = 3; }                          << 468     if(minSubRange < 1.0) { theSubCuts[i] = subcut; }
469   }                                            << 469   }
470                                                << 470 
471   theCuts =                                    << 471   // initialize models
472     static_cast<const G4DataVector*>(theCouple << 472   G4int nn = 0;
473                                                << 473   severalModels = true;
474   // for the second run the check on cuts shou << 474   for(G4int jj=0; jj<nEmModels; ++jj) {
475   if(nullptr != theCutsNew) { *theCutsNew = *t << 475     if(1 == isUsed[jj]) {
476                                                << 476       ++nn;
477   //  G4cout << "========Start define cuts" << << 477       currModel = models[jj];
478   // define cut values                         << 478       currModel->Initialise(particle, theCuts);
479   for(std::size_t i=0; i<numOfCouples; ++i) {  << 479       if(flucModels[jj]) { flucModels[jj]->InitialiseMe(particle); }
480                                                << 480     }
481     const G4MaterialCutsCouple* couple =       << 481   }
482       theCoupleTable->GetMaterialCutsCouple((G << 482   if(1 == nn) { severalModels = false; }
483     const G4Material* material = couple->GetMa << 483 
484     const G4ProductionCuts* pcuts = couple->Ge << 484   if(1 < verboseLevel) {
485                                                << 485     G4cout << "G4EmModelManager for " << particle->GetParticleName() 
486     G4int reg = 0;                             << 486            << " is initialised; nRegions=  " << nRegions
487     if(nRegions > 1 && nEmModels > 1) {        << 487            << G4endl;
488       reg = nRegions;                          << 488   }
489       // Loop checking, 03-Aug-2015, Vladimir  << 489 
490       do {--reg;} while (reg>0 && pcuts != (se << 490   return &theCuts;
491       idxOfRegionModels[i] = reg;              << 491 }
492     }                                          << 492 
493     if(1 < verboseLevel) {                     << 493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
494       G4cout << "G4EmModelManager::Initialise( << 494 
495              << material->GetName()            << 495 void G4EmModelManager::FillDEDXVector(G4PhysicsVector* aVector,
496              << " indexOfCouple= " << i        << 496               const G4MaterialCutsCouple* couple,
497              << " indexOfRegion= " << reg      << 497                                       G4EmTableType tType)
498              << G4endl;                        << 498 {
499     }                                          << 499   size_t i = couple->GetIndex();
500                                                << 500   G4double cut = theCuts[i];
501     G4double cut = (*theCuts)[i];              << 501   G4double emin = 0.0;
502     if(nullptr != secondaryParticle) {         << 502 
503                                                << 503   if(fTotal == tType) cut = DBL_MAX;
504       // note that idxOfRegionModels[] not alw << 504   else if(fSubRestricted == tType) {
505       G4int inn = 0;                           << 505     emin = cut;
506       G4int nnm = 1;                           << 506     if(theSubCuts.size() > 0) { emin = theSubCuts[i]; }
507       if(nRegions > 1 && nEmModels > 1) {      << 507   }
508         inn = idxOfRegionModels[i];            << 508 
509       }                                        << 509   if(1 < verboseLevel) {
510       // check cuts and introduce upper limits << 510     G4cout << "G4EmModelManager::FillDEDXVector() for "
511       //G4cout << "idx= " << i << " cut(keV)=  << 511            << couple->GetMaterial()->GetName()
512       currRegionModel = setOfRegionModels[inn] << 512      << "  cut(MeV)= " << cut
513       nnm = currRegionModel->NumberOfModels(); << 513      << "  emin(MeV)= " << emin
514                                                << 514      << "  Type " << tType
515       //G4cout << "idx= " << i << " Nmod= " << << 515      << "  for " << particle->GetParticleName()
516                                                << 516            << G4endl;
517       for(G4int jj=0; jj<nnm; ++jj) {          << 517   }
518         //G4cout << "jj= " << jj << "  modidx= << 518 
519         //       << currRegionModel->ModelInde << 519   G4int reg  = 0;
520         currModel = models[currRegionModel->Mo << 520   if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
521         G4double cutlim = currModel->MinEnergy << 521   const G4RegionModels* regModels = setOfRegionModels[reg];
522         if(cutlim > cut) {                     << 522   G4int nmod = regModels->NumberOfModels();
523           if(nullptr == theCutsNew) { theCutsN << 523 
524           (*theCutsNew)[i] = cutlim;           << 524   // Calculate energy losses vector
525           /*                                   << 525 
526           G4cout << "### " << partname << " en << 526   //G4cout << "nmod= " << nmod << G4endl;
527                  << material->GetName()        << 527   size_t totBinsLoss = aVector->GetVectorLength();
528                  << "  Cut was changed from "  << 528   G4double del = 0.0;
529                  << cutlim/keV << " keV " << " << 529   G4int    k0  = 0;
530                  << currModel->GetName() << G4 << 530 
531           */                                   << 531   for(size_t j=0; j<totBinsLoss; ++j) {
532         }                                      << 532 
533       }                                        << 533     G4double e = aVector->Energy(j);
534     }                                          << 534 
535   }                                            << 535     // Choose a model of energy losses
536   if(nullptr != theCutsNew) { theCuts = theCut << 536     G4int k = 0;
537                                                << 537     if (nmod > 1) {
538   // initialize models                         << 538       k = nmod;
539   G4int nn = 0;                                << 539       do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
540   severalModels = true;                        << 540       //G4cout << "k= " << k << G4endl;
541   for(G4int jj=0; jj<nEmModels; ++jj) {        << 541       if(k > 0 && k != k0) {
542     if(1 == isUsed[jj]) {                      << 542         k0 = k;
543       ++nn;                                    << 543         G4double elow = regModels->LowEdgeEnergy(k);
544       currModel = models[jj];                  << 544   G4double dedx1 = ComputeDEDX(models[regModels->ModelIndex(k-1)],
545       currModel->Initialise(particle, *theCuts << 545              couple,elow,cut,emin);
546       if(nullptr != flucModels[jj]) { flucMode << 546   G4double dedx2 = ComputeDEDX(models[regModels->ModelIndex(k)],
547     }                                          << 547              couple,elow,cut,emin);
548   }                                            << 548         if(dedx2 > 0.0) { del = (dedx1/dedx2 - 1.0)*elow; }
549   if(1 == nn) { severalModels = false; }       << 549   //G4cout << "elow= " << elow 
550                                                << 550   //       << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl;
551   if(1 < verboseLevel) {                       << 551       }
552     G4cout << "G4EmModelManager for " << parti << 552     }
553            << " is initialised; nRegions=  " < << 553     G4double dedx = 
554            << " severalModels: " << severalMod << 554       ComputeDEDX(models[regModels->ModelIndex(k)],couple,e,cut,emin);
555            << G4endl;                          << 555     if(del > 0.0) { dedx *= (1.0 + del/e); }
556   }                                            << 556 
557   return theCuts;                              << 557     if(dedx < 0.0) { dedx = 0.0; }
558 }                                              << 558     if(2 < verboseLevel) {
559                                                << 559       G4cout << "Material= " << couple->GetMaterial()->GetName()
560 //....oooOO0OOooo........oooOO0OOooo........oo << 560        << "   E(MeV)= " << e/MeV
561                                                << 561        << "  dEdx(MeV/mm)= " << dedx*mm/MeV
562 void G4EmModelManager::FillDEDXVector(G4Physic << 562        << "  del= " << del*mm/MeV<< " k= " << k 
563                                       const G4 << 563        << " modelIdx= " << regModels->ModelIndex(k)
564                                       G4EmTabl << 564        << G4endl;
565 {                                              << 565     }
566   std::size_t i = couple->GetIndex();          << 566     aVector->PutValue(j, dedx);
567   G4double cut  = (fTotal == tType) ? DBL_MAX  << 567   }
568                                                << 568 }
569   if(1 < verboseLevel) {                       << 569 
570     G4cout << "G4EmModelManager::FillDEDXVecto << 570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
571            << couple->GetMaterial()->GetName() << 571 
572            << "  cut(MeV)= " << cut            << 572 void G4EmModelManager::FillLambdaVector(G4PhysicsVector* aVector,
573            << "  Type " << tType               << 573           const G4MaterialCutsCouple* couple,
574            << "  for " << particle->GetParticl << 574                 G4bool startFromNull,
575            << G4endl;                          << 575           G4EmTableType tType)
576   }                                            << 576 {
577                                                << 577   size_t i = couple->GetIndex();
578   G4int reg  = 0;                              << 578   G4double cut  = theCuts[i];
579   if(nRegions > 1 && nEmModels > 1) { reg = id << 579   G4double tmax = DBL_MAX;
580   const G4RegionModels* regModels = setOfRegio << 580   if (fSubRestricted == tType) {
581   G4int nmod = regModels->NumberOfModels();    << 581     tmax = cut;
582                                                << 582     if(theSubCuts.size() > 0) { cut  = theSubCuts[i]; }
583   // Calculate energy losses vector            << 583   }
584   std::size_t totBinsLoss = aVector->GetVector << 584 
585   G4double del = 0.0;                          << 585   G4int reg  = 0;
586   G4int    k0  = 0;                            << 586   if(nRegions > 1 && nEmModels > 1) { reg  = idxOfRegionModels[i]; }
587                                                << 587   const G4RegionModels* regModels = setOfRegionModels[reg];
588   for(std::size_t j=0; j<totBinsLoss; ++j) {   << 588   G4int nmod = regModels->NumberOfModels();
589     G4double e = aVector->Energy(j);           << 589   if(1 < verboseLevel) {
590                                                << 590     G4cout << "G4EmModelManager::FillLambdaVector() for "
591     // Choose a model of energy losses         << 591            << particle->GetParticleName()
592     G4int k = 0;                               << 592            << " in " << couple->GetMaterial()->GetName()
593     if (nmod > 1) {                            << 593      << " Ecut(MeV)= " << cut
594       k = nmod;                                << 594      << " Emax(MeV)= " << tmax
595       // Loop checking, 03-Aug-2015, Vladimir  << 595      << " Type " << tType   
596       do {--k;} while (k>0 && e <= regModels-> << 596      << " nmod= " << nmod
597       //G4cout << "k= " << k << G4endl;        << 597            << G4endl;
598       if(k > 0 && k != k0) {                   << 598   }
599         k0 = k;                                << 599 
600         G4double elow = regModels->LowEdgeEner << 600 
601   G4double dedx1 =                             << 601   // Calculate lambda vector
602     models[regModels->ModelIndex(k-1)]->Comput << 602   size_t totBinsLambda = aVector->GetVectorLength();
603         G4double dedx2 =                       << 603   G4double del = 0.0;
604     models[regModels->ModelIndex(k)]->ComputeD << 604   G4int    k0  = 0;
605         del = (dedx2 > 0.0) ? (dedx1/dedx2 - 1 << 605   for(size_t j=0; j<totBinsLambda; ++j) {
606         //G4cout << "elow= " << elow           << 606 
607         //       << " dedx1= " << dedx1 << " d << 607     G4double e = aVector->Energy(j);
608       }                                        << 608 
609     }                                          << 609     // Choose a model 
610     G4double dedx = (1.0 + del/e)*             << 610     G4int k = 0;
611     models[regModels->ModelIndex(k)]->ComputeD << 611     G4VEmModel* mod = models[regModels->ModelIndex(0)]; 
612                                                << 612     if (nmod > 1) {
613     if(2 < verboseLevel) {                     << 613       k = nmod;
614       G4cout << "Material= " << couple->GetMat << 614       do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
615              << "   E(MeV)= " << e/MeV         << 615       if(k > 0 && k != k0) {
616              << "  dEdx(MeV/mm)= " << dedx*mm/ << 616         k0 = k;
617              << "  del= " << del*mm/MeV<< " k= << 617         G4double elow = regModels->LowEdgeEnergy(k);
618              << " modelIdx= " << regModels->Mo << 618         G4VEmModel* m = models[regModels->ModelIndex(k-1)]; 
619              << G4endl;                        << 619   G4double xs1  = m->CrossSection(couple,particle,elow,cut,tmax);
620     }                                          << 620         mod = models[regModels->ModelIndex(k)]; 
621     dedx = std::max(dedx, 0.0);                << 621   G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax);
622     aVector->PutValue(j, dedx);                << 622         if(xs2 > 0.0) { del = (xs1 - xs2)*elow; }
623   }                                            << 623       }
624 }                                              << 624     }
625                                                << 625     G4double cross = mod->CrossSection(couple,particle,e,cut,tmax);
626 //....oooOO0OOooo........oooOO0OOooo........oo << 626     if(del > 0.0) { cross *= (1.0 + del/e); }
627                                                << 627     
628 void G4EmModelManager::FillLambdaVector(G4Phys << 628     if(j==0 && startFromNull) { cross = 0.0; }
629                                         const  << 629 
630                                         G4bool << 630     if(2 < verboseLevel) {
631                                         G4EmTa << 631       G4cout << "FillLambdaVector: " << j << ".   e(MeV)= " << e/MeV
632 {                                              << 632        << "  cross(1/mm)= " << cross*mm
633   std::size_t i = couple->GetIndex();          << 633        << " del= " << del*mm << " k= " << k 
634   G4double cut  = (*theCuts)[i];               << 634        << " modelIdx= " << regModels->ModelIndex(k)
635   G4double tmax = DBL_MAX;                     << 635        << G4endl;
636                                                << 636     }
637   G4int reg  = 0;                              << 637     if(cross < 0.0) cross = 0.0;
638   if(nRegions > 1 && nEmModels > 1) { reg  = i << 638 
639   const G4RegionModels* regModels = setOfRegio << 639     aVector->PutValue(j, cross);
640   G4int nmod = regModels->NumberOfModels();    << 640   }
641   if(1 < verboseLevel) {                       << 641 }
642     G4cout << "G4EmModelManager::FillLambdaVec << 642 
643            << particle->GetParticleName()      << 643 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
644            << " in " << couple->GetMaterial()- << 644 
645            << " Emin(MeV)= " << aVector->Energ << 645 void G4EmModelManager::DumpModelList(G4int verb)
646            << " Emax(MeV)= " << aVector->GetMa << 646 {
647            << " cut= " << cut                  << 647   if(verb == 0) { return; }
648            << " Type " << tType                << 648   for(G4int i=0; i<nRegions; ++i) {
649            << " nmod= " << nmod                << 649     G4RegionModels* r = setOfRegionModels[i];
650            << G4endl;                          << 650     const G4Region* reg = r->Region();
651   }                                            << 651     G4int n = r->NumberOfModels();  
652                                                << 652     if(n > 0) {
653   // Calculate lambda vector                   << 653       G4cout << "      ===== EM models for the G4Region  " << reg->GetName()
654   std::size_t totBinsLambda = aVector->GetVect << 654        << " ======" << G4endl;;
655   G4double del = 0.0;                          << 655       for(G4int j=0; j<n; ++j) {
656   G4int    k0  = 0;                            << 656   G4VEmModel* m = models[r->ModelIndex(j)];
657   G4int     k  = 0;                            << 657   G4cout << std::setw(20);
658   G4VEmModel* mod = models[regModels->ModelInd << 658   G4cout << m->GetName() << " :   Emin= " 
659   for(std::size_t j=0; j<totBinsLambda; ++j) { << 659          << std::setw(8) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
660                                                << 660          << "      Emax=   " 
661     G4double e = aVector->Energy(j);           << 661          << std::setw(8) << G4BestUnit(r->LowEdgeEnergy(j+1),"Energy");
662                                                << 662   G4VEmAngularDistribution* an = m->GetAngularDistribution();
663     // Choose a model                          << 663         if(an) { G4cout << "  " << an->GetName(); }
664     if (nmod > 1) {                            << 664   G4cout << G4endl;
665       k = nmod;                                << 665       }  
666       // Loop checking, 03-Aug-2015, Vladimir  << 666     }
667       do {--k;} while (k>0 && e <= regModels-> << 667     if(1 == nEmModels) { break; }
668       if(k > 0 && k != k0) {                   << 668   }
669         k0 = k;                                << 669 }
670         G4double elow = regModels->LowEdgeEner << 670 
671         G4VEmModel* mod1 = models[regModels->M << 671 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
672         G4double xs1  = mod1->CrossSection(cou << 
673         mod = models[regModels->ModelIndex(k)] << 
674         G4double xs2 = mod->CrossSection(coupl << 
675         del = (xs2 > 0.0) ? (xs1/xs2 - 1.0)*el << 
676         //G4cout << "New model k=" << k << " E << 
677         //       << " Elow(MeV)= " << elow/MeV << 
678       }                                        << 
679     }                                          << 
680     G4double cross = (1.0 + del/e)*mod->CrossS << 
681     if(fIsCrossSectionPrim == tType) { cross * << 
682                                                << 
683     if(j==0 && startFromNull) { cross = 0.0; } << 
684                                                << 
685     if(2 < verboseLevel) {                     << 
686       G4cout << "FillLambdaVector: " << j << " << 
687              << "  cross(1/mm)= " << cross*mm  << 
688              << " del= " << del*mm << " k= " < << 
689              << " modelIdx= " << regModels->Mo << 
690              << G4endl;                        << 
691     }                                          << 
692     cross = std::max(cross, 0.0);              << 
693     aVector->PutValue(j, cross);               << 
694   }                                            << 
695 }                                              << 
696                                                << 
697 //....oooOO0OOooo........oooOO0OOooo........oo << 
698                                                << 
699 void G4EmModelManager::DumpModelList(std::ostr << 
700 {                                              << 
701   if(verb == 0) { return; }                    << 
702   for(G4int i=0; i<nRegions; ++i) {            << 
703     G4RegionModels* r = setOfRegionModels[i];  << 
704     const G4Region* reg = r->Region();         << 
705     G4int n = r->NumberOfModels();             << 
706     if(n > 0) {                                << 
707       out << "      ===== EM models for the G4 << 
708     << " ======" << G4endl;                    << 
709       for(G4int j=0; j<n; ++j) {               << 
710         G4VEmModel* model = models[r->ModelInd << 
711         G4double emin =                        << 
712           std::max(r->LowEdgeEnergy(j),model-> << 
713         G4double emax =                        << 
714           std::min(r->LowEdgeEnergy(j+1),model << 
715         if(emax > emin) {                      << 
716     out << std::setw(20);                      << 
717     out << model->GetName() << " : Emin="      << 
718         << std::setw(5) << G4BestUnit(emin,"En << 
719         << " Emax="                            << 
720         << std::setw(5) << G4BestUnit(emax,"En << 
721     G4PhysicsTable* table = model->GetCrossSec << 
722     if(table) {                                << 
723       std::size_t kk = table->size();          << 
724       for(std::size_t k=0; k<kk; ++k) {        << 
725         const G4PhysicsVector* v = (*table)[k] << 
726         if(v) {                                << 
727     G4int nn = G4int(v->GetVectorLength() - 1) << 
728     out << " Nbins=" << nn << " "              << 
729         << std::setw(3) << G4BestUnit(v->Energ << 
730         << " - "                               << 
731         << std::setw(3) << G4BestUnit(v->Energ << 
732     break;                                     << 
733         }                                      << 
734       }                                        << 
735           }                                    << 
736     G4VEmAngularDistribution* an = model->GetA << 
737     if(an) { out << "  " << an->GetName(); }   << 
738     if(fluoFlag && model->DeexcitationFlag())  << 
739       out << " Fluo";                          << 
740     }                                          << 
741     out << G4endl;                             << 
742           auto msc = dynamic_cast<G4VMscModel* << 
743           if(msc != nullptr) msc->DumpParamete << 
744   }                                            << 
745       }                                        << 
746     }                                          << 
747     if(1 == nEmModels) { break; }              << 
748   }                                            << 
749   if(theCutsNew) {                             << 
750     out << "      ===== Limit on energy thresh << 
751   }                                            << 
752 }                                              << 
753                                                << 
754 //....oooOO0OOooo........oooOO0OOooo........oo << 
755                                                   672