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