| Geant4 Cross Reference |
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: not supported by cvs2svn $
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 del = 0.0;
549 if(1 == nn) { severalModels = false; } << 549 if(dedx2 > 0.0) { del = (dedx1/dedx2 - 1.0)*elow; }
550 << 550 //G4cout << "elow= " << elow
551 if(1 < verboseLevel) { << 551 // << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl;
552 G4cout << "G4EmModelManager for " << parti << 552 }
553 << " is initialised; nRegions= " < << 553 }
554 << " severalModels: " << severalMod << 554 G4double dedx =
555 << G4endl; << 555 ComputeDEDX(models[regModels->ModelIndex(k)],couple,e,cut,emin);
556 } << 556 dedx *= (1.0 + del/e);
557 return theCuts; << 557
558 } << 558 if(dedx < 0.0) { dedx = 0.0; }
559 << 559 if(2 < verboseLevel) {
560 //....oooOO0OOooo........oooOO0OOooo........oo << 560 G4cout << "Material= " << couple->GetMaterial()->GetName()
561 << 561 << " E(MeV)= " << e/MeV
562 void G4EmModelManager::FillDEDXVector(G4Physic << 562 << " dEdx(MeV/mm)= " << dedx*mm/MeV
563 const G4 << 563 << " del= " << del*mm/MeV<< " k= " << k
564 G4EmTabl << 564 << " modelIdx= " << regModels->ModelIndex(k)
565 { << 565 << G4endl;
566 std::size_t i = couple->GetIndex(); << 566 }
567 G4double cut = (fTotal == tType) ? DBL_MAX << 567 aVector->PutValue(j, dedx);
568 << 568 }
569 if(1 < verboseLevel) { << 569 }
570 G4cout << "G4EmModelManager::FillDEDXVecto << 570
571 << couple->GetMaterial()->GetName() << 571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
572 << " cut(MeV)= " << cut << 572
573 << " Type " << tType << 573 void G4EmModelManager::FillLambdaVector(G4PhysicsVector* aVector,
574 << " for " << particle->GetParticl << 574 const G4MaterialCutsCouple* couple,
575 << G4endl; << 575 G4bool startFromNull,
576 } << 576 G4EmTableType tType)
577 << 577 {
578 G4int reg = 0; << 578 size_t i = couple->GetIndex();
579 if(nRegions > 1 && nEmModels > 1) { reg = id << 579 G4double cut = theCuts[i];
580 const G4RegionModels* regModels = setOfRegio << 580 G4double tmax = DBL_MAX;
581 G4int nmod = regModels->NumberOfModels(); << 581 if (fSubRestricted == tType) {
582 << 582 tmax = cut;
583 // Calculate energy losses vector << 583 if(theSubCuts.size() > 0) { cut = theSubCuts[i]; }
584 std::size_t totBinsLoss = aVector->GetVector << 584 }
585 G4double del = 0.0; << 585
586 G4int k0 = 0; << 586 G4int reg = 0;
587 << 587 if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
588 for(std::size_t j=0; j<totBinsLoss; ++j) { << 588 const G4RegionModels* regModels = setOfRegionModels[reg];
589 G4double e = aVector->Energy(j); << 589 G4int nmod = regModels->NumberOfModels();
590 << 590 if(1 < verboseLevel) {
591 // Choose a model of energy losses << 591 G4cout << "G4EmModelManager::FillLambdaVector() for "
592 G4int k = 0; << 592 << particle->GetParticleName()
593 if (nmod > 1) { << 593 << " in " << couple->GetMaterial()->GetName()
594 k = nmod; << 594 << " Ecut(MeV)= " << cut
595 // Loop checking, 03-Aug-2015, Vladimir << 595 << " Emax(MeV)= " << tmax
596 do {--k;} while (k>0 && e <= regModels-> << 596 << " Type " << tType
597 //G4cout << "k= " << k << G4endl; << 597 << " nmod= " << nmod
598 if(k > 0 && k != k0) { << 598 << G4endl;
599 k0 = k; << 599 }
600 G4double elow = regModels->LowEdgeEner << 600
601 G4double dedx1 = << 601
602 models[regModels->ModelIndex(k-1)]->Comput << 602 // Calculate lambda vector
603 G4double dedx2 = << 603 size_t totBinsLambda = aVector->GetVectorLength();
604 models[regModels->ModelIndex(k)]->ComputeD << 604 G4double del = 0.0;
605 del = (dedx2 > 0.0) ? (dedx1/dedx2 - 1 << 605 G4int k0 = 0;
606 //G4cout << "elow= " << elow << 606 G4int k = 0;
607 // << " dedx1= " << dedx1 << " d << 607 G4VEmModel* mod = models[regModels->ModelIndex(0)];
608 } << 608
609 } << 609 for(size_t j=0; j<totBinsLambda; ++j) {
610 G4double dedx = (1.0 + del/e)* << 610
611 models[regModels->ModelIndex(k)]->ComputeD << 611 G4double e = aVector->Energy(j);
612 << 612
613 if(2 < verboseLevel) { << 613 // Choose a model
614 G4cout << "Material= " << couple->GetMat << 614 if (nmod > 1) {
615 << " E(MeV)= " << e/MeV << 615 k = nmod;
616 << " dEdx(MeV/mm)= " << dedx*mm/ << 616 do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
617 << " del= " << del*mm/MeV<< " k= << 617 if(k > 0 && k != k0) {
618 << " modelIdx= " << regModels->Mo << 618 k0 = k;
619 << G4endl; << 619 G4double elow = regModels->LowEdgeEnergy(k);
620 } << 620 G4VEmModel* mod1 = models[regModels->ModelIndex(k-1)];
621 dedx = std::max(dedx, 0.0); << 621 G4double xs1 = mod1->CrossSection(couple,particle,elow,cut,tmax);
622 aVector->PutValue(j, dedx); << 622 mod = models[regModels->ModelIndex(k)];
623 } << 623 G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax);
624 } << 624 del = 0.0;
625 << 625 if(xs2 > 0.0) { del = (xs1/xs2 - 1.0)*elow; }
626 //....oooOO0OOooo........oooOO0OOooo........oo << 626 }
627 << 627 }
628 void G4EmModelManager::FillLambdaVector(G4Phys << 628 G4double cross = mod->CrossSection(couple,particle,e,cut,tmax);
629 const << 629 cross *= (1.0 + del/e);
630 G4bool << 630
631 G4EmTa << 631 if(j==0 && startFromNull) { cross = 0.0; }
632 { << 632
633 std::size_t i = couple->GetIndex(); << 633 if(2 < verboseLevel) {
634 G4double cut = (*theCuts)[i]; << 634 G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV
635 G4double tmax = DBL_MAX; << 635 << " cross(1/mm)= " << cross*mm
636 << 636 << " del= " << del*mm << " k= " << k
637 G4int reg = 0; << 637 << " modelIdx= " << regModels->ModelIndex(k)
638 if(nRegions > 1 && nEmModels > 1) { reg = i << 638 << G4endl;
639 const G4RegionModels* regModels = setOfRegio << 639 }
640 G4int nmod = regModels->NumberOfModels(); << 640 if(cross < 0.0) cross = 0.0;
641 if(1 < verboseLevel) { << 641
642 G4cout << "G4EmModelManager::FillLambdaVec << 642 aVector->PutValue(j, cross);
643 << particle->GetParticleName() << 643 }
644 << " in " << couple->GetMaterial()- << 644 }
645 << " Emin(MeV)= " << aVector->Energ << 645
646 << " Emax(MeV)= " << aVector->GetMa << 646 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
647 << " cut= " << cut << 647
648 << " Type " << tType << 648 void G4EmModelManager::DumpModelList(G4int verb)
649 << " nmod= " << nmod << 649 {
650 << G4endl; << 650 if(verb == 0) { return; }
651 } << 651 for(G4int i=0; i<nRegions; ++i) {
652 << 652 G4RegionModels* r = setOfRegionModels[i];
653 // Calculate lambda vector << 653 const G4Region* reg = r->Region();
654 std::size_t totBinsLambda = aVector->GetVect << 654 G4int n = r->NumberOfModels();
655 G4double del = 0.0; << 655 if(n > 0) {
656 G4int k0 = 0; << 656 G4cout << " ===== EM models for the G4Region " << reg->GetName()
657 G4int k = 0; << 657 << " ======" << G4endl;;
658 G4VEmModel* mod = models[regModels->ModelInd << 658 for(G4int j=0; j<n; ++j) {
659 for(std::size_t j=0; j<totBinsLambda; ++j) { << 659 G4VEmModel* m = models[r->ModelIndex(j)];
660 << 660 G4cout << std::setw(20);
661 G4double e = aVector->Energy(j); << 661 G4cout << m->GetName() << " : Emin= "
662 << 662 << std::setw(8) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
663 // Choose a model << 663 << " Emax= "
664 if (nmod > 1) { << 664 << std::setw(8) << G4BestUnit(r->LowEdgeEnergy(j+1),"Energy");
665 k = nmod; << 665 G4VEmAngularDistribution* an = m->GetAngularDistribution();
666 // Loop checking, 03-Aug-2015, Vladimir << 666 if(an) { G4cout << " " << an->GetName(); }
667 do {--k;} while (k>0 && e <= regModels-> << 667 G4cout << G4endl;
668 if(k > 0 && k != k0) { << 668 }
669 k0 = k; << 669 }
670 G4double elow = regModels->LowEdgeEner << 670 if(1 == nEmModels) { break; }
671 G4VEmModel* mod1 = models[regModels->M << 671 }
672 G4double xs1 = mod1->CrossSection(cou << 672 }
673 mod = models[regModels->ModelIndex(k)] << 673
674 G4double xs2 = mod->CrossSection(coupl << 674 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 675