| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4VRangeToEnergyConverter class implementat 26 // G4VRangeToEnergyConverter class implementation
27 // 27 //
28 // Author: H.Kurashige, 05 October 2002 - Firs 28 // Author: H.Kurashige, 05 October 2002 - First implementation
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4VRangeToEnergyConverter.hh" 31 #include "G4VRangeToEnergyConverter.hh"
32 #include "G4ParticleTable.hh" 32 #include "G4ParticleTable.hh"
33 #include "G4Element.hh" 33 #include "G4Element.hh"
34 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh"
35 #include "G4Log.hh" 35 #include "G4Log.hh"
36 #include "G4Exp.hh" 36 #include "G4Exp.hh"
37 #include "G4AutoLock.hh" 37 #include "G4AutoLock.hh"
38 38
39 namespace 39 namespace
40 { 40 {
41 G4Mutex theREMutex = G4MUTEX_INITIALIZER; 41 G4Mutex theREMutex = G4MUTEX_INITIALIZER;
42 } 42 }
43 43
44 G4double G4VRangeToEnergyConverter::sEmin = CL 44 G4double G4VRangeToEnergyConverter::sEmin = CLHEP::keV;
45 G4double G4VRangeToEnergyConverter::sEmax = 10 45 G4double G4VRangeToEnergyConverter::sEmax = 10.*CLHEP::GeV;
46 46
47 std::vector<G4double>* G4VRangeToEnergyConvert 47 std::vector<G4double>* G4VRangeToEnergyConverter::sEnergy = nullptr;
48 48
49 G4int G4VRangeToEnergyConverter::sNbinPerDecad 49 G4int G4VRangeToEnergyConverter::sNbinPerDecade = 50;
50 G4int G4VRangeToEnergyConverter::sNbin = 350; 50 G4int G4VRangeToEnergyConverter::sNbin = 350;
51 51
52 // ------------------------------------------- 52 // --------------------------------------------------------------------
53 G4VRangeToEnergyConverter::G4VRangeToEnergyCon 53 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter()
54 { 54 {
55 if(nullptr == sEnergy) 55 if(nullptr == sEnergy)
56 { 56 {
57 G4AutoLock l(&theREMutex); 57 G4AutoLock l(&theREMutex);
58 if(nullptr == sEnergy) 58 if(nullptr == sEnergy)
59 { 59 {
60 isFirstInstance = true; 60 isFirstInstance = true;
61 } 61 }
62 l.unlock(); 62 l.unlock();
63 } 63 }
64 // this method defines lock itself 64 // this method defines lock itself
65 if(isFirstInstance) 65 if(isFirstInstance)
66 { 66 {
67 FillEnergyVector(CLHEP::keV, 10.0*CLHEP::G 67 FillEnergyVector(CLHEP::keV, 10.0*CLHEP::GeV);
68 } 68 }
69 } 69 }
70 70
71 // ------------------------------------------- 71 // --------------------------------------------------------------------
72 G4VRangeToEnergyConverter::~G4VRangeToEnergyCo 72 G4VRangeToEnergyConverter::~G4VRangeToEnergyConverter()
73 { 73 {
74 if(isFirstInstance) 74 if(isFirstInstance)
75 { 75 {
76 delete sEnergy; 76 delete sEnergy;
77 sEnergy = nullptr; 77 sEnergy = nullptr;
78 sEmin = CLHEP::keV; 78 sEmin = CLHEP::keV;
79 sEmax = 10.*CLHEP::GeV; 79 sEmax = 10.*CLHEP::GeV;
80 } 80 }
81 } 81 }
82 82
83 // ------------------------------------------- 83 // --------------------------------------------------------------------
84 G4double G4VRangeToEnergyConverter::Convert(co 84 G4double G4VRangeToEnergyConverter::Convert(const G4double rangeCut,
85 co 85 const G4Material* material)
86 { 86 {
87 #ifdef G4VERBOSE 87 #ifdef G4VERBOSE
88 if (GetVerboseLevel()>3) 88 if (GetVerboseLevel()>3)
89 { 89 {
90 G4cout << "G4VRangeToEnergyConverter::Conv 90 G4cout << "G4VRangeToEnergyConverter::Convert() - ";
91 G4cout << "Convert for " << material->GetN 91 G4cout << "Convert for " << material->GetName()
92 << " with Range Cut " << rangeCut/mm << " 92 << " with Range Cut " << rangeCut/mm << "[mm]" << G4endl;
93 } 93 }
94 #endif 94 #endif
95 95
96 G4double cut = 0.0; 96 G4double cut = 0.0;
97 if(fPDG == 22) 97 if(fPDG == 22)
98 { 98 {
99 cut = ConvertForGamma(rangeCut, material); 99 cut = ConvertForGamma(rangeCut, material);
100 } 100 }
101 else 101 else
102 { 102 {
103 cut = ConvertForElectron(rangeCut, materia 103 cut = ConvertForElectron(rangeCut, material);
104 104
105 const G4double tune = 0.025*CLHEP::mm*CLHE 105 const G4double tune = 0.025*CLHEP::mm*CLHEP::g/CLHEP::cm3;
106 const G4double lowen = 30.*CLHEP::keV; 106 const G4double lowen = 30.*CLHEP::keV;
107 if(cut < lowen) 107 if(cut < lowen)
108 { 108 {
109 // corr. should be switched on smoothly 109 // corr. should be switched on smoothly
110 cut /= (1.+(1.-cut/lowen)*tune/(rangeCut 110 cut /= (1.+(1.-cut/lowen)*tune/(rangeCut*material->GetDensity()));
111 } 111 }
112 } 112 }
113 113
114 cut = std::max(sEmin, std::min(cut, sEmax)); 114 cut = std::max(sEmin, std::min(cut, sEmax));
115 return cut; 115 return cut;
116 } 116 }
117 117
118 // ------------------------------------------- 118 // --------------------------------------------------------------------
119 void G4VRangeToEnergyConverter::SetEnergyRange 119 void G4VRangeToEnergyConverter::SetEnergyRange(const G4double lowedge,
120 120 const G4double highedge)
121 { 121 {
122 G4double ehigh = std::min(10.*CLHEP::GeV, hi 122 G4double ehigh = std::min(10.*CLHEP::GeV, highedge);
123 if(ehigh > lowedge) 123 if(ehigh > lowedge)
124 { 124 {
125 FillEnergyVector(lowedge, ehigh); 125 FillEnergyVector(lowedge, ehigh);
126 } 126 }
127 } 127 }
128 128
129 // ------------------------------------------- 129 // --------------------------------------------------------------------
130 G4double G4VRangeToEnergyConverter::GetLowEdge 130 G4double G4VRangeToEnergyConverter::GetLowEdgeEnergy()
131 { 131 {
132 return sEmin; 132 return sEmin;
133 } 133 }
134 134
135 // ------------------------------------------- 135 // --------------------------------------------------------------------
136 G4double G4VRangeToEnergyConverter::GetHighEdg 136 G4double G4VRangeToEnergyConverter::GetHighEdgeEnergy()
137 { 137 {
138 return sEmax; 138 return sEmax;
139 } 139 }
140 140
141 // ------------------------------------------- 141 // --------------------------------------------------------------------
142 142
143 G4double G4VRangeToEnergyConverter::GetMaxEner 143 G4double G4VRangeToEnergyConverter::GetMaxEnergyCut()
144 { 144 {
145 return sEmax; 145 return sEmax;
146 } 146 }
147 147
148 // ------------------------------------------- 148 // --------------------------------------------------------------------
149 void G4VRangeToEnergyConverter::SetMaxEnergyCu 149 void G4VRangeToEnergyConverter::SetMaxEnergyCut(const G4double value)
150 { 150 {
151 G4double ehigh = std::min(10.*CLHEP::GeV, va 151 G4double ehigh = std::min(10.*CLHEP::GeV, value);
152 if(ehigh > sEmin) 152 if(ehigh > sEmin)
153 { 153 {
154 FillEnergyVector(sEmin, ehigh); 154 FillEnergyVector(sEmin, ehigh);
155 } 155 }
156 } 156 }
157 157
158 // ------------------------------------------- 158 // --------------------------------------------------------------------
159 void G4VRangeToEnergyConverter::FillEnergyVect 159 void G4VRangeToEnergyConverter::FillEnergyVector(const G4double emin,
160 160 const G4double emax)
161 { 161 {
162 if(emin != sEmin || emax != sEmax || nullptr 162 if(emin != sEmin || emax != sEmax || nullptr == sEnergy)
163 { 163 {
164 sEmin = emin; 164 sEmin = emin;
165 sEmax = emax; 165 sEmax = emax;
166 sNbin = sNbinPerDecade*G4lrint(std::log10( 166 sNbin = sNbinPerDecade*G4lrint(std::log10(emax/emin));
167 if(nullptr == sEnergy) { sEnergy = new std 167 if(nullptr == sEnergy) { sEnergy = new std::vector<G4double>; }
168 sEnergy->resize(sNbin + 1); 168 sEnergy->resize(sNbin + 1);
169 (*sEnergy)[0] = emin; 169 (*sEnergy)[0] = emin;
170 (*sEnergy)[sNbin] = emax; 170 (*sEnergy)[sNbin] = emax;
171 G4double fact = G4Log(emax/emin)/sNbin; 171 G4double fact = G4Log(emax/emin)/sNbin;
172 for(G4int i=1; i<sNbin; ++i) { (*sEnergy)[ 172 for(G4int i=1; i<sNbin; ++i) { (*sEnergy)[i] = emin*G4Exp(i * fact); }
173 } 173 }
174 } 174 }
175 175
176 // ------------------------------------------- 176 // --------------------------------------------------------------------
177 G4double 177 G4double
178 G4VRangeToEnergyConverter::ConvertForGamma(con 178 G4VRangeToEnergyConverter::ConvertForGamma(const G4double rangeCut,
179 con 179 const G4Material* material)
180 { 180 {
181 const G4ElementVector* elm = material->GetEl 181 const G4ElementVector* elm = material->GetElementVector();
182 const G4double* dens = material->GetAtomicNu 182 const G4double* dens = material->GetAtomicNumDensityVector();
183 183
184 // fill absorption length vector 184 // fill absorption length vector
185 G4int nelm = (G4int)material->GetNumberOfEle 185 G4int nelm = (G4int)material->GetNumberOfElements();
186 G4double range1 = 0.0; 186 G4double range1 = 0.0;
187 G4double range2 = 0.0; 187 G4double range2 = 0.0;
188 G4double e1 = 0.0; 188 G4double e1 = 0.0;
189 G4double e2 = 0.0; 189 G4double e2 = 0.0;
190 for (G4int i=0; i<sNbin; ++i) 190 for (G4int i=0; i<sNbin; ++i)
191 { 191 {
192 e2 = (*sEnergy)[i]; 192 e2 = (*sEnergy)[i];
193 G4double sig = 0.; 193 G4double sig = 0.;
194 194
195 for (G4int j=0; j<nelm; ++j) 195 for (G4int j=0; j<nelm; ++j)
196 { 196 {
197 sig += dens[j]*ComputeValue((*elm)[j]->G 197 sig += dens[j]*ComputeValue((*elm)[j]->GetZasInt(), e2);
198 } 198 }
199 range2 = (sig > 0.0) ? 5./sig : DBL_MAX; 199 range2 = (sig > 0.0) ? 5./sig : DBL_MAX;
200 if(i == 0 || range2 < rangeCut) 200 if(i == 0 || range2 < rangeCut)
201 { 201 {
202 e1 = e2; 202 e1 = e2;
203 range1 = range2; 203 range1 = range2;
204 } 204 }
205 else 205 else
206 { 206 {
207 break; 207 break;
208 } 208 }
209 } 209 }
210 return LiniearInterpolation(e1, e2, range1, 210 return LiniearInterpolation(e1, e2, range1, range2, rangeCut);
211 } 211 }
212 212
213 // ------------------------------------------- 213 // --------------------------------------------------------------------
214 G4double 214 G4double
215 G4VRangeToEnergyConverter::ConvertForElectron( 215 G4VRangeToEnergyConverter::ConvertForElectron(const G4double rangeCut,
216 216 const G4Material* material)
217 { 217 {
218 const G4ElementVector* elm = material->GetEl 218 const G4ElementVector* elm = material->GetElementVector();
219 const G4double* dens = material->GetAtomicNu 219 const G4double* dens = material->GetAtomicNumDensityVector();
220 220
221 // fill absorption length vector 221 // fill absorption length vector
222 G4int nelm = (G4int)material->GetNumberOfEle 222 G4int nelm = (G4int)material->GetNumberOfElements();
223 G4double dedx1 = 0.0; 223 G4double dedx1 = 0.0;
224 G4double dedx2 = 0.0; 224 G4double dedx2 = 0.0;
225 G4double range1 = 0.0; 225 G4double range1 = 0.0;
226 G4double range2 = 0.0; 226 G4double range2 = 0.0;
227 G4double e1 = 0.0; 227 G4double e1 = 0.0;
228 G4double e2 = 0.0; 228 G4double e2 = 0.0;
229 G4double range = 0.; 229 G4double range = 0.;
230 for (G4int i=0; i<sNbin; ++i) 230 for (G4int i=0; i<sNbin; ++i)
231 { 231 {
232 e2 = (*sEnergy)[i]; 232 e2 = (*sEnergy)[i];
233 dedx2 = 0.0; 233 dedx2 = 0.0;
234 for (G4int j=0; j<nelm; ++j) 234 for (G4int j=0; j<nelm; ++j)
235 { 235 {
236 dedx2 += dens[j]*ComputeValue((*elm)[j]- 236 dedx2 += dens[j]*ComputeValue((*elm)[j]->GetZasInt(), e2);
237 } 237 }
238 range += (dedx1 + dedx2 > 0.0) ? 2*(e2 - e 238 range += (dedx1 + dedx2 > 0.0) ? 2*(e2 - e1)/(dedx1 + dedx2) : 0.0;
239 range2 = range; 239 range2 = range;
240 if(range2 < rangeCut) 240 if(range2 < rangeCut)
241 { 241 {
242 e1 = e2; 242 e1 = e2;
243 dedx1 = dedx2; 243 dedx1 = dedx2;
244 range1 = range2; 244 range1 = range2;
245 } 245 }
246 else 246 else
247 { 247 {
248 break; 248 break;
249 } 249 }
250 } 250 }
251 return LiniearInterpolation(e1, e2, range1, 251 return LiniearInterpolation(e1, e2, range1, range2, rangeCut);
252 } 252 }
253 253
254 // ------------------------------------------- 254 // --------------------------------------------------------------------
255 255