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 8.0.p1)


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