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