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