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