| 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 // G4VRangeToEnergyConverter class implementat
27 //
28 // Author: H.Kurashige, 05 October 2002 - Firs
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 = CL
45 G4double G4VRangeToEnergyConverter::sEmax = 10
46
47 std::vector<G4double>* G4VRangeToEnergyConvert
48
49 G4int G4VRangeToEnergyConverter::sNbinPerDecad
50 G4int G4VRangeToEnergyConverter::sNbin = 350;
51
52 // -------------------------------------------
53 G4VRangeToEnergyConverter::G4VRangeToEnergyCon
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::G
68 }
69 }
70
71 // -------------------------------------------
72 G4VRangeToEnergyConverter::~G4VRangeToEnergyCo
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(co
85 co
86 {
87 #ifdef G4VERBOSE
88 if (GetVerboseLevel()>3)
89 {
90 G4cout << "G4VRangeToEnergyConverter::Conv
91 G4cout << "Convert for " << material->GetN
92 << " with Range Cut " << rangeCut/mm << "
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, materia
104
105 const G4double tune = 0.025*CLHEP::mm*CLHE
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
111 }
112 }
113
114 cut = std::max(sEmin, std::min(cut, sEmax));
115 return cut;
116 }
117
118 // -------------------------------------------
119 void G4VRangeToEnergyConverter::SetEnergyRange
120
121 {
122 G4double ehigh = std::min(10.*CLHEP::GeV, hi
123 if(ehigh > lowedge)
124 {
125 FillEnergyVector(lowedge, ehigh);
126 }
127 }
128
129 // -------------------------------------------
130 G4double G4VRangeToEnergyConverter::GetLowEdge
131 {
132 return sEmin;
133 }
134
135 // -------------------------------------------
136 G4double G4VRangeToEnergyConverter::GetHighEdg
137 {
138 return sEmax;
139 }
140
141 // -------------------------------------------
142
143 G4double G4VRangeToEnergyConverter::GetMaxEner
144 {
145 return sEmax;
146 }
147
148 // -------------------------------------------
149 void G4VRangeToEnergyConverter::SetMaxEnergyCu
150 {
151 G4double ehigh = std::min(10.*CLHEP::GeV, va
152 if(ehigh > sEmin)
153 {
154 FillEnergyVector(sEmin, ehigh);
155 }
156 }
157
158 // -------------------------------------------
159 void G4VRangeToEnergyConverter::FillEnergyVect
160
161 {
162 if(emin != sEmin || emax != sEmax || nullptr
163 {
164 sEmin = emin;
165 sEmax = emax;
166 sNbin = sNbinPerDecade*G4lrint(std::log10(
167 if(nullptr == sEnergy) { sEnergy = new std
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)[
173 }
174 }
175
176 // -------------------------------------------
177 G4double
178 G4VRangeToEnergyConverter::ConvertForGamma(con
179 con
180 {
181 const G4ElementVector* elm = material->GetEl
182 const G4double* dens = material->GetAtomicNu
183
184 // fill absorption length vector
185 G4int nelm = (G4int)material->GetNumberOfEle
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]->G
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,
211 }
212
213 // -------------------------------------------
214 G4double
215 G4VRangeToEnergyConverter::ConvertForElectron(
216
217 {
218 const G4ElementVector* elm = material->GetEl
219 const G4double* dens = material->GetAtomicNu
220
221 // fill absorption length vector
222 G4int nelm = (G4int)material->GetNumberOfEle
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]-
237 }
238 range += (dedx1 + dedx2 > 0.0) ? 2*(e2 - e
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,
252 }
253
254 // -------------------------------------------
255