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