Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4EmUtility.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 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........oooOO0OOooo........oooOO0OOooo......
 43 
 44 static const G4double g4log10 = G4Log(10.);
 45 
 46 const G4Region* 
 47 G4EmUtility::FindRegion(const G4String& regionName, const G4int verbose)
 48 {
 49   const G4Region* reg = nullptr;
 50   G4RegionStore* regStore = G4RegionStore::GetInstance();
 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 to find a region <"
 56            << r << G4endl;
 57   } else if(verbose > 1) {
 58     G4cout << "### G4EmUtility finds out G4Region <" << r << ">" 
 59            << G4endl;
 60   } 
 61   return reg;  
 62 }
 63 
 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 65 
 66 const G4Element* G4EmUtility::SampleRandomElement(const G4Material* mat)
 67 {
 68   const G4Element* elm = mat->GetElement(0);
 69   std::size_t nElements = mat->GetNumberOfElements();
 70   if(1 < nElements) {
 71     G4double x = mat->GetTotNbOfElectPerVolume()*G4UniformRand();
 72     const G4double* y = mat->GetVecNbOfAtomsPerVolume();
 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........oooOO0OOooo........oooOO0OOooo.....
 83 
 84 const G4Isotope* G4EmUtility::SampleRandomIsotope(const G4Element* elm)
 85 {
 86   const std::size_t ni = elm->GetNumberOfIsotopes();
 87   const G4Isotope* iso = elm->GetIsotope(0);
 88   if(ni > 1) {
 89     const G4double* ab = elm->GetRelativeAbundanceVector();
 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........oooOO0OOooo........oooOO0OOooo.....
103 
104 std::vector<G4double>* G4EmUtility::FindCrossSectionMax(G4PhysicsTable* p)
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........oooOO0OOooo........oooOO0OOooo.....
147 
148 std::vector<G4double>* 
149 G4EmUtility::FindCrossSectionMax(G4VDiscreteProcess* p,
150                                  const G4ParticleDefinition* part)
151 {
152   std::vector<G4double>* ptr = nullptr;
153   if (nullptr == p || nullptr == part) { return ptr; }
154 
155   G4EmParameters* theParameters = G4EmParameters::Instance();
156   const G4double tmin = theParameters->MinKinEnergy();
157   const G4double tmax = theParameters->MaxKinEnergy();
158   const G4double ee = G4Log(tmax/tmin);
159   const G4double scale = theParameters->NumberOfBinsPerDecade()/g4log10;
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::GetProductionCutsTable();
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, theCoupleTable->GetMaterialCutsCouple(i));
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........oooOO0OOooo........oooOO0OOooo.....
200 
201 std::vector<G4TwoPeaksXS*>* 
202 G4EmUtility::FillPeaksStructure(G4PhysicsTable* p, G4LossTableBuilder* bld)
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, e3peak;
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 = DBL_MAX;
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........oooOO0OOooo........oooOO0OOooo.....
322 
323 void G4EmUtility::InitialiseElementSelectors(G4VEmModel* mod,
324                const G4ParticleDefinition* part,
325                const G4DataVector& cuts,
326                                              const G4double elow,
327                                              const G4double ehigh)
328 {
329   // using spline for element selectors should be investigated in details
330   // because small number of points may provide biased results
331   // large number of points requires significant increase of memory
332   G4bool spline = false;
333 
334   G4int nbinsPerDec = G4EmParameters::Instance()->NumberOfBinsPerDecade();
335 
336   G4ProductionCutsTable* theCoupleTable=
337     G4ProductionCutsTable::GetProductionCutsTable();
338   std::size_t numOfCouples = theCoupleTable->GetTableSize();
339 
340   // prepare vector
341   auto elmSelectors = mod->GetElementSelectors();
342   if(nullptr == elmSelectors) {
343     elmSelectors = new std::vector<G4EmElementSelector*>;
344   }
345   std::size_t nSelectors = elmSelectors->size();
346   if(numOfCouples > nSelectors) { 
347     for(std::size_t i=nSelectors; i<numOfCouples; ++i) { 
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 infinite cuts
357     if(cuts[i] == DBL_MAX) { continue; }
358    
359     auto couple = theCoupleTable->GetMaterialCutsCouple((G4int)i); 
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->MinPrimaryEnergy(mat, part, cuts[i]));
367     G4double emax = std::max(ehigh, 10*emin);
368     static const G4double invlog106 = 1.0/(6*G4Log(10.));
369     G4int nbins = G4lrint(nbinsPerDec*G4Log(emax/emin)*invlog106);
370     nbins = std::max(nbins, 3);
371 
372     (*elmSelectors)[i] = new G4EmElementSelector(mod,mat,nbins,
373              emin,emax,spline);
374     ((*elmSelectors)[i])->Initialise(part, cuts[i]);
375     /*      
376       G4cout << "G4VEmModel::InitialiseElmSelectors i= " << i 
377              << "  "  << part->GetParticleName() 
378              << " for " << mod->GetName() << "  cut= " << cuts[i] 
379              << "  " << (*elmSelectors)[i] << G4endl;      
380       ((*elmSelectors)[i])->Dump(part);
381     */
382   }
383   mod->SetElementSelectors(elmSelectors);
384 }
385 
386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
387