ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EmElementSelector.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EmElementSelector.hh
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 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 //
32 // File name: G4EmElementSelector
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 29.05.2008
37 //
38 // Modifications:
39 //
40 //
41 // Class Description:
42 //
43 // Generic helper class for the random selection of an element
44 
45 // -------------------------------------------------------------------
46 //
47 
48 #ifndef G4EmElementSelector_h
49 #define G4EmElementSelector_h 1
50 
51 #include "globals.hh"
52 #include "G4ParticleDefinition.hh"
53 #include "G4Material.hh"
54 #include "G4Element.hh"
55 #include "G4ElementVector.hh"
56 #include "G4PhysicsLogVector.hh"
57 #include "Randomize.hh"
58 #include <vector>
59 
60 class G4VEmModel;
61 
63 {
64 
65 public:
66 
69  G4bool spline = true);
70 
72 
73  void Initialise(const G4ParticleDefinition*, G4double cut = 0.0);
74 
75  void Dump(const G4ParticleDefinition* p = nullptr);
76 
77  inline const G4Element* SelectRandomAtom(G4double kineticEnergy) const;
78  inline const G4Element* SelectRandomAtom(const G4double kineticEnergy,
79  const G4double logEKin) const;
80 
81  inline const G4Material* GetMaterial() const;
82 
83 private:
84 
85  // hide assignment operator
88 
92 
95 
99 
100  std::vector<G4PhysicsLogVector*> xSections;
101 
102 };
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
106 
108 {
109  const G4Element* element = (*theElementVector)[nElmMinusOne];
110  if (nElmMinusOne > 0) {
112  size_t idx(0);
113  for(G4int i=0; i<nElmMinusOne; ++i) {
114  if (x <= (xSections[i])->Value(e, idx)) {
115  element = (*theElementVector)[i];
116  break;
117  }
118  }
119  }
120  return element;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
124 inline const G4Element*
126 {
127  const G4Element* element = (*theElementVector)[nElmMinusOne];
128  if (nElmMinusOne > 0) {
129  // 1. Determine energy index (only once)
130  const size_t nBins = (xSections[0])->GetVectorLength();
131  // handle cases below/above the enrgy grid (by ekin, idx that gives a=0/1)
132  // ekin = x[0] if e<=x[0] and idx will be 0 ^ a=0 => so y=y0
133  // ekin = x[N-1] if e>=x[N-1] and idx will be N-2 ^ a=1 => so y=y_{N-1}
134  const G4double ekin = std::max((xSections[0])->Energy(0),
135  std::min((xSections[0])->Energy(nBins-1),e));
136  // compute the lower index of the bin (idx \in [0,N-2] will be guaranted)
137  const size_t idx = (xSections[0])->ComputeLogVectorBin(loge);
138  // 2. Do the linear interp.(robust for corner cases through ekin, idx and a)
139  const G4double x1 = (xSections[0])->Energy(idx);
140  const G4double x2 = (xSections[0])->Energy(idx+1);
141  // note: all corner cases of the previous methods are covered and eventually
142  // gives a=0/1 that results in y=y0\y_{N-1} if e<=x[0]/e>=x[N-1] or
143  // y=y_i/y_{i+1} if e<x[i]/e>=x[i+1] due to small numerical errors
144  const G4double a = std::max(0., std::min(1., (ekin - x1)/(x2 - x1)));
145  const G4double urnd = G4UniformRand();
146  for (G4int i = 0; i < nElmMinusOne; ++i) {
147  const G4double y1 = (*xSections[i])[idx];
148  const G4double y2 = (*xSections[i])[idx+1];
149  if (urnd <= y1 + a*(y2-y1)) {
150  element = (*theElementVector)[i];
151  break;
152  }
153  }
154  }
155  return element;
156 }
157 
158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
159 
161 {
162  return material;
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
166 
167 #endif
168