| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Geant4 class G4EmUtility 27 // Geant4 class G4EmUtility
28 // 28 //
29 // Author V.Ivanchenko 14.03.2022 29 // Author V.Ivanchenko 14.03.2022
30 // 30 //
31 31
32 #include "G4EmUtility.hh" 32 #include "G4EmUtility.hh"
33 #include "G4RegionStore.hh" 33 #include "G4RegionStore.hh"
34 #include "G4ProductionCutsTable.hh" 34 #include "G4ProductionCutsTable.hh"
35 #include "G4VEmProcess.hh" 35 #include "G4VEmProcess.hh"
36 #include "G4EmParameters.hh" 36 #include "G4EmParameters.hh"
37 #include "G4PhysicsVector.hh" 37 #include "G4PhysicsVector.hh"
38 #include "Randomize.hh" 38 #include "Randomize.hh"
39 #include "G4Log.hh" 39 #include "G4Log.hh"
40 #include "G4Exp.hh" 40 #include "G4Exp.hh"
41 41
42 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 43
44 static const G4double g4log10 = G4Log(10.); 44 static const G4double g4log10 = G4Log(10.);
45 45
46 const G4Region* 46 const G4Region*
47 G4EmUtility::FindRegion(const G4String& region 47 G4EmUtility::FindRegion(const G4String& regionName, const G4int verbose)
48 { 48 {
49 const G4Region* reg = nullptr; 49 const G4Region* reg = nullptr;
50 G4RegionStore* regStore = G4RegionStore::Get 50 G4RegionStore* regStore = G4RegionStore::GetInstance();
51 G4String r = regionName; 51 G4String r = regionName;
52 if(r == "") { r = "DefaultRegionForTheWorld" 52 if(r == "") { r = "DefaultRegionForTheWorld"; }
53 reg = regStore->GetRegion(r, true); 53 reg = regStore->GetRegion(r, true);
54 if(nullptr == reg && verbose > 0) { 54 if(nullptr == reg && verbose > 0) {
55 G4cout << "### G4EmUtility WARNING: fails 55 G4cout << "### G4EmUtility WARNING: fails to find a region <"
56 << r << G4endl; 56 << r << G4endl;
57 } else if(verbose > 1) { 57 } else if(verbose > 1) {
58 G4cout << "### G4EmUtility finds out G4Reg 58 G4cout << "### G4EmUtility finds out G4Region <" << r << ">"
59 << G4endl; 59 << G4endl;
60 } 60 }
61 return reg; 61 return reg;
62 } 62 }
63 63
64 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
65 65
66 const G4Element* G4EmUtility::SampleRandomElem 66 const G4Element* G4EmUtility::SampleRandomElement(const G4Material* mat)
67 { 67 {
68 const G4Element* elm = mat->GetElement(0); 68 const G4Element* elm = mat->GetElement(0);
69 std::size_t nElements = mat->GetNumberOfElem 69 std::size_t nElements = mat->GetNumberOfElements();
70 if(1 < nElements) { 70 if(1 < nElements) {
71 G4double x = mat->GetTotNbOfElectPerVolume 71 G4double x = mat->GetTotNbOfElectPerVolume()*G4UniformRand();
72 const G4double* y = mat->GetVecNbOfAtomsPe 72 const G4double* y = mat->GetVecNbOfAtomsPerVolume();
73 for(std::size_t i=0; i<nElements; ++i) { 73 for(std::size_t i=0; i<nElements; ++i) {
74 elm = mat->GetElement((G4int)i); 74 elm = mat->GetElement((G4int)i);
75 x -= y[i]*elm->GetZ(); 75 x -= y[i]*elm->GetZ();
76 if(x <= 0.0) { break; } 76 if(x <= 0.0) { break; }
77 } 77 }
78 } 78 }
79 return elm; 79 return elm;
80 } 80 }
81 81
82 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
83 83
84 const G4Isotope* G4EmUtility::SampleRandomIsot 84 const G4Isotope* G4EmUtility::SampleRandomIsotope(const G4Element* elm)
85 { 85 {
86 const std::size_t ni = elm->GetNumberOfIsoto 86 const std::size_t ni = elm->GetNumberOfIsotopes();
87 const G4Isotope* iso = elm->GetIsotope(0); 87 const G4Isotope* iso = elm->GetIsotope(0);
88 if(ni > 1) { 88 if(ni > 1) {
89 const G4double* ab = elm->GetRelativeAbund 89 const G4double* ab = elm->GetRelativeAbundanceVector();
90 G4double x = G4UniformRand(); 90 G4double x = G4UniformRand();
91 for(std::size_t idx=0; idx<ni; ++idx) { 91 for(std::size_t idx=0; idx<ni; ++idx) {
92 x -= ab[idx]; 92 x -= ab[idx];
93 if (x <= 0.0) { 93 if (x <= 0.0) {
94 iso = elm->GetIsotope((G4int)idx); 94 iso = elm->GetIsotope((G4int)idx);
95 break; 95 break;
96 } 96 }
97 } 97 }
98 } 98 }
99 return iso; 99 return iso;
100 } 100 }
101 101
102 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
103 103
104 std::vector<G4double>* G4EmUtility::FindCrossS 104 std::vector<G4double>* G4EmUtility::FindCrossSectionMax(G4PhysicsTable* p)
105 { 105 {
106 std::vector<G4double>* ptr = nullptr; 106 std::vector<G4double>* ptr = nullptr;
107 if(nullptr == p) { return ptr; } 107 if(nullptr == p) { return ptr; }
108 108
109 const std::size_t n = p->length(); 109 const std::size_t n = p->length();
110 ptr = new std::vector<G4double>; 110 ptr = new std::vector<G4double>;
111 ptr->resize(n, DBL_MAX); 111 ptr->resize(n, DBL_MAX);
112 112
113 G4bool isPeak = false; 113 G4bool isPeak = false;
114 G4double e, ss, ee, xs; 114 G4double e, ss, ee, xs;
115 115
116 // first loop on existing vectors 116 // first loop on existing vectors
117 for (std::size_t i=0; i<n; ++i) { 117 for (std::size_t i=0; i<n; ++i) {
118 const G4PhysicsVector* pv = (*p)[i]; 118 const G4PhysicsVector* pv = (*p)[i];
119 xs = ee = 0.0; 119 xs = ee = 0.0;
120 if(nullptr != pv) { 120 if(nullptr != pv) {
121 G4int nb = (G4int)pv->GetVectorLength(); 121 G4int nb = (G4int)pv->GetVectorLength();
122 for (G4int j=0; j<nb; ++j) { 122 for (G4int j=0; j<nb; ++j) {
123 e = pv->Energy(j); 123 e = pv->Energy(j);
124 ss = (*pv)(j); 124 ss = (*pv)(j);
125 if(ss >= xs) { 125 if(ss >= xs) {
126 xs = ss; 126 xs = ss;
127 ee = e; 127 ee = e;
128 continue; 128 continue;
129 } else { 129 } else {
130 isPeak = true; 130 isPeak = true;
131 (*ptr)[i] = ee; 131 (*ptr)[i] = ee;
132 break; 132 break;
133 } 133 }
134 } 134 }
135 } 135 }
136 } 136 }
137 137
138 // there is no peak for any material 138 // there is no peak for any material
139 if(!isPeak) { 139 if(!isPeak) {
140 delete ptr; 140 delete ptr;
141 ptr = nullptr; 141 ptr = nullptr;
142 } 142 }
143 return ptr; 143 return ptr;
144 } 144 }
145 145
146 //....oooOO0OOooo........oooOO0OOooo........oo 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
147 147
148 std::vector<G4double>* 148 std::vector<G4double>*
149 G4EmUtility::FindCrossSectionMax(G4VDiscretePr 149 G4EmUtility::FindCrossSectionMax(G4VDiscreteProcess* p,
150 const G4Parti 150 const G4ParticleDefinition* part)
151 { 151 {
152 std::vector<G4double>* ptr = nullptr; 152 std::vector<G4double>* ptr = nullptr;
153 if (nullptr == p || nullptr == part) { retur << 153 if(nullptr == p || nullptr == part) { return ptr; }
154 << 154 /*
>> 155 G4cout << "G4EmUtility::FindCrossSectionMax for "
>> 156 << p->GetProcessName() << " and " << part->GetParticleName() << G4endl;
>> 157 */
155 G4EmParameters* theParameters = G4EmParamete 158 G4EmParameters* theParameters = G4EmParameters::Instance();
156 const G4double tmin = theParameters->MinKinE << 159 G4double tmin = theParameters->MinKinEnergy();
157 const G4double tmax = theParameters->MaxKinE << 160 G4double tmax = theParameters->MaxKinEnergy();
158 const G4double ee = G4Log(tmax/tmin); <<
159 const G4double scale = theParameters->Number <<
160 G4int nbin = static_cast<G4int>(ee*scale); <<
161 nbin = std::max(nbin, 4); <<
162 G4double x = G4Exp(ee/(G4double)nbin); <<
163 161
164 const G4ProductionCutsTable* theCoupleTable= 162 const G4ProductionCutsTable* theCoupleTable=
165 G4ProductionCutsTable::GetProductionCu 163 G4ProductionCutsTable::GetProductionCutsTable();
166 std::size_t n = theCoupleTable->GetTableSize 164 std::size_t n = theCoupleTable->GetTableSize();
167 ptr = new std::vector<G4double>; 165 ptr = new std::vector<G4double>;
168 ptr->resize(n, DBL_MAX); 166 ptr->resize(n, DBL_MAX);
169 167
170 G4bool isPeak = false; 168 G4bool isPeak = false;
>> 169 G4double scale = theParameters->NumberOfBinsPerDecade()/g4log10;
>> 170
>> 171 G4double e, sig, ee, x, sm, em, emin, emax;
171 172
172 // first loop on existing vectors 173 // first loop on existing vectors
173 const G4int nn = static_cast<G4int>(n); << 174 for (std::size_t i=0; i<n; ++i) {
174 for (G4int i=0; i<nn; ++i) { << 175 auto couple = theCoupleTable->GetMaterialCutsCouple((G4int)i);
175 G4double sm = 0.0; << 176 emin = std::max(p->MinPrimaryEnergy(part, couple->GetMaterial()), tmin);
176 G4double em = 0.0; << 177 emax = std::max(tmax, 2*emin);
177 G4double e = tmin; << 178 ee = G4Log(emax/emin);
178 for (G4int j=0; j<=nbin; ++j) { << 179
179 G4double sig = p->GetCrossSection(e, the << 180 G4int nbin = G4lrint(ee*scale);
180 if (sig >= sm) { << 181 if(nbin < 4) { nbin = 4; }
>> 182 x = G4Exp(ee/nbin);
>> 183 sm = 0.0;
>> 184 em = 0.0;
>> 185 e = emin;
>> 186 for(G4int j=0; j<=nbin; ++j) {
>> 187 sig = p->GetCrossSection(e, couple);
>> 188 if(sig >= sm) {
181 em = e; 189 em = e;
182 sm = sig; 190 sm = sig;
183 e = (j+1 < nbin) ? e*x : tmax; << 191 e = (j+1 < nbin) ? e*x : emax;
184 } else { 192 } else {
185 isPeak = true; 193 isPeak = true;
186 (*ptr)[i] = em; 194 (*ptr)[i] = em;
187 break; 195 break;
188 } 196 }
189 } 197 }
>> 198 //G4cout << i << ". em=" << em << " sm=" << sm << G4endl;
190 } 199 }
191 // there is no peak for any couple 200 // there is no peak for any couple
192 if (!isPeak) { << 201 if(!isPeak) {
193 delete ptr; 202 delete ptr;
194 ptr = nullptr; 203 ptr = nullptr;
195 } 204 }
196 return ptr; 205 return ptr;
197 } 206 }
198 207
199 //....oooOO0OOooo........oooOO0OOooo........oo 208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
200 209
201 std::vector<G4TwoPeaksXS*>* 210 std::vector<G4TwoPeaksXS*>*
202 G4EmUtility::FillPeaksStructure(G4PhysicsTable 211 G4EmUtility::FillPeaksStructure(G4PhysicsTable* p, G4LossTableBuilder* bld)
203 { 212 {
204 std::vector<G4TwoPeaksXS*>* ptr = nullptr; 213 std::vector<G4TwoPeaksXS*>* ptr = nullptr;
205 if(nullptr == p) { return ptr; } 214 if(nullptr == p) { return ptr; }
206 215
207 const G4int n = (G4int)p->length(); 216 const G4int n = (G4int)p->length();
208 ptr = new std::vector<G4TwoPeaksXS*>; 217 ptr = new std::vector<G4TwoPeaksXS*>;
209 ptr->resize(n, nullptr); 218 ptr->resize(n, nullptr);
210 219
211 G4double e, ss, xs, ee; 220 G4double e, ss, xs, ee;
212 G4double e1peak, e1deep, e2peak, e2deep, e3p 221 G4double e1peak, e1deep, e2peak, e2deep, e3peak;
213 G4bool isDeep = false; 222 G4bool isDeep = false;
214 223
215 // first loop on existing vectors 224 // first loop on existing vectors
216 for (G4int i=0; i<n; ++i) { 225 for (G4int i=0; i<n; ++i) {
217 const G4PhysicsVector* pv = (*p)[i]; 226 const G4PhysicsVector* pv = (*p)[i];
218 ee = xs = 0.0; 227 ee = xs = 0.0;
219 e1peak = e1deep = e2peak = e2deep = e3peak 228 e1peak = e1deep = e2peak = e2deep = e3peak = DBL_MAX;
220 if(nullptr != pv) { 229 if(nullptr != pv) {
221 G4int nb = (G4int)pv->GetVectorLength(); 230 G4int nb = (G4int)pv->GetVectorLength();
222 for (G4int j=0; j<nb; ++j) { 231 for (G4int j=0; j<nb; ++j) {
223 e = pv->Energy(j); 232 e = pv->Energy(j);
224 ss = (*pv)(j); 233 ss = (*pv)(j);
225 // find out 1st peak 234 // find out 1st peak
226 if(e1peak == DBL_MAX) { 235 if(e1peak == DBL_MAX) {
227 if(ss >= xs) { 236 if(ss >= xs) {
228 xs = ss; 237 xs = ss;
229 ee = e; 238 ee = e;
230 continue; 239 continue;
231 } else { 240 } else {
232 e1peak = ee; 241 e1peak = ee;
233 } 242 }
234 } 243 }
235 // find out the deep 244 // find out the deep
236 if(e1deep == DBL_MAX) { 245 if(e1deep == DBL_MAX) {
237 if(ss <= xs) { 246 if(ss <= xs) {
238 xs = ss; 247 xs = ss;
239 ee = e; 248 ee = e;
240 continue; 249 continue;
241 } else { 250 } else {
242 e1deep = ee; 251 e1deep = ee;
243 isDeep = true; 252 isDeep = true;
244 } 253 }
245 } 254 }
246 // find out 2nd peak 255 // find out 2nd peak
247 if(e2peak == DBL_MAX) { 256 if(e2peak == DBL_MAX) {
248 if(ss >= xs) { 257 if(ss >= xs) {
249 xs = ss; 258 xs = ss;
250 ee = e; 259 ee = e;
251 continue; 260 continue;
252 } else { 261 } else {
253 e2peak = ee; 262 e2peak = ee;
254 } 263 }
255 } 264 }
256 if(e2deep == DBL_MAX) { 265 if(e2deep == DBL_MAX) {
257 if(ss <= xs) { 266 if(ss <= xs) {
258 xs = ss; 267 xs = ss;
259 ee = e; 268 ee = e;
260 continue; 269 continue;
261 } else { 270 } else {
262 e2deep = ee; 271 e2deep = ee;
263 break; 272 break;
264 } 273 }
265 } 274 }
266 // find out 3d peak 275 // find out 3d peak
267 if(e3peak == DBL_MAX) { 276 if(e3peak == DBL_MAX) {
268 if(ss >= xs) { 277 if(ss >= xs) {
269 xs = ss; 278 xs = ss;
270 ee = e; 279 ee = e;
271 continue; 280 continue;
272 } else { 281 } else {
273 e3peak = ee; 282 e3peak = ee;
274 } 283 }
275 } 284 }
276 } 285 }
277 } 286 }
278 G4TwoPeaksXS* x = (*ptr)[i]; 287 G4TwoPeaksXS* x = (*ptr)[i];
279 if(nullptr == x) { 288 if(nullptr == x) {
280 x = new G4TwoPeaksXS(); 289 x = new G4TwoPeaksXS();
281 (*ptr)[i] = x; 290 (*ptr)[i] = x;
282 } 291 }
283 x->e1peak = e1peak; 292 x->e1peak = e1peak;
284 x->e1deep = e1deep; 293 x->e1deep = e1deep;
285 x->e2peak = e2peak; 294 x->e2peak = e2peak;
286 x->e2deep = e2deep; 295 x->e2deep = e2deep;
287 x->e3peak = e3peak; 296 x->e3peak = e3peak;
288 } 297 }
289 // case of no 1st peak in all vectors 298 // case of no 1st peak in all vectors
290 if(!isDeep) { 299 if(!isDeep) {
291 delete ptr; 300 delete ptr;
292 ptr = nullptr; 301 ptr = nullptr;
293 return ptr; 302 return ptr;
294 } 303 }
295 // check base particles 304 // check base particles
296 if(!bld->GetBaseMaterialFlag()) { return ptr 305 if(!bld->GetBaseMaterialFlag()) { return ptr; }
297 306
298 auto theDensityIdx = bld->GetCoupleIndexes() 307 auto theDensityIdx = bld->GetCoupleIndexes();
299 // second loop using base materials 308 // second loop using base materials
300 for (G4int i=0; i<n; ++i) { 309 for (G4int i=0; i<n; ++i) {
301 const G4PhysicsVector* pv = (*p)[i]; 310 const G4PhysicsVector* pv = (*p)[i];
302 if (nullptr == pv) { 311 if (nullptr == pv) {
303 G4int j = (*theDensityIdx)[i]; 312 G4int j = (*theDensityIdx)[i];
304 if(j == i) { continue; } 313 if(j == i) { continue; }
305 G4TwoPeaksXS* x = (*ptr)[i]; 314 G4TwoPeaksXS* x = (*ptr)[i];
306 G4TwoPeaksXS* y = (*ptr)[j]; 315 G4TwoPeaksXS* y = (*ptr)[j];
307 if(nullptr == x) { 316 if(nullptr == x) {
308 x = new G4TwoPeaksXS(); 317 x = new G4TwoPeaksXS();
309 (*ptr)[i] = x; 318 (*ptr)[i] = x;
310 } 319 }
311 x->e1peak = y->e1peak; 320 x->e1peak = y->e1peak;
312 x->e1deep = y->e1deep; 321 x->e1deep = y->e1deep;
313 x->e2peak = y->e2peak; 322 x->e2peak = y->e2peak;
314 x->e2deep = y->e2deep; 323 x->e2deep = y->e2deep;
315 x->e3peak = y->e3peak; 324 x->e3peak = y->e3peak;
316 } 325 }
317 } 326 }
318 return ptr; 327 return ptr;
319 } 328 }
320 329
321 //....oooOO0OOooo........oooOO0OOooo........oo 330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
322 331
323 void G4EmUtility::InitialiseElementSelectors(G 332 void G4EmUtility::InitialiseElementSelectors(G4VEmModel* mod,
324 const G4ParticleDefinition* par 333 const G4ParticleDefinition* part,
325 const G4DataVector& cuts, 334 const G4DataVector& cuts,
326 c 335 const G4double elow,
327 c 336 const G4double ehigh)
328 { 337 {
329 // using spline for element selectors should 338 // using spline for element selectors should be investigated in details
330 // because small number of points may provid 339 // because small number of points may provide biased results
331 // large number of points requires significa 340 // large number of points requires significant increase of memory
332 G4bool spline = false; 341 G4bool spline = false;
333 342
334 G4int nbinsPerDec = G4EmParameters::Instance 343 G4int nbinsPerDec = G4EmParameters::Instance()->NumberOfBinsPerDecade();
335 344
336 G4ProductionCutsTable* theCoupleTable= 345 G4ProductionCutsTable* theCoupleTable=
337 G4ProductionCutsTable::GetProductionCutsTa 346 G4ProductionCutsTable::GetProductionCutsTable();
338 std::size_t numOfCouples = theCoupleTable->G 347 std::size_t numOfCouples = theCoupleTable->GetTableSize();
339 348
340 // prepare vector 349 // prepare vector
341 auto elmSelectors = mod->GetElementSelectors 350 auto elmSelectors = mod->GetElementSelectors();
342 if(nullptr == elmSelectors) { 351 if(nullptr == elmSelectors) {
343 elmSelectors = new std::vector<G4EmElement 352 elmSelectors = new std::vector<G4EmElementSelector*>;
344 } 353 }
345 std::size_t nSelectors = elmSelectors->size( 354 std::size_t nSelectors = elmSelectors->size();
346 if(numOfCouples > nSelectors) { 355 if(numOfCouples > nSelectors) {
347 for(std::size_t i=nSelectors; i<numOfCoupl 356 for(std::size_t i=nSelectors; i<numOfCouples; ++i) {
348 elmSelectors->push_back(nullptr); 357 elmSelectors->push_back(nullptr);
349 } 358 }
350 nSelectors = numOfCouples; 359 nSelectors = numOfCouples;
351 } 360 }
352 361
353 // initialise vector 362 // initialise vector
354 for(std::size_t i=0; i<numOfCouples; ++i) { 363 for(std::size_t i=0; i<numOfCouples; ++i) {
355 364
356 // no need in element selectors for infini 365 // no need in element selectors for infinite cuts
357 if(cuts[i] == DBL_MAX) { continue; } 366 if(cuts[i] == DBL_MAX) { continue; }
358 367
359 auto couple = theCoupleTable->GetMaterialC 368 auto couple = theCoupleTable->GetMaterialCutsCouple((G4int)i);
360 auto mat = couple->GetMaterial(); 369 auto mat = couple->GetMaterial();
361 mod->SetCurrentCouple(couple); 370 mod->SetCurrentCouple(couple);
362 371
363 // selector already exist then delete 372 // selector already exist then delete
364 delete (*elmSelectors)[i]; 373 delete (*elmSelectors)[i];
365 374
366 G4double emin = std::max(elow, mod->MinPri 375 G4double emin = std::max(elow, mod->MinPrimaryEnergy(mat, part, cuts[i]));
367 G4double emax = std::max(ehigh, 10*emin); 376 G4double emax = std::max(ehigh, 10*emin);
368 static const G4double invlog106 = 1.0/(6*G 377 static const G4double invlog106 = 1.0/(6*G4Log(10.));
369 G4int nbins = G4lrint(nbinsPerDec*G4Log(em 378 G4int nbins = G4lrint(nbinsPerDec*G4Log(emax/emin)*invlog106);
370 nbins = std::max(nbins, 3); 379 nbins = std::max(nbins, 3);
371 380
372 (*elmSelectors)[i] = new G4EmElementSelect 381 (*elmSelectors)[i] = new G4EmElementSelector(mod,mat,nbins,
373 emin,emax,spline); 382 emin,emax,spline);
374 ((*elmSelectors)[i])->Initialise(part, cut 383 ((*elmSelectors)[i])->Initialise(part, cuts[i]);
375 /* 384 /*
376 G4cout << "G4VEmModel::InitialiseElmSele 385 G4cout << "G4VEmModel::InitialiseElmSelectors i= " << i
377 << " " << part->GetParticleName 386 << " " << part->GetParticleName()
378 << " for " << mod->GetName() << " 387 << " for " << mod->GetName() << " cut= " << cuts[i]
379 << " " << (*elmSelectors)[i] << 388 << " " << (*elmSelectors)[i] << G4endl;
380 ((*elmSelectors)[i])->Dump(part); 389 ((*elmSelectors)[i])->Dump(part);
381 */ 390 */
382 } 391 }
383 mod->SetElementSelectors(elmSelectors); 392 mod->SetElementSelectors(elmSelectors);
384 } 393 }
385 394
386 //....oooOO0OOooo........oooOO0OOooo........oo 395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
387 396