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