| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
>> 24 // $Id: G4EnergyLossTables.cc,v 1.17 2001/10/29 09:40:52 maire Exp $
>> 25 // GEANT4 tag $Name: geant4-04-00 $
27 // 26 //
28 // ------------------------------------------- << 27 // -------------------------------------------------------------------
29 // first version created by P.Urban , 06/04/19 28 // first version created by P.Urban , 06/04/1998
30 // modifications + "precise" functions added b 29 // modifications + "precise" functions added by L.Urban , 27/05/98
31 // modifications , TOF functions , 26/10/98, L 30 // modifications , TOF functions , 26/10/98, L.Urban
32 // cache mechanism in order to gain time, 11/0 31 // cache mechanism in order to gain time, 11/02/99, L.Urban
33 // bug fixed , 12/04/99 , L.Urban 32 // bug fixed , 12/04/99 , L.Urban
34 // 10.11.99: moved from RWT hash dictionary to 33 // 10.11.99: moved from RWT hash dictionary to STL map, G.Barrand, M.Maire
35 // 27.09.01 L.Urban , bug fixed (negative ener 34 // 27.09.01 L.Urban , bug fixed (negative energy deposit)
36 // 26.10.01 all static functions moved from .i 35 // 26.10.01 all static functions moved from .icc files (mma)
37 // 15.01.03 Add interfaces required for "cut p <<
38 // 12.03.03 Add warnings to obsolete interface <<
39 // 10.04.03 Add call to G4LossTableManager is <<
40 // <<
41 // ------------------------------------------- 36 // -------------------------------------------------------------------
42 37
43 #include "G4EnergyLossTables.hh" 38 #include "G4EnergyLossTables.hh"
44 #include "G4SystemOfUnits.hh" <<
45 #include "G4MaterialCutsCouple.hh" <<
46 #include "G4RegionStore.hh" <<
47 #include "G4LossTableManager.hh" <<
48 <<
49 39
50 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 41
52 G4EnergyLossTablesHelper *G4EnergyLossTables:: << 42 G4EnergyLossTablesHelper G4EnergyLossTables::t ;
53 G4EnergyLossTablesHelper *G4EnergyLossTables:: << 43 const G4ParticleDefinition* G4EnergyLossTables::lastParticle = 0;
54 G4ParticleDefinition* G4EnergyLossTables::last << 44 G4double G4EnergyLossTables::QQPositron = eplus*eplus ;
55 G4double G4EnergyLossTables::QQPositron = 1.0; <<
56 G4double G4EnergyLossTables::Chargesquare ; 45 G4double G4EnergyLossTables::Chargesquare ;
57 G4int G4EnergyLossTables::oldIndex = -1 ; 46 G4int G4EnergyLossTables::oldIndex = -1 ;
58 G4double G4EnergyLossTables::rmin = 0. ; 47 G4double G4EnergyLossTables::rmin = 0. ;
59 G4double G4EnergyLossTables::rmax = 0. ; 48 G4double G4EnergyLossTables::rmax = 0. ;
60 G4double G4EnergyLossTables::Thigh = 0. ; 49 G4double G4EnergyLossTables::Thigh = 0. ;
61 G4int G4EnergyLossTables::let_counter = 0; <<
62 G4int G4EnergyLossTables::let_max_num_warni <<
63 G4bool G4EnergyLossTables::first_loss = true <<
64 50
65 G4EnergyLossTables::helper_map *G4EnergyLossTa << 51 G4EnergyLossTables::helper_map G4EnergyLossTables::dict;
66 52
67 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 54
69 G4EnergyLossTablesHelper::G4EnergyLossTablesHe 55 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper(
70 const G4PhysicsTable* aDEDXTable, 56 const G4PhysicsTable* aDEDXTable,
71 const G4PhysicsTable* aRangeTable, 57 const G4PhysicsTable* aRangeTable,
72 const G4PhysicsTable* anInverseRangeTable, 58 const G4PhysicsTable* anInverseRangeTable,
73 const G4PhysicsTable* aLabTimeTable, 59 const G4PhysicsTable* aLabTimeTable,
74 const G4PhysicsTable* aProperTimeTable, 60 const G4PhysicsTable* aProperTimeTable,
75 G4double aLowestKineticEnergy, 61 G4double aLowestKineticEnergy,
76 G4double aHighestKineticEnergy, 62 G4double aHighestKineticEnergy,
77 G4double aMassRatio, 63 G4double aMassRatio,
78 G4int aNumberOfBins) 64 G4int aNumberOfBins)
79 : 65 :
80 theDEDXTable(aDEDXTable), theRangeTable(aRan 66 theDEDXTable(aDEDXTable), theRangeTable(aRangeTable),
81 theInverseRangeTable(anInverseRangeTable), 67 theInverseRangeTable(anInverseRangeTable),
82 theLabTimeTable(aLabTimeTable), 68 theLabTimeTable(aLabTimeTable),
83 theProperTimeTable(aProperTimeTable), 69 theProperTimeTable(aProperTimeTable),
84 theLowestKineticEnergy(aLowestKineticEnergy) 70 theLowestKineticEnergy(aLowestKineticEnergy),
85 theHighestKineticEnergy(aHighestKineticEnerg 71 theHighestKineticEnergy(aHighestKineticEnergy),
86 theMassRatio(aMassRatio), 72 theMassRatio(aMassRatio),
87 theNumberOfBins(aNumberOfBins) 73 theNumberOfBins(aNumberOfBins)
88 { } 74 { }
89 75
90 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 77
92 G4EnergyLossTablesHelper::G4EnergyLossTablesHe 78 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper()
93 { << 79 { }
94 theLowestKineticEnergy = 0.0; <<
95 theHighestKineticEnergy= 0.0; <<
96 theMassRatio = 0.0; <<
97 theNumberOfBins = 0; <<
98 theDEDXTable = theRangeTable = theInverseRan <<
99 = theProperTimeTable = nullptr; <<
100 } <<
101 80
102 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 82
104 void G4EnergyLossTables::Register( 83 void G4EnergyLossTables::Register(
105 const G4ParticleDefinition* p, 84 const G4ParticleDefinition* p,
106 const G4PhysicsTable* tDEDX, 85 const G4PhysicsTable* tDEDX,
107 const G4PhysicsTable* tRange, 86 const G4PhysicsTable* tRange,
108 const G4PhysicsTable* tInverseRange, 87 const G4PhysicsTable* tInverseRange,
109 const G4PhysicsTable* tLabTime, 88 const G4PhysicsTable* tLabTime,
110 const G4PhysicsTable* tProperTime, 89 const G4PhysicsTable* tProperTime,
111 G4double lowestKineticEnergy, 90 G4double lowestKineticEnergy,
112 G4double highestKineticEnergy, 91 G4double highestKineticEnergy,
113 G4double massRatio, 92 G4double massRatio,
114 G4int NumberOfBins) 93 G4int NumberOfBins)
115 { 94 {
116 if (!dict) dict = new G4EnergyLossTables::he << 95 dict[p]= G4EnergyLossTablesHelper(tDEDX, tRange,tInverseRange,
117 if (!null_loss) null_loss = new G4EnergyLoss <<
118 if (!t) t = new G4EnergyLossTablesHelper; <<
119 <<
120 (*dict)[p]= G4EnergyLossTablesHelper(tDEDX, <<
121 tLabTime,tProperTime,lowes 96 tLabTime,tProperTime,lowestKineticEnergy,
122 highestKineticEnergy, massRatio,Number 97 highestKineticEnergy, massRatio,NumberOfBins);
123 << 98
124 *t = GetTables(p) ; // important for cach << 99 t = GetTables(p) ; // important for cache !!!!!
125 lastParticle = (G4ParticleDefinition*) p ; << 100 lastParticle = p ;
126 Chargesquare = (p->GetPDGCharge())*(p->GetPD 101 Chargesquare = (p->GetPDGCharge())*(p->GetPDGCharge())/
127 QQPositron ; 102 QQPositron ;
128 if (first_loss ) { <<
129 *null_loss = G4EnergyLossTablesHelper( <<
130 nullptr, nullptr, nullptr, nu <<
131 first_loss = false; <<
132 } <<
133 } 103 }
134 104
135 //....oooOO0OOooo........oooOO0OOooo........oo 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136 106
137 const G4PhysicsTable* G4EnergyLossTables::GetD 107 const G4PhysicsTable* G4EnergyLossTables::GetDEDXTable(
138 const G4ParticleDefinition* p) 108 const G4ParticleDefinition* p)
139 { 109 {
140 if (!dict) dict = new G4EnergyLossTables::he <<
141 helper_map::iterator it; 110 helper_map::iterator it;
142 if((it=dict->find(p))==dict->end()) return n << 111 if((it=dict.find(p))==dict.end()) return 0;
143 return (*it).second.theDEDXTable; << 112 return (*it).second.theDEDXTable;
144 } 113 }
145 114
146 //....oooOO0OOooo........oooOO0OOooo........oo 115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 116
148 const G4PhysicsTable* G4EnergyLossTables::GetR 117 const G4PhysicsTable* G4EnergyLossTables::GetRangeTable(
149 const G4ParticleDefinition* p) 118 const G4ParticleDefinition* p)
150 { 119 {
151 if (!dict) dict = new G4EnergyLossTables::he <<
152 helper_map::iterator it; 120 helper_map::iterator it;
153 if((it=dict->find(p))==dict->end()) return n << 121 if((it=dict.find(p))==dict.end()) return 0;
154 return (*it).second.theRangeTable; << 122 return (*it).second.theRangeTable;
155 } 123 }
156 124
157 //....oooOO0OOooo........oooOO0OOooo........oo 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158 126
159 const G4PhysicsTable* G4EnergyLossTables::GetI 127 const G4PhysicsTable* G4EnergyLossTables::GetInverseRangeTable(
160 const G4ParticleDefinition* p) 128 const G4ParticleDefinition* p)
161 { 129 {
162 if (!dict) dict = new G4EnergyLossTables::he <<
163 helper_map::iterator it; 130 helper_map::iterator it;
164 if((it=dict->find(p))==dict->end()) return n << 131 if((it=dict.find(p))==dict.end()) return 0;
165 return (*it).second.theInverseRangeTable; << 132 return (*it).second.theInverseRangeTable;
166 } 133 }
167 134
168 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
169 136
170 const G4PhysicsTable* G4EnergyLossTables::GetL 137 const G4PhysicsTable* G4EnergyLossTables::GetLabTimeTable(
171 const G4ParticleDefinition* p) 138 const G4ParticleDefinition* p)
172 { 139 {
173 if (!dict) dict = new G4EnergyLossTables::he <<
174 helper_map::iterator it; 140 helper_map::iterator it;
175 if((it=dict->find(p))==dict->end()) return n << 141 if((it=dict.find(p))==dict.end()) return 0;
176 return (*it).second.theLabTimeTable; << 142 return (*it).second.theLabTimeTable;
177 } 143 }
178 144
179 //....oooOO0OOooo........oooOO0OOooo........oo 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180 146
181 const G4PhysicsTable* G4EnergyLossTables::GetP 147 const G4PhysicsTable* G4EnergyLossTables::GetProperTimeTable(
182 const G4ParticleDefinition* p) 148 const G4ParticleDefinition* p)
183 { 149 {
184 if (!dict) dict = new G4EnergyLossTables::he <<
185 helper_map::iterator it; 150 helper_map::iterator it;
186 if((it=dict->find(p))==dict->end()) return n << 151 if((it=dict.find(p))==dict.end()) return 0;
187 return (*it).second.theProperTimeTable; << 152 return (*it).second.theProperTimeTable;
188 } 153 }
189 154
190 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 156
192 G4EnergyLossTablesHelper G4EnergyLossTables::G 157 G4EnergyLossTablesHelper G4EnergyLossTables::GetTables(
193 const G4ParticleDefinition* p) 158 const G4ParticleDefinition* p)
194 { 159 {
195 if (!dict) dict = new G4EnergyLossTables::he <<
196 if (!null_loss) null_loss = new G4EnergyLoss <<
197 <<
198 helper_map::iterator it; 160 helper_map::iterator it;
199 if ((it=dict->find(p))==dict->end()) { << 161 if((it=dict.find(p))==dict.end()) {
200 return *null_loss; << 162 G4Exception("G4EnergyLossTables::GetTables: table not found!");
>> 163 exit(1);
201 } 164 }
202 return (*it).second; 165 return (*it).second;
203 } 166 }
204 167
205 //....oooOO0OOooo........oooOO0OOooo........oo 168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
206 169
207 G4double G4EnergyLossTables::GetDEDX( 170 G4double G4EnergyLossTables::GetDEDX(
208 const G4ParticleDefinition *aParticle, 171 const G4ParticleDefinition *aParticle,
209 G4double KineticEnergy, 172 G4double KineticEnergy,
210 const G4Material *aMaterial) 173 const G4Material *aMaterial)
211 { 174 {
212 if (!t) t = new G4EnergyLossTablesHelper; << 175 if(aParticle != lastParticle)
213 <<
214 CPRWarning(); <<
215 if(aParticle != (const G4ParticleDefinition* <<
216 { 176 {
217 *t= GetTables(aParticle); << 177 t= GetTables(aParticle);
218 lastParticle = (G4ParticleDefinition*) aPa << 178 lastParticle = aParticle ;
219 Chargesquare = (aParticle->GetPDGCharge()) 179 Chargesquare = (aParticle->GetPDGCharge())*
220 (aParticle->GetPDGCharge()) 180 (aParticle->GetPDGCharge())/
221 QQPositron ; 181 QQPositron ;
222 oldIndex = -1 ; 182 oldIndex = -1 ;
223 } 183 }
224 const G4PhysicsTable* dEdxTable= t->theDEDX << 184 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
225 if (!dEdxTable) { <<
226 ParticleHaveNoLoss(aParticle,"dEdx"); <<
227 return 0.0; <<
228 } <<
229 185
230 G4int materialIndex = (G4int)aMaterial->GetI << 186 G4int materialIndex = aMaterial->GetIndex();
231 G4double scaledKineticEnergy = KineticEnergy << 187 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
232 G4double dEdx; 188 G4double dEdx;
233 G4bool isOut; 189 G4bool isOut;
234 190
235 if (scaledKineticEnergy<t->theLowestKineticE << 191 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
236 192
237 dEdx =(*dEdxTable)(materialIndex)->GetVal 193 dEdx =(*dEdxTable)(materialIndex)->GetValue(
238 t->theLowestKineticEnergy,isOut) << 194 t.theLowestKineticEnergy,isOut)
239 *std::sqrt(scaledKineticEnergy/t->t << 195 *sqrt(scaledKineticEnergy/t.theLowestKineticEnergy);
240 196
241 } else if (scaledKineticEnergy>t->theHighest << 197 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
242 198
243 dEdx = (*dEdxTable)(materialIndex)->GetVa 199 dEdx = (*dEdxTable)(materialIndex)->GetValue(
244 t->theHighestKineticEnergy,isOut); << 200 t.theHighestKineticEnergy,isOut);
245 201
246 } else { 202 } else {
247 << 203
248 dEdx = (*dEdxTable)(materialIndex)->GetVal 204 dEdx = (*dEdxTable)(materialIndex)->GetValue(
249 scaledKineticEnergy,isOut); 205 scaledKineticEnergy,isOut);
250 206
251 } 207 }
252 208
253 return dEdx*Chargesquare; 209 return dEdx*Chargesquare;
254 } 210 }
255 211
256 //....oooOO0OOooo........oooOO0OOooo........oo 212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
257 213
258 G4double G4EnergyLossTables::GetLabTime( 214 G4double G4EnergyLossTables::GetLabTime(
259 const G4ParticleDefinition *aParticle, 215 const G4ParticleDefinition *aParticle,
260 G4double KineticEnergy, 216 G4double KineticEnergy,
261 const G4Material *aMaterial) 217 const G4Material *aMaterial)
262 { 218 {
263 if (!t) t = new G4EnergyLossTablesHelper; << 219 if(aParticle != lastParticle)
264 <<
265 CPRWarning(); <<
266 if(aParticle != (const G4ParticleDefinition* <<
267 { 220 {
268 *t= GetTables(aParticle); << 221 t= GetTables(aParticle);
269 lastParticle = (G4ParticleDefinition*) aPa << 222 lastParticle = aParticle ;
270 oldIndex = -1 ; 223 oldIndex = -1 ;
271 } 224 }
272 const G4PhysicsTable* labtimeTable= t->theLa << 225 const G4PhysicsTable* labtimeTable= t.theLabTimeTable;
273 if (!labtimeTable) { <<
274 ParticleHaveNoLoss(aParticle,"LabTime"); <<
275 return 0.0; <<
276 } <<
277 226
278 const G4double parlowen=0.4 , ppar=0.5-parlo << 227 const G4double parlowen=0.4 , ppar=0.5-parlowen ;
279 G4int materialIndex = (G4int)aMaterial->GetI << 228 G4int materialIndex = aMaterial->GetIndex();
280 G4double scaledKineticEnergy = KineticEnergy << 229 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
281 G4double time; 230 G4double time;
282 G4bool isOut; 231 G4bool isOut;
283 232
284 if (scaledKineticEnergy<t->theLowestKineticE << 233 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
285 234
286 time = std::exp(ppar*std::log(scaledKinet << 235 time = exp(ppar*log(scaledKineticEnergy/t.theLowestKineticEnergy))*
287 (*labtimeTable)(materialIndex)->Ge 236 (*labtimeTable)(materialIndex)->GetValue(
288 t->theLowestKineticEnergy,isOut) << 237 t.theLowestKineticEnergy,isOut);
289 238
290 239
291 } else if (scaledKineticEnergy>t->theHighest << 240 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
292 241
293 time = (*labtimeTable)(materialIndex)->Ge 242 time = (*labtimeTable)(materialIndex)->GetValue(
294 t->theHighestKineticEnergy,isOut << 243 t.theHighestKineticEnergy,isOut);
295 244
296 } else { 245 } else {
297 << 246
298 time = (*labtimeTable)(materialIndex)->Get 247 time = (*labtimeTable)(materialIndex)->GetValue(
299 scaledKineticEnergy,isOut); 248 scaledKineticEnergy,isOut);
300 249
301 } 250 }
302 251
303 return time/t->theMassRatio ; << 252 return time/t.theMassRatio ;
304 } 253 }
305 254
306 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
307 256
308 G4double G4EnergyLossTables::GetDeltaLabTime( 257 G4double G4EnergyLossTables::GetDeltaLabTime(
309 const G4ParticleDefinition *aParticle, 258 const G4ParticleDefinition *aParticle,
310 G4double KineticEnergyStart, 259 G4double KineticEnergyStart,
311 G4double KineticEnergyEnd, 260 G4double KineticEnergyEnd,
312 const G4Material *aMaterial) 261 const G4Material *aMaterial)
313 { 262 {
314 if (!t) t = new G4EnergyLossTablesHelper; << 263 if(aParticle != lastParticle)
315 <<
316 CPRWarning(); <<
317 if(aParticle != (const G4ParticleDefinition* <<
318 { 264 {
319 *t= GetTables(aParticle); << 265 t= GetTables(aParticle);
320 lastParticle = (G4ParticleDefinition*) aPa << 266 lastParticle = aParticle ;
321 oldIndex = -1 ; 267 oldIndex = -1 ;
322 } 268 }
323 const G4PhysicsTable* labtimeTable= t->theLa << 269 const G4PhysicsTable* labtimeTable= t.theLabTimeTable;
324 if (!labtimeTable) { <<
325 ParticleHaveNoLoss(aParticle,"LabTime"); <<
326 return 0.0; <<
327 } <<
328 270
329 const G4double parlowen=0.4 , ppar=0.5-parlo 271 const G4double parlowen=0.4 , ppar=0.5-parlowen ;
330 const G4double dToverT = 0.05 , facT = 1. -d << 272 const G4double dToverT = 0.05 , facT = 1. -dToverT ;
331 G4double timestart,timeend,deltatime,dTT; 273 G4double timestart,timeend,deltatime,dTT;
332 G4bool isOut; 274 G4bool isOut;
333 275
334 G4int materialIndex = (G4int)aMaterial->GetI << 276 G4int materialIndex = aMaterial->GetIndex();
335 G4double scaledKineticEnergy = KineticEnergy << 277 G4double scaledKineticEnergy = KineticEnergyStart*t.theMassRatio;
336 278
337 if (scaledKineticEnergy<t->theLowestKineticE << 279 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
338 280
339 timestart = std::exp(ppar*std::log(scaled << 281 timestart = exp(ppar*log(scaledKineticEnergy/t.theLowestKineticEnergy))*
340 (*labtimeTable)(materialIndex) 282 (*labtimeTable)(materialIndex)->GetValue(
341 t->theLowestKineticEnergy,isOu << 283 t.theLowestKineticEnergy,isOut);
342 284
343 285
344 } else if (scaledKineticEnergy>t->theHighest << 286 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
345 287
346 timestart = (*labtimeTable)(materialIndex 288 timestart = (*labtimeTable)(materialIndex)->GetValue(
347 t->theHighestKineticEnergy,isO << 289 t.theHighestKineticEnergy,isOut);
348 290
349 } else { 291 } else {
350 292
351 timestart = (*labtimeTable)(materialIndex) 293 timestart = (*labtimeTable)(materialIndex)->GetValue(
352 scaledKineticEnergy,isOut); 294 scaledKineticEnergy,isOut);
353 295
354 } 296 }
355 297
356 dTT = (KineticEnergyStart - KineticEnergyEnd 298 dTT = (KineticEnergyStart - KineticEnergyEnd)/KineticEnergyStart ;
357 << 299
358 if( dTT < dToverT ) 300 if( dTT < dToverT )
359 scaledKineticEnergy = facT*KineticEnergySt << 301 scaledKineticEnergy = facT*KineticEnergyStart*t.theMassRatio;
360 else 302 else
361 scaledKineticEnergy = KineticEnergyEnd*t-> << 303 scaledKineticEnergy = KineticEnergyEnd*t.theMassRatio;
>> 304
>> 305 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
362 306
363 if (scaledKineticEnergy<t->theLowestKineticE << 307 timeend = exp(ppar*log(scaledKineticEnergy/t.theLowestKineticEnergy))*
364 <<
365 timeend = std::exp(ppar*std::log(scaledKi <<
366 (*labtimeTable)(materialIndex) 308 (*labtimeTable)(materialIndex)->GetValue(
367 t->theLowestKineticEnergy,isOu << 309 t.theLowestKineticEnergy,isOut);
368 310
369 311
370 } else if (scaledKineticEnergy>t->theHighest << 312 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
371 313
372 timeend = (*labtimeTable)(materialIndex)- 314 timeend = (*labtimeTable)(materialIndex)->GetValue(
373 t->theHighestKineticEnergy,isO << 315 t.theHighestKineticEnergy,isOut);
374 316
375 } else { 317 } else {
376 318
377 timeend = (*labtimeTable)(materialIndex)-> 319 timeend = (*labtimeTable)(materialIndex)->GetValue(
378 scaledKineticEnergy,isOut); 320 scaledKineticEnergy,isOut);
379 321
380 } 322 }
381 323
382 deltatime = timestart - timeend ; 324 deltatime = timestart - timeend ;
383 325
384 if( dTT < dToverT ) 326 if( dTT < dToverT )
385 deltatime *= dTT/dToverT; 327 deltatime *= dTT/dToverT;
386 328
387 return deltatime/t->theMassRatio ; << 329 return deltatime/t.theMassRatio ;
388 } 330 }
389 331
390 //....oooOO0OOooo........oooOO0OOooo........oo 332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
391 333
392 G4double G4EnergyLossTables::GetProperTime( 334 G4double G4EnergyLossTables::GetProperTime(
393 const G4ParticleDefinition *aParticle, 335 const G4ParticleDefinition *aParticle,
394 G4double KineticEnergy, 336 G4double KineticEnergy,
395 const G4Material *aMaterial) 337 const G4Material *aMaterial)
396 { 338 {
397 if (!t) t = new G4EnergyLossTablesHelper; << 339 if(aParticle != lastParticle)
398 <<
399 CPRWarning(); <<
400 if(aParticle != (const G4ParticleDefinition* <<
401 { 340 {
402 *t= GetTables(aParticle); << 341 t= GetTables(aParticle);
403 lastParticle = (G4ParticleDefinition*) aPa << 342 lastParticle = aParticle ;
404 oldIndex = -1 ; 343 oldIndex = -1 ;
405 } 344 }
406 const G4PhysicsTable* propertimeTable= t->th << 345 const G4PhysicsTable* propertimeTable= t.theProperTimeTable;
407 if (!propertimeTable) { <<
408 ParticleHaveNoLoss(aParticle,"ProperTime") <<
409 return 0.0; <<
410 } <<
411 346
412 const G4double parlowen=0.4 , ppar=0.5-parlo 347 const G4double parlowen=0.4 , ppar=0.5-parlowen ;
413 G4int materialIndex = (G4int)aMaterial->GetI << 348 G4int materialIndex = aMaterial->GetIndex();
414 G4double scaledKineticEnergy = KineticEnergy << 349 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
415 G4double time; 350 G4double time;
416 G4bool isOut; 351 G4bool isOut;
417 352
418 if (scaledKineticEnergy<t->theLowestKineticE << 353 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
419 354
420 time = std::exp(ppar*std::log(scaledKinet << 355 time = exp(ppar*log(scaledKineticEnergy/t.theLowestKineticEnergy))*
421 (*propertimeTable)(materialIndex)- 356 (*propertimeTable)(materialIndex)->GetValue(
422 t->theLowestKineticEnergy,isOut) << 357 t.theLowestKineticEnergy,isOut);
423 358
424 359
425 } else if (scaledKineticEnergy>t->theHighest << 360 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
426 361
427 time = (*propertimeTable)(materialIndex)- 362 time = (*propertimeTable)(materialIndex)->GetValue(
428 t->theHighestKineticEnergy,isOut << 363 t.theHighestKineticEnergy,isOut);
429 364
430 } else { 365 } else {
431 << 366
432 time = (*propertimeTable)(materialIndex)-> 367 time = (*propertimeTable)(materialIndex)->GetValue(
433 scaledKineticEnergy,isOut); 368 scaledKineticEnergy,isOut);
434 369
435 } 370 }
436 371
437 return time/t->theMassRatio ; << 372 return time/t.theMassRatio ;
438 } 373 }
439 374
440 //....oooOO0OOooo........oooOO0OOooo........oo 375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
441 376
442 G4double G4EnergyLossTables::GetDeltaProperTim 377 G4double G4EnergyLossTables::GetDeltaProperTime(
443 const G4ParticleDefinition *aParticle, 378 const G4ParticleDefinition *aParticle,
444 G4double KineticEnergyStart, 379 G4double KineticEnergyStart,
445 G4double KineticEnergyEnd, 380 G4double KineticEnergyEnd,
446 const G4Material *aMaterial) 381 const G4Material *aMaterial)
447 { 382 {
448 if (!t) t = new G4EnergyLossTablesHelper; << 383 if(aParticle != lastParticle)
449 <<
450 CPRWarning(); <<
451 if(aParticle != (const G4ParticleDefinition* <<
452 { 384 {
453 *t= GetTables(aParticle); << 385 t= GetTables(aParticle);
454 lastParticle = (G4ParticleDefinition*) aPa << 386 lastParticle = aParticle ;
455 oldIndex = -1 ; 387 oldIndex = -1 ;
456 } 388 }
457 const G4PhysicsTable* propertimeTable= t->th << 389 const G4PhysicsTable* propertimeTable= t.theProperTimeTable;
458 if (!propertimeTable) { <<
459 ParticleHaveNoLoss(aParticle,"ProperTime") <<
460 return 0.0; <<
461 } <<
462 390
463 const G4double parlowen=0.4 , ppar=0.5-parlo 391 const G4double parlowen=0.4 , ppar=0.5-parlowen ;
464 const G4double dToverT = 0.05 , facT = 1. -d << 392 const G4double dToverT = 0.05 , facT = 1. -dToverT ;
465 G4double timestart,timeend,deltatime,dTT; 393 G4double timestart,timeend,deltatime,dTT;
466 G4bool isOut; 394 G4bool isOut;
467 395
468 G4int materialIndex = (G4int)aMaterial->GetI << 396 G4int materialIndex = aMaterial->GetIndex();
469 G4double scaledKineticEnergy = KineticEnergy << 397 G4double scaledKineticEnergy = KineticEnergyStart*t.theMassRatio;
470 398
471 if (scaledKineticEnergy<t->theLowestKineticE << 399 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
472 400
473 timestart = std::exp(ppar*std::log(scaled << 401 timestart = exp(ppar*log(scaledKineticEnergy/t.theLowestKineticEnergy))*
474 (*propertimeTable)(materialInd 402 (*propertimeTable)(materialIndex)->GetValue(
475 t->theLowestKineticEnergy,isOu << 403 t.theLowestKineticEnergy,isOut);
476 404
477 405
478 } else if (scaledKineticEnergy>t->theHighest << 406 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
479 407
480 timestart = (*propertimeTable)(materialIn 408 timestart = (*propertimeTable)(materialIndex)->GetValue(
481 t->theHighestKineticEnergy,isO << 409 t.theHighestKineticEnergy,isOut);
482 410
483 } else { 411 } else {
484 412
485 timestart = (*propertimeTable)(materialInd 413 timestart = (*propertimeTable)(materialIndex)->GetValue(
486 scaledKineticEnergy,isOut); 414 scaledKineticEnergy,isOut);
487 415
488 } 416 }
489 417
490 dTT = (KineticEnergyStart - KineticEnergyEnd 418 dTT = (KineticEnergyStart - KineticEnergyEnd)/KineticEnergyStart ;
491 << 419
492 if( dTT < dToverT ) 420 if( dTT < dToverT )
493 scaledKineticEnergy = facT*KineticEnergySt << 421 scaledKineticEnergy = facT*KineticEnergyStart*t.theMassRatio;
494 else 422 else
495 scaledKineticEnergy = KineticEnergyEnd*t-> << 423 scaledKineticEnergy = KineticEnergyEnd*t.theMassRatio;
>> 424
>> 425 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
496 426
497 if (scaledKineticEnergy<t->theLowestKineticE << 427 timeend = exp(ppar*log(scaledKineticEnergy/t.theLowestKineticEnergy))*
498 <<
499 timeend = std::exp(ppar*std::log(scaledKi <<
500 (*propertimeTable)(materialInd 428 (*propertimeTable)(materialIndex)->GetValue(
501 t->theLowestKineticEnergy,isOu << 429 t.theLowestKineticEnergy,isOut);
502 430
503 431
504 } else if (scaledKineticEnergy>t->theHighest << 432 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
505 433
506 timeend = (*propertimeTable)(materialInde 434 timeend = (*propertimeTable)(materialIndex)->GetValue(
507 t->theHighestKineticEnergy,isO << 435 t.theHighestKineticEnergy,isOut);
508 436
509 } else { 437 } else {
510 438
511 timeend = (*propertimeTable)(materialIndex 439 timeend = (*propertimeTable)(materialIndex)->GetValue(
512 scaledKineticEnergy,isOut); 440 scaledKineticEnergy,isOut);
513 441
514 } 442 }
515 443
516 deltatime = timestart - timeend ; 444 deltatime = timestart - timeend ;
517 445
518 if( dTT < dToverT ) 446 if( dTT < dToverT )
519 deltatime *= dTT/dToverT ; 447 deltatime *= dTT/dToverT ;
520 448
521 return deltatime/t->theMassRatio ; << 449 return deltatime/t.theMassRatio ;
522 } 450 }
523 451
524 //....oooOO0OOooo........oooOO0OOooo........oo 452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
525 453
526 G4double G4EnergyLossTables::GetRange( 454 G4double G4EnergyLossTables::GetRange(
527 const G4ParticleDefinition *aParticle, 455 const G4ParticleDefinition *aParticle,
528 G4double KineticEnergy, 456 G4double KineticEnergy,
529 const G4Material *aMaterial) 457 const G4Material *aMaterial)
530 { 458 {
531 if (!t) t = new G4EnergyLossTablesHelper; << 459 if(aParticle != lastParticle)
532 <<
533 CPRWarning(); <<
534 if(aParticle != (const G4ParticleDefinition* <<
535 { 460 {
536 *t= GetTables(aParticle); << 461 t= GetTables(aParticle);
537 lastParticle = (G4ParticleDefinition*) aPa << 462 lastParticle = aParticle ;
538 Chargesquare = (aParticle->GetPDGCharge()) 463 Chargesquare = (aParticle->GetPDGCharge())*
539 (aParticle->GetPDGCharge()) 464 (aParticle->GetPDGCharge())/
540 QQPositron ; 465 QQPositron ;
541 oldIndex = -1 ; 466 oldIndex = -1 ;
542 } 467 }
543 const G4PhysicsTable* rangeTable= t->theRang << 468 const G4PhysicsTable* rangeTable= t.theRangeTable;
544 const G4PhysicsTable* dEdxTable= t->theDEDX << 469 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
545 if (!rangeTable) { <<
546 ParticleHaveNoLoss(aParticle,"Range"); <<
547 return 0.0; <<
548 } <<
549 470
550 G4int materialIndex = (G4int)aMaterial->GetI << 471 G4int materialIndex = aMaterial->GetIndex();
551 G4double scaledKineticEnergy = KineticEnergy << 472 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
552 G4double Range; 473 G4double Range;
553 G4bool isOut; 474 G4bool isOut;
554 475
555 if (scaledKineticEnergy<t->theLowestKineticE << 476 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
556 477
557 Range = std::sqrt(scaledKineticEnergy/t->t << 478 Range = sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)*
558 (*rangeTable)(materialIndex)->GetV 479 (*rangeTable)(materialIndex)->GetValue(
559 t->theLowestKineticEnergy,isOut) << 480 t.theLowestKineticEnergy,isOut);
560 481
561 } else if (scaledKineticEnergy>t->theHighest << 482 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
562 483
563 Range = (*rangeTable)(materialIndex)->GetV 484 Range = (*rangeTable)(materialIndex)->GetValue(
564 t->theHighestKineticEnergy,isOut)+ << 485 t.theHighestKineticEnergy,isOut)+
565 (scaledKineticEnergy-t->theHighest << 486 (scaledKineticEnergy-t.theHighestKineticEnergy)/
566 (*dEdxTable)(materialIndex)->GetVa 487 (*dEdxTable)(materialIndex)->GetValue(
567 t->theHighestKineticEnergy,isOut << 488 t.theHighestKineticEnergy,isOut);
568 489
569 } else { 490 } else {
570 << 491
571 Range = (*rangeTable)(materialIndex)->GetV 492 Range = (*rangeTable)(materialIndex)->GetValue(
572 scaledKineticEnergy,isOut); 493 scaledKineticEnergy,isOut);
573 494
574 } 495 }
575 496
576 return Range/(Chargesquare*t->theMassRatio); << 497 return Range/(Chargesquare*t.theMassRatio);
577 } 498 }
578 499
579 //....oooOO0OOooo........oooOO0OOooo........oo 500 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
580 501
581 G4double G4EnergyLossTables::GetPreciseEnergyF 502 G4double G4EnergyLossTables::GetPreciseEnergyFromRange(
582 const G4P 503 const G4ParticleDefinition *aParticle,
583 G4d 504 G4double range,
584 const G4M 505 const G4Material *aMaterial)
585 // it returns the value of the kinetic energy 506 // it returns the value of the kinetic energy for a given range
586 { 507 {
587 if (!t) t = new G4EnergyLossTablesHelper; << 508 if( aParticle != lastParticle)
588 <<
589 CPRWarning(); <<
590 if( aParticle != (const G4ParticleDefinition <<
591 { 509 {
592 *t= GetTables(aParticle); << 510 t= GetTables(aParticle);
593 lastParticle = (G4ParticleDefinition*) aPa << 511 lastParticle = aParticle;
594 Chargesquare = (aParticle->GetPDGCharge()) 512 Chargesquare = (aParticle->GetPDGCharge())*
595 (aParticle->GetPDGCharge()) 513 (aParticle->GetPDGCharge())/
596 QQPositron ; 514 QQPositron ;
597 oldIndex = -1 ; 515 oldIndex = -1 ;
598 } 516 }
599 const G4PhysicsTable* dEdxTable= t->theDEDX << 517 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
600 const G4PhysicsTable* inverseRangeTable= t- << 518 const G4PhysicsTable* inverseRangeTable= t.theInverseRangeTable;
601 if (!inverseRangeTable) { <<
602 ParticleHaveNoLoss(aParticle,"InverseRange <<
603 return 0.0; <<
604 } <<
605 519
606 G4double scaledrange,scaledKineticEnergy ; 520 G4double scaledrange,scaledKineticEnergy ;
607 G4bool isOut ; 521 G4bool isOut ;
608 522
609 G4int materialIndex = (G4int)aMaterial->GetI << 523 G4int materialIndex = aMaterial->GetIndex() ;
610 524
611 if(materialIndex != oldIndex) 525 if(materialIndex != oldIndex)
612 { 526 {
613 oldIndex = materialIndex ; 527 oldIndex = materialIndex ;
614 rmin = (*inverseRangeTable)(materialIndex) 528 rmin = (*inverseRangeTable)(materialIndex)->
615 GetLowEdgeEnergy 529 GetLowEdgeEnergy(0) ;
616 rmax = (*inverseRangeTable)(materialIndex) 530 rmax = (*inverseRangeTable)(materialIndex)->
617 GetLowEdgeEnergy(t->theNumb << 531 GetLowEdgeEnergy(t.theNumberOfBins-2) ;
618 Thigh = (*inverseRangeTable)(materialIndex 532 Thigh = (*inverseRangeTable)(materialIndex)->
619 GetValue(rmax,is 533 GetValue(rmax,isOut) ;
620 } 534 }
621 535
622 scaledrange = range*Chargesquare*t->theMassR << 536 scaledrange = range*Chargesquare*t.theMassRatio ;
623 537
624 if(scaledrange < rmin) 538 if(scaledrange < rmin)
625 { 539 {
626 scaledKineticEnergy = t->theLowestKineticE << 540 scaledKineticEnergy = t.theLowestKineticEnergy*
627 scaledrange*scaledrange/(rm 541 scaledrange*scaledrange/(rmin*rmin) ;
628 } 542 }
629 else 543 else
630 { 544 {
631 if(scaledrange < rmax) 545 if(scaledrange < rmax)
632 { 546 {
633 scaledKineticEnergy = (*inverseRangeTabl 547 scaledKineticEnergy = (*inverseRangeTable)(materialIndex)->
634 GetValue( scaled 548 GetValue( scaledrange,isOut) ;
635 } 549 }
636 else 550 else
637 { 551 {
638 scaledKineticEnergy = Thigh + 552 scaledKineticEnergy = Thigh +
639 (scaledrange-rmax)* 553 (scaledrange-rmax)*
640 (*dEdxTable)(materialInd 554 (*dEdxTable)(materialIndex)->
641 GetValue(Thig 555 GetValue(Thigh,isOut) ;
642 } 556 }
643 } 557 }
644 558
645 return scaledKineticEnergy/t->theMassRatio ; << 559 return scaledKineticEnergy/t.theMassRatio ;
646 } 560 }
647 561
648 //....oooOO0OOooo........oooOO0OOooo........oo 562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
649 563
650 G4double G4EnergyLossTables::GetPreciseDEDX( 564 G4double G4EnergyLossTables::GetPreciseDEDX(
651 const G4ParticleDefinition *aParticle, 565 const G4ParticleDefinition *aParticle,
652 G4double KineticEnergy, 566 G4double KineticEnergy,
653 const G4Material *aMaterial) 567 const G4Material *aMaterial)
654 { 568 {
655 if (!t) t = new G4EnergyLossTablesHelper; <<
656 569
657 CPRWarning(); << 570 if( aParticle != lastParticle)
658 if( aParticle != (const G4ParticleDefinition <<
659 { 571 {
660 *t= GetTables(aParticle); << 572 t= GetTables(aParticle);
661 lastParticle = (G4ParticleDefinition*) aPa << 573 lastParticle = aParticle;
662 Chargesquare = (aParticle->GetPDGCharge()) 574 Chargesquare = (aParticle->GetPDGCharge())*
663 (aParticle->GetPDGCharge()) 575 (aParticle->GetPDGCharge())/
664 QQPositron ; 576 QQPositron ;
665 oldIndex = -1 ; 577 oldIndex = -1 ;
666 } 578 }
667 const G4PhysicsTable* dEdxTable= t->theDEDX << 579 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
668 if (!dEdxTable) { <<
669 ParticleHaveNoLoss(aParticle,"dEdx"); <<
670 return 0.0; <<
671 } <<
672 580
673 G4int materialIndex = (G4int)aMaterial->GetI << 581 G4int materialIndex = aMaterial->GetIndex();
674 G4double scaledKineticEnergy = KineticEnergy << 582 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
675 G4double dEdx; 583 G4double dEdx;
676 G4bool isOut; 584 G4bool isOut;
677 585
678 if (scaledKineticEnergy<t->theLowestKineticE << 586 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
679 587
680 dEdx = std::sqrt(scaledKineticEnergy/t->t << 588 dEdx = sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)
681 *(*dEdxTable)(materialIndex)->GetV 589 *(*dEdxTable)(materialIndex)->GetValue(
682 t->theLowestKineticEnergy,isOut) << 590 t.theLowestKineticEnergy,isOut);
683 591
684 } else if (scaledKineticEnergy>t->theHighest << 592 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
685 593
686 dEdx = (*dEdxTable)(materialIndex)->GetVa 594 dEdx = (*dEdxTable)(materialIndex)->GetValue(
687 t->theHighestKineticEnergy,isOut); << 595 t.theHighestKineticEnergy,isOut);
688 596
689 } else { 597 } else {
690 << 598
691 dEdx = (*dEdxTable)(materialIndex)->GetV 599 dEdx = (*dEdxTable)(materialIndex)->GetValue(
692 scaledKineticEnergy, 600 scaledKineticEnergy,isOut) ;
693 601
694 } 602 }
695 603
696 return dEdx*Chargesquare; 604 return dEdx*Chargesquare;
697 } 605 }
698 606
699 //....oooOO0OOooo........oooOO0OOooo........oo 607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
700 608
701 G4double G4EnergyLossTables::GetPreciseRangeF 609 G4double G4EnergyLossTables::GetPreciseRangeFromEnergy(
702 const G4ParticleDefinition *aParticle, 610 const G4ParticleDefinition *aParticle,
703 G4double KineticEnergy, 611 G4double KineticEnergy,
704 const G4Material *aMaterial) 612 const G4Material *aMaterial)
705 { 613 {
706 if (!t) t = new G4EnergyLossTablesHelper; << 614 if( aParticle != lastParticle)
707 <<
708 CPRWarning(); <<
709 if( aParticle != (const G4ParticleDefinition <<
710 { <<
711 *t= GetTables(aParticle); <<
712 lastParticle = (G4ParticleDefinition*) aPa <<
713 Chargesquare = (aParticle->GetPDGCharge()) <<
714 (aParticle->GetPDGCharge()) <<
715 QQPositron ; <<
716 oldIndex = -1 ; <<
717 } <<
718 const G4PhysicsTable* rangeTable= t->theRang <<
719 const G4PhysicsTable* dEdxTable= t->theDEDX <<
720 if (!rangeTable) { <<
721 ParticleHaveNoLoss(aParticle,"Range"); <<
722 return 0.0; <<
723 } <<
724 G4int materialIndex = (G4int)aMaterial->GetI <<
725 <<
726 G4double Thighr = t->theHighestKineticEnergy <<
727 (*rangeTable)(materialIndex <<
728 GetLowEdgeEnergy(1) ; <<
729 <<
730 G4double scaledKineticEnergy = KineticEnergy <<
731 G4double Range; <<
732 G4bool isOut; <<
733 <<
734 if (scaledKineticEnergy<t->theLowestKineticE <<
735 <<
736 Range = std::sqrt(scaledKineticEnergy/t->t <<
737 (*rangeTable)(materialIndex)->GetV <<
738 t->theLowestKineticEnergy,isOut) <<
739 <<
740 } else if (scaledKineticEnergy>Thighr) { <<
741 <<
742 Range = (*rangeTable)(materialIndex)->GetV <<
743 Thighr,isOut)+ <<
744 (scaledKineticEnergy-Thighr)/ <<
745 (*dEdxTable)(materialIndex)->GetVa <<
746 Thighr,isOut); <<
747 <<
748 } else { <<
749 <<
750 Range = (*rangeTable)(materialIndex)->Get <<
751 scaledKineticEnergy,isO <<
752 <<
753 } <<
754 <<
755 return Range/(Chargesquare*t->theMassRatio); <<
756 } <<
757 <<
758 //....oooOO0OOooo........oooOO0OOooo........oo <<
759 <<
760 G4double G4EnergyLossTables::GetDEDX( <<
761 const G4ParticleDefinition *aParticle, <<
762 G4double KineticEnergy, <<
763 const G4MaterialCutsCouple *couple, <<
764 G4bool check) <<
765 { <<
766 if (!t) t = new G4EnergyLossTablesHelper; <<
767 <<
768 if(aParticle != (const G4ParticleDefinition* <<
769 { 615 {
770 *t= GetTables(aParticle); << 616 t= GetTables(aParticle);
771 lastParticle = (G4ParticleDefinition*) aPa << 617 lastParticle = aParticle;
772 Chargesquare = (aParticle->GetPDGCharge()) <<
773 (aParticle->GetPDGCharge()) <<
774 QQPositron ; <<
775 oldIndex = -1 ; <<
776 } <<
777 const G4PhysicsTable* dEdxTable= t->theDEDX <<
778 <<
779 if (!dEdxTable ) { <<
780 if (check) return G4LossTableManager::Inst <<
781 else ParticleHaveNoLoss(aParticle, " <<
782 return 0.0; <<
783 } <<
784 <<
785 G4int materialIndex = couple->GetIndex(); <<
786 G4double scaledKineticEnergy = KineticEnergy <<
787 G4double dEdx; <<
788 G4bool isOut; <<
789 <<
790 if (scaledKineticEnergy<t->theLowestKineticE <<
791 <<
792 dEdx =(*dEdxTable)(materialIndex)->GetVal <<
793 t->theLowestKineticEnergy,isOut) <<
794 *std::sqrt(scaledKineticEnergy/t->t <<
795 <<
796 } else if (scaledKineticEnergy>t->theHighest <<
797 <<
798 dEdx = (*dEdxTable)(materialIndex)->GetVa <<
799 t->theHighestKineticEnergy,isOut); <<
800 <<
801 } else { <<
802 <<
803 dEdx = (*dEdxTable)(materialIndex)->GetVal <<
804 scaledKineticEnergy,isOut); <<
805 <<
806 } <<
807 <<
808 return dEdx*Chargesquare; <<
809 } <<
810 <<
811 //....oooOO0OOooo........oooOO0OOooo........oo <<
812 <<
813 G4double G4EnergyLossTables::GetRange( <<
814 const G4ParticleDefinition *aParticle, <<
815 G4double KineticEnergy, <<
816 const G4MaterialCutsCouple *couple, <<
817 G4bool check) <<
818 { <<
819 if (!t) t = new G4EnergyLossTablesHelper; <<
820 <<
821 if(aParticle != (const G4ParticleDefinition* <<
822 { <<
823 *t= GetTables(aParticle); <<
824 lastParticle = (G4ParticleDefinition*) aPa <<
825 Chargesquare = (aParticle->GetPDGCharge()) 618 Chargesquare = (aParticle->GetPDGCharge())*
826 (aParticle->GetPDGCharge()) 619 (aParticle->GetPDGCharge())/
827 QQPositron ; 620 QQPositron ;
828 oldIndex = -1 ; 621 oldIndex = -1 ;
829 } 622 }
830 const G4PhysicsTable* rangeTable= t->theRang << 623 const G4PhysicsTable* rangeTable= t.theRangeTable;
831 const G4PhysicsTable* dEdxTable= t->theDEDX << 624 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
832 if (!rangeTable) { <<
833 if(check) return G4LossTableManager::Insta <<
834 else return DBL_MAX; <<
835 //ParticleHaveNoLoss(aParticle,"Range"); <<
836 } <<
837 625
838 G4int materialIndex = couple->GetIndex(); << 626 G4int materialIndex = aMaterial->GetIndex();
839 G4double scaledKineticEnergy = KineticEnergy <<
840 G4double Range; <<
841 G4bool isOut; <<
842 <<
843 if (scaledKineticEnergy<t->theLowestKineticE <<
844 <<
845 Range = std::sqrt(scaledKineticEnergy/t->t <<
846 (*rangeTable)(materialIndex)->GetV <<
847 t->theLowestKineticEnergy,isOut) <<
848 627
849 } else if (scaledKineticEnergy>t->theHighest << 628 G4double Thighr = t.theHighestKineticEnergy*t.theLowestKineticEnergy/
850 <<
851 Range = (*rangeTable)(materialIndex)->GetV <<
852 t->theHighestKineticEnergy,isOut)+ <<
853 (scaledKineticEnergy-t->theHighest <<
854 (*dEdxTable)(materialIndex)->GetVa <<
855 t->theHighestKineticEnergy,isOut <<
856 <<
857 } else { <<
858 <<
859 Range = (*rangeTable)(materialIndex)->GetV <<
860 scaledKineticEnergy,isOut); <<
861 <<
862 } <<
863 <<
864 return Range/(Chargesquare*t->theMassRatio); <<
865 } <<
866 <<
867 //....oooOO0OOooo........oooOO0OOooo........oo <<
868 <<
869 G4double G4EnergyLossTables::GetPreciseEnergyF <<
870 const G4P <<
871 G4d <<
872 const G4M <<
873 G4bool check) <<
874 // it returns the value of the kinetic energy <<
875 { <<
876 if (!t) t = new G4EnergyLossTablesHelper; <<
877 <<
878 if( aParticle != (const G4ParticleDefinition <<
879 { <<
880 *t= GetTables(aParticle); <<
881 lastParticle = (G4ParticleDefinition*) aPa <<
882 Chargesquare = (aParticle->GetPDGCharge()) <<
883 (aParticle->GetPDGCharge()) <<
884 QQPositron ; <<
885 oldIndex = -1 ; <<
886 } <<
887 const G4PhysicsTable* dEdxTable= t->theDEDX <<
888 const G4PhysicsTable* inverseRangeTable= t- <<
889 <<
890 if (!inverseRangeTable) { <<
891 if(check) return G4LossTableManager::Insta <<
892 else return DBL_MAX; <<
893 // else ParticleHaveNoLoss(aPartic <<
894 } <<
895 <<
896 G4double scaledrange,scaledKineticEnergy ; <<
897 G4bool isOut ; <<
898 <<
899 G4int materialIndex = couple->GetIndex() ; <<
900 <<
901 if(materialIndex != oldIndex) <<
902 { <<
903 oldIndex = materialIndex ; <<
904 rmin = (*inverseRangeTable)(materialIndex) <<
905 GetLowEdgeEnergy <<
906 rmax = (*inverseRangeTable)(materialIndex) <<
907 GetLowEdgeEnergy(t->theNumb <<
908 Thigh = (*inverseRangeTable)(materialIndex <<
909 GetValue(rmax,is <<
910 } <<
911 <<
912 scaledrange = range*Chargesquare*t->theMassR <<
913 <<
914 if(scaledrange < rmin) <<
915 { <<
916 scaledKineticEnergy = t->theLowestKineticE <<
917 scaledrange*scaledrange/(rm <<
918 } <<
919 else <<
920 { <<
921 if(scaledrange < rmax) <<
922 { <<
923 scaledKineticEnergy = (*inverseRangeTabl <<
924 GetValue( scaled <<
925 } <<
926 else <<
927 { <<
928 scaledKineticEnergy = Thigh + <<
929 (scaledrange-rmax)* <<
930 (*dEdxTable)(materialInd <<
931 GetValue(Thig <<
932 } <<
933 } <<
934 <<
935 return scaledKineticEnergy/t->theMassRatio ; <<
936 } <<
937 <<
938 //....oooOO0OOooo........oooOO0OOooo........oo <<
939 <<
940 G4double G4EnergyLossTables::GetPreciseDEDX( <<
941 const G4ParticleDefinition *aParticle, <<
942 G4double KineticEnergy, <<
943 const G4MaterialCutsCouple *couple) <<
944 { <<
945 if (!t) t = new G4EnergyLossTablesHelper; <<
946 <<
947 if( aParticle != (const G4ParticleDefinition <<
948 { <<
949 *t= GetTables(aParticle); <<
950 lastParticle = (G4ParticleDefinition*) aPa <<
951 Chargesquare = (aParticle->GetPDGCharge()) <<
952 (aParticle->GetPDGCharge()) <<
953 QQPositron ; <<
954 oldIndex = -1 ; <<
955 } <<
956 const G4PhysicsTable* dEdxTable= t->theDEDX <<
957 if ( !dEdxTable ) <<
958 return G4LossTableManager::Instance()->Get <<
959 <<
960 G4int materialIndex = couple->GetIndex(); <<
961 G4double scaledKineticEnergy = KineticEnergy <<
962 G4double dEdx; <<
963 G4bool isOut; <<
964 <<
965 if (scaledKineticEnergy<t->theLowestKineticE <<
966 <<
967 dEdx = std::sqrt(scaledKineticEnergy/t->t <<
968 *(*dEdxTable)(materialIndex)->GetV <<
969 t->theLowestKineticEnergy,isOut) <<
970 <<
971 } else if (scaledKineticEnergy>t->theHighest <<
972 <<
973 dEdx = (*dEdxTable)(materialIndex)->GetVa <<
974 t->theHighestKineticEnergy,isOut); <<
975 <<
976 } else { <<
977 <<
978 dEdx = (*dEdxTable)(materialIndex)->GetV <<
979 scaledKineticEnergy, <<
980 <<
981 } <<
982 <<
983 return dEdx*Chargesquare; <<
984 } <<
985 <<
986 //....oooOO0OOooo........oooOO0OOooo........oo <<
987 <<
988 G4double G4EnergyLossTables::GetPreciseRangeFr <<
989 const G4ParticleDefinition *aParticle, <<
990 G4double KineticEnergy, <<
991 const G4MaterialCutsCouple *couple) <<
992 { <<
993 if (!t) t = new G4EnergyLossTablesHelper; <<
994 <<
995 if( aParticle != (const G4ParticleDefinition <<
996 { <<
997 *t= GetTables(aParticle); <<
998 lastParticle = (G4ParticleDefinition*) aPa <<
999 Chargesquare = (aParticle->GetPDGCharge()) <<
1000 (aParticle->GetPDGCharge() <<
1001 QQPositron ; <<
1002 oldIndex = -1 ; <<
1003 } <<
1004 const G4PhysicsTable* rangeTable= t->theRan <<
1005 const G4PhysicsTable* dEdxTable= t->theDED <<
1006 if ( !dEdxTable || !rangeTable) <<
1007 return G4LossTableManager::Instance()->Ge <<
1008 <<
1009 G4int materialIndex = couple->GetIndex(); <<
1010 <<
1011 G4double Thighr = t->theHighestKineticEnerg <<
1012 (*rangeTable)(materialInde 629 (*rangeTable)(materialIndex)->
1013 GetLowEdgeEnergy(1) ; 630 GetLowEdgeEnergy(1) ;
1014 631
1015 G4double scaledKineticEnergy = KineticEnerg << 632 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
1016 G4double Range; 633 G4double Range;
1017 G4bool isOut; 634 G4bool isOut;
1018 635
1019 if (scaledKineticEnergy<t->theLowestKinetic << 636 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
1020 637
1021 Range = std::sqrt(scaledKineticEnergy/t-> << 638 Range = sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)*
1022 (*rangeTable)(materialIndex)->Get 639 (*rangeTable)(materialIndex)->GetValue(
1023 t->theLowestKineticEnergy,isOut << 640 t.theLowestKineticEnergy,isOut);
1024 641
1025 } else if (scaledKineticEnergy>Thighr) { 642 } else if (scaledKineticEnergy>Thighr) {
1026 643
1027 Range = (*rangeTable)(materialIndex)->Get 644 Range = (*rangeTable)(materialIndex)->GetValue(
1028 Thighr,isOut)+ 645 Thighr,isOut)+
1029 (scaledKineticEnergy-Thighr)/ 646 (scaledKineticEnergy-Thighr)/
1030 (*dEdxTable)(materialIndex)->GetV 647 (*dEdxTable)(materialIndex)->GetValue(
1031 Thighr,isOut); 648 Thighr,isOut);
1032 649
1033 } else { 650 } else {
1034 << 651
1035 Range = (*rangeTable)(materialIndex)->Ge 652 Range = (*rangeTable)(materialIndex)->GetValue(
1036 scaledKineticEnergy,is 653 scaledKineticEnergy,isOut) ;
1037 654
1038 } 655 }
1039 656
1040 return Range/(Chargesquare*t->theMassRatio) << 657 return Range/(Chargesquare*t.theMassRatio);
1041 } 658 }
1042 659
1043 //....oooOO0OOooo........oooOO0OOooo........o 660 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1044 661
1045 void G4EnergyLossTables::CPRWarning() <<
1046 { <<
1047 if (let_counter < let_max_num_warnings) { <<
1048 662
1049 G4cout << G4endl; <<
1050 G4cout << "##### G4EnergyLossTable WARNIN <<
1051 G4cout << "##### RESULTS ARE NOT GARANTEE <<
1052 G4cout << "##### Please, substitute G4Mat <<
1053 G4cout << "##### Obsolete interface will <<
1054 G4cout << G4endl; <<
1055 let_counter++; <<
1056 663
1057 } else if (let_counter == let_max_num_warni <<
1058 664
1059 G4cout << "##### G4EnergyLossTable WARNIN <<
1060 let_counter++; <<
1061 } <<
1062 } <<
1063 <<
1064 //....oooOO0OOooo........oooOO0OOooo........o <<
1065 <<
1066 void <<
1067 G4EnergyLossTables::ParticleHaveNoLoss(const <<
1068 const G4String& /*q*/) <<
1069 { <<
1070 } <<
1071 <<
1072 //....oooOO0OOooo........oooOO0OOooo........o <<
1073 665