| Geant4 Cross Reference |
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 // G4VRangeToEnergyConverter class implementation
27 //
28 // Author: H.Kurashige, 05 October 2002 - First implementation
29 // --------------------------------------------------------------------
30
31 #include "G4VRangeToEnergyConverter.hh"
32 #include "G4ParticleTable.hh"
33 #include "G4Element.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4Log.hh"
36 #include "G4Exp.hh"
37 #include "G4AutoLock.hh"
38
39 namespace
40 {
41 G4Mutex theREMutex = G4MUTEX_INITIALIZER;
42 }
43
44 G4double G4VRangeToEnergyConverter::sEmin = CLHEP::keV;
45 G4double G4VRangeToEnergyConverter::sEmax = 10.*CLHEP::GeV;
46
47 std::vector<G4double>* G4VRangeToEnergyConverter::sEnergy = nullptr;
48
49 G4int G4VRangeToEnergyConverter::sNbinPerDecade = 50;
50 G4int G4VRangeToEnergyConverter::sNbin = 350;
51
52 // --------------------------------------------------------------------
53 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter()
54 {
55 if(nullptr == sEnergy)
56 {
57 G4AutoLock l(&theREMutex);
58 if(nullptr == sEnergy)
59 {
60 isFirstInstance = true;
61 }
62 l.unlock();
63 }
64 // this method defines lock itself
65 if(isFirstInstance)
66 {
67 FillEnergyVector(CLHEP::keV, 10.0*CLHEP::GeV);
68 }
69 }
70
71 // --------------------------------------------------------------------
72 G4VRangeToEnergyConverter::~G4VRangeToEnergyConverter()
73 {
74 if(isFirstInstance)
75 {
76 delete sEnergy;
77 sEnergy = nullptr;
78 sEmin = CLHEP::keV;
79 sEmax = 10.*CLHEP::GeV;
80 }
81 }
82
83 // --------------------------------------------------------------------
84 G4double G4VRangeToEnergyConverter::Convert(const G4double rangeCut,
85 const G4Material* material)
86 {
87 #ifdef G4VERBOSE
88 if (GetVerboseLevel()>3)
89 {
90 G4cout << "G4VRangeToEnergyConverter::Convert() - ";
91 G4cout << "Convert for " << material->GetName()
92 << " with Range Cut " << rangeCut/mm << "[mm]" << G4endl;
93 }
94 #endif
95
96 G4double cut = 0.0;
97 if(fPDG == 22)
98 {
99 cut = ConvertForGamma(rangeCut, material);
100 }
101 else
102 {
103 cut = ConvertForElectron(rangeCut, material);
104
105 const G4double tune = 0.025*CLHEP::mm*CLHEP::g/CLHEP::cm3;
106 const G4double lowen = 30.*CLHEP::keV;
107 if(cut < lowen)
108 {
109 // corr. should be switched on smoothly
110 cut /= (1.+(1.-cut/lowen)*tune/(rangeCut*material->GetDensity()));
111 }
112 }
113
114 cut = std::max(sEmin, std::min(cut, sEmax));
115 return cut;
116 }
117
118 // --------------------------------------------------------------------
119 void G4VRangeToEnergyConverter::SetEnergyRange(const G4double lowedge,
120 const G4double highedge)
121 {
122 G4double ehigh = std::min(10.*CLHEP::GeV, highedge);
123 if(ehigh > lowedge)
124 {
125 FillEnergyVector(lowedge, ehigh);
126 }
127 }
128
129 // --------------------------------------------------------------------
130 G4double G4VRangeToEnergyConverter::GetLowEdgeEnergy()
131 {
132 return sEmin;
133 }
134
135 // --------------------------------------------------------------------
136 G4double G4VRangeToEnergyConverter::GetHighEdgeEnergy()
137 {
138 return sEmax;
139 }
140
141 // --------------------------------------------------------------------
142
143 G4double G4VRangeToEnergyConverter::GetMaxEnergyCut()
144 {
145 return sEmax;
146 }
147
148 // --------------------------------------------------------------------
149 void G4VRangeToEnergyConverter::SetMaxEnergyCut(const G4double value)
150 {
151 G4double ehigh = std::min(10.*CLHEP::GeV, value);
152 if(ehigh > sEmin)
153 {
154 FillEnergyVector(sEmin, ehigh);
155 }
156 }
157
158 // --------------------------------------------------------------------
159 void G4VRangeToEnergyConverter::FillEnergyVector(const G4double emin,
160 const G4double emax)
161 {
162 if(emin != sEmin || emax != sEmax || nullptr == sEnergy)
163 {
164 sEmin = emin;
165 sEmax = emax;
166 sNbin = sNbinPerDecade*G4lrint(std::log10(emax/emin));
167 if(nullptr == sEnergy) { sEnergy = new std::vector<G4double>; }
168 sEnergy->resize(sNbin + 1);
169 (*sEnergy)[0] = emin;
170 (*sEnergy)[sNbin] = emax;
171 G4double fact = G4Log(emax/emin)/sNbin;
172 for(G4int i=1; i<sNbin; ++i) { (*sEnergy)[i] = emin*G4Exp(i * fact); }
173 }
174 }
175
176 // --------------------------------------------------------------------
177 G4double
178 G4VRangeToEnergyConverter::ConvertForGamma(const G4double rangeCut,
179 const G4Material* material)
180 {
181 const G4ElementVector* elm = material->GetElementVector();
182 const G4double* dens = material->GetAtomicNumDensityVector();
183
184 // fill absorption length vector
185 G4int nelm = (G4int)material->GetNumberOfElements();
186 G4double range1 = 0.0;
187 G4double range2 = 0.0;
188 G4double e1 = 0.0;
189 G4double e2 = 0.0;
190 for (G4int i=0; i<sNbin; ++i)
191 {
192 e2 = (*sEnergy)[i];
193 G4double sig = 0.;
194
195 for (G4int j=0; j<nelm; ++j)
196 {
197 sig += dens[j]*ComputeValue((*elm)[j]->GetZasInt(), e2);
198 }
199 range2 = (sig > 0.0) ? 5./sig : DBL_MAX;
200 if(i == 0 || range2 < rangeCut)
201 {
202 e1 = e2;
203 range1 = range2;
204 }
205 else
206 {
207 break;
208 }
209 }
210 return LiniearInterpolation(e1, e2, range1, range2, rangeCut);
211 }
212
213 // --------------------------------------------------------------------
214 G4double
215 G4VRangeToEnergyConverter::ConvertForElectron(const G4double rangeCut,
216 const G4Material* material)
217 {
218 const G4ElementVector* elm = material->GetElementVector();
219 const G4double* dens = material->GetAtomicNumDensityVector();
220
221 // fill absorption length vector
222 G4int nelm = (G4int)material->GetNumberOfElements();
223 G4double dedx1 = 0.0;
224 G4double dedx2 = 0.0;
225 G4double range1 = 0.0;
226 G4double range2 = 0.0;
227 G4double e1 = 0.0;
228 G4double e2 = 0.0;
229 G4double range = 0.;
230 for (G4int i=0; i<sNbin; ++i)
231 {
232 e2 = (*sEnergy)[i];
233 dedx2 = 0.0;
234 for (G4int j=0; j<nelm; ++j)
235 {
236 dedx2 += dens[j]*ComputeValue((*elm)[j]->GetZasInt(), e2);
237 }
238 range += (dedx1 + dedx2 > 0.0) ? 2*(e2 - e1)/(dedx1 + dedx2) : 0.0;
239 range2 = range;
240 if(range2 < rangeCut)
241 {
242 e1 = e2;
243 dedx1 = dedx2;
244 range1 = range2;
245 }
246 else
247 {
248 break;
249 }
250 }
251 return LiniearInterpolation(e1, e2, range1, range2, rangeCut);
252 }
253
254 // --------------------------------------------------------------------
255