| 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 // 26 //
>> 27 // $Id$
>> 28 // GEANT4 tag $Name: geant4-09-01-ref-00 $
27 // 29 //
28 // ------------------------------------------- 30 // -----------------------------------------------------------
29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
30 // 32 //
31 // History: based on object model of 33 // History: based on object model of
32 // 2nd December 1995, G.Cosmo 34 // 2nd December 1995, G.Cosmo
33 // ---------- G4eLowEnergyLoss physics pr 35 // ---------- G4eLowEnergyLoss physics process -----------
34 // by Laszlo Urban, 20 March 19 36 // by Laszlo Urban, 20 March 1997
35 // ******************************************* 37 // **************************************************************
36 // It calculates the energy loss of e+/e-. 38 // It calculates the energy loss of e+/e-.
37 // ------------------------------------------- 39 // --------------------------------------------------------------
38 // 40 //
39 // 08-05-97: small changes by L.Urban 41 // 08-05-97: small changes by L.Urban
40 // 27-05-98: several bugs and inconsistencies 42 // 27-05-98: several bugs and inconsistencies are corrected,
41 // new table (the inverse of the ran 43 // new table (the inverse of the range table) added ,
42 // AlongStepDoit uses now this new t 44 // AlongStepDoit uses now this new table. L.Urban
43 // 08-09-98: cleanup 45 // 08-09-98: cleanup
44 // 24-09-98: rndmStepFlag false by default (no 46 // 24-09-98: rndmStepFlag false by default (no randomization of the step)
45 // 14-10-98: messenger file added. 47 // 14-10-98: messenger file added.
46 // 16-10-98: public method SetStepFunction() 48 // 16-10-98: public method SetStepFunction()
47 // 20-01-99: important correction in AlongStep 49 // 20-01-99: important correction in AlongStepDoIt , L.Urban
48 // 10/02/00 modifications , new e.m. structur 50 // 10/02/00 modifications , new e.m. structure, L.Urban
49 // 11/04/00: Bug fix in dE/dx fluctuation simu 51 // 11/04/00: Bug fix in dE/dx fluctuation simulation, Veronique Lefebure
50 // 19-09-00 change of fluctuation sampling V. 52 // 19-09-00 change of fluctuation sampling V.Ivanchenko
51 // 20/09/00 update fluctuations V.Ivanchenko 53 // 20/09/00 update fluctuations V.Ivanchenko
52 // 18/10/01 add fluorescence AlongStepDoIt V. 54 // 18/10/01 add fluorescence AlongStepDoIt V.Ivanchenko
53 // 18/10/01 Revision to improve code quality 55 // 18/10/01 Revision to improve code quality and consistency with design, MGP
54 // 19/10/01 update according to new design, V 56 // 19/10/01 update according to new design, V.Ivanchenko
55 // 24/10/01 MGP - Protection against negative 57 // 24/10/01 MGP - Protection against negative energy loss in AlongStepDoIt
56 // 26/10/01 VI Clean up access to deexcitatio 58 // 26/10/01 VI Clean up access to deexcitation
57 // 23/11/01 VI Move static member-functions f 59 // 23/11/01 VI Move static member-functions from header to source
58 // 28/05/02 VI Remove flag fStopAndKill 60 // 28/05/02 VI Remove flag fStopAndKill
59 // 03/06/02 MGP - Restore fStopAndKill 61 // 03/06/02 MGP - Restore fStopAndKill
60 // 28/10/02 VI Optimal binning for dE/dx 62 // 28/10/02 VI Optimal binning for dE/dx
61 // 21/01/03 VI cut per region 63 // 21/01/03 VI cut per region
62 // 01/06/04 VI check if stopped particle has 64 // 01/06/04 VI check if stopped particle has AtRest processes
63 // 65 //
64 // ------------------------------------------- 66 // --------------------------------------------------------------
65 67
66 #include "G4eLowEnergyLoss.hh" 68 #include "G4eLowEnergyLoss.hh"
67 #include "G4SystemOfUnits.hh" 69 #include "G4SystemOfUnits.hh"
68 #include "" 70 #include "G4EnergyLossMessenger.hh"
69 #include "G4Poisson.hh" 71 #include "G4Poisson.hh"
70 #include "G4ProductionCutsTable.hh" 72 #include "G4ProductionCutsTable.hh"
71 73
72 // 74 //
73 75
74 // Initialisation of static data members 76 // Initialisation of static data members
75 // ------------------------------------- 77 // -------------------------------------
76 // Contributing processes : ion.loss + soft br 78 // Contributing processes : ion.loss + soft brems->NbOfProcesses is initialized
77 // to 2 . YOU DO NOT HAVE TO CHANGE this varia 79 // to 2 . YOU DO NOT HAVE TO CHANGE this variable for a 'normal' run.
78 // 80 //
79 // You have to change NbOfProcesses if you inv 81 // You have to change NbOfProcesses if you invent a new process contributing
80 // to the continuous energy loss. 82 // to the continuous energy loss.
81 // The NbOfProcesses data member can be change 83 // The NbOfProcesses data member can be changed using the (public static)
82 // functions Get/Set/Plus/MinusNbOfProcesses ( 84 // functions Get/Set/Plus/MinusNbOfProcesses (see G4eLowEnergyLoss.hh)
83 85
84 G4int G4eLowEnergyLoss::NbOfProcess 86 G4int G4eLowEnergyLoss::NbOfProcesses = 2;
85 87
86 G4int G4eLowEnergyLoss::CounterOfEl 88 G4int G4eLowEnergyLoss::CounterOfElectronProcess = 0;
87 G4int G4eLowEnergyLoss::CounterOfPo 89 G4int G4eLowEnergyLoss::CounterOfPositronProcess = 0;
88 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfE 90 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfElectronProcess =
89 new 91 new G4PhysicsTable*[10];
90 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfP 92 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfPositronProcess =
91 new 93 new G4PhysicsTable*[10];
92 94
93 95
94 G4PhysicsTable* G4eLowEnergyLoss::theDEDXElec 96 G4PhysicsTable* G4eLowEnergyLoss::theDEDXElectronTable = 0;
95 G4PhysicsTable* G4eLowEnergyLoss::theDEDXPosi 97 G4PhysicsTable* G4eLowEnergyLoss::theDEDXPositronTable = 0;
96 G4PhysicsTable* G4eLowEnergyLoss::theRangeEle 98 G4PhysicsTable* G4eLowEnergyLoss::theRangeElectronTable = 0;
97 G4PhysicsTable* G4eLowEnergyLoss::theRangePos 99 G4PhysicsTable* G4eLowEnergyLoss::theRangePositronTable = 0;
98 G4PhysicsTable* G4eLowEnergyLoss::theInverseR 100 G4PhysicsTable* G4eLowEnergyLoss::theInverseRangeElectronTable = 0;
99 G4PhysicsTable* G4eLowEnergyLoss::theInverseR 101 G4PhysicsTable* G4eLowEnergyLoss::theInverseRangePositronTable = 0;
100 G4PhysicsTable* G4eLowEnergyLoss::theLabTimeE 102 G4PhysicsTable* G4eLowEnergyLoss::theLabTimeElectronTable = 0;
101 G4PhysicsTable* G4eLowEnergyLoss::theLabTimeP 103 G4PhysicsTable* G4eLowEnergyLoss::theLabTimePositronTable = 0;
102 G4PhysicsTable* G4eLowEnergyLoss::theProperTi 104 G4PhysicsTable* G4eLowEnergyLoss::theProperTimeElectronTable = 0;
103 G4PhysicsTable* G4eLowEnergyLoss::theProperTi 105 G4PhysicsTable* G4eLowEnergyLoss::theProperTimePositronTable = 0;
104 106
105 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo 107 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCoeffATable = 0;
106 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo 108 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCoeffBTable = 0;
107 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo 109 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCoeffCTable = 0;
108 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo 110 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCoeffATable = 0;
109 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo 111 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCoeffBTable = 0;
110 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo 112 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCoeffCTable = 0;
111 113
112 G4double G4eLowEnergyLoss::LowerBoundE 114 G4double G4eLowEnergyLoss::LowerBoundEloss = 10.*eV ;
113 G4double G4eLowEnergyLoss::UpperBoundE 115 G4double G4eLowEnergyLoss::UpperBoundEloss = 100.*GeV ;
114 G4int G4eLowEnergyLoss::NbinEloss = 116 G4int G4eLowEnergyLoss::NbinEloss = 360 ;
115 G4double G4eLowEnergyLoss::RTable ; 117 G4double G4eLowEnergyLoss::RTable ;
116 G4double G4eLowEnergyLoss::LOGRTable ; 118 G4double G4eLowEnergyLoss::LOGRTable ;
117 119
118 120
119 G4EnergyLossMessenger* G4eLowEnergyLoss::eLoss 121 G4EnergyLossMessenger* G4eLowEnergyLoss::eLossMessenger = 0;
120 122
121 // 123 //
122 124
123 // constructor and destructor 125 // constructor and destructor
124 126
125 G4eLowEnergyLoss::G4eLowEnergyLoss(const G4Str 127 G4eLowEnergyLoss::G4eLowEnergyLoss(const G4String& processName)
126 : G4RDVeLowEnergyLoss (processName), 128 : G4RDVeLowEnergyLoss (processName),
127 theLossTable(0), 129 theLossTable(0),
128 MinKineticEnergy(1.*eV), 130 MinKineticEnergy(1.*eV),
129 Charge(-1.), 131 Charge(-1.),
130 lastCharge(0.), 132 lastCharge(0.),
131 theDEDXTable(0), 133 theDEDXTable(0),
132 CounterOfProcess(0), 134 CounterOfProcess(0),
133 RecorderOfProcess(0), 135 RecorderOfProcess(0),
134 fdEdx(0), 136 fdEdx(0),
135 fRangeNow(0), 137 fRangeNow(0),
136 linLossLimit(0.05), 138 linLossLimit(0.05),
137 theFluo(false) 139 theFluo(false)
138 { 140 {
139 141
140 //create (only once) EnergyLoss messenger 142 //create (only once) EnergyLoss messenger
141 if(!eLossMessenger) eLossMessenger = new G4En 143 if(!eLossMessenger) eLossMessenger = new G4EnergyLossMessenger();
142 } 144 }
143 145
144 // 146 //
145 147
146 G4eLowEnergyLoss::~G4eLowEnergyLoss() 148 G4eLowEnergyLoss::~G4eLowEnergyLoss()
147 { 149 {
148 if (theLossTable) 150 if (theLossTable)
149 { 151 {
150 theLossTable->clearAndDestroy(); 152 theLossTable->clearAndDestroy();
151 delete theLossTable; 153 delete theLossTable;
152 } 154 }
153 } 155 }
154 156
155 void G4eLowEnergyLoss::SetNbOfProcesses(G4int 157 void G4eLowEnergyLoss::SetNbOfProcesses(G4int nb)
156 { 158 {
157 NbOfProcesses=nb; 159 NbOfProcesses=nb;
158 } 160 }
159 161
160 void G4eLowEnergyLoss::PlusNbOfProcesses() 162 void G4eLowEnergyLoss::PlusNbOfProcesses()
161 { 163 {
162 NbOfProcesses++; 164 NbOfProcesses++;
163 } 165 }
164 166
165 void G4eLowEnergyLoss::MinusNbOfProcesses() 167 void G4eLowEnergyLoss::MinusNbOfProcesses()
166 { 168 {
167 NbOfProcesses--; 169 NbOfProcesses--;
168 } 170 }
169 171
170 G4int G4eLowEnergyLoss::GetNbOfProcesses() 172 G4int G4eLowEnergyLoss::GetNbOfProcesses()
171 { 173 {
172 return NbOfProcesses; 174 return NbOfProcesses;
173 } 175 }
174 176
175 void G4eLowEnergyLoss::SetLowerBoundEloss(G4do 177 void G4eLowEnergyLoss::SetLowerBoundEloss(G4double val)
176 { 178 {
177 LowerBoundEloss=val; 179 LowerBoundEloss=val;
178 } 180 }
179 181
180 void G4eLowEnergyLoss::SetUpperBoundEloss(G4do 182 void G4eLowEnergyLoss::SetUpperBoundEloss(G4double val)
181 { 183 {
182 UpperBoundEloss=val; 184 UpperBoundEloss=val;
183 } 185 }
184 186
185 void G4eLowEnergyLoss::SetNbinEloss(G4int nb) 187 void G4eLowEnergyLoss::SetNbinEloss(G4int nb)
186 { 188 {
187 NbinEloss=nb; 189 NbinEloss=nb;
188 } 190 }
189 191
190 G4double G4eLowEnergyLoss::GetLowerBoundEloss( 192 G4double G4eLowEnergyLoss::GetLowerBoundEloss()
191 { 193 {
192 return LowerBoundEloss; 194 return LowerBoundEloss;
193 } 195 }
194 196
195 G4double G4eLowEnergyLoss::GetUpperBoundEloss( 197 G4double G4eLowEnergyLoss::GetUpperBoundEloss()
196 { 198 {
197 return UpperBoundEloss; 199 return UpperBoundEloss;
198 } 200 }
199 201
200 G4int G4eLowEnergyLoss::GetNbinEloss() 202 G4int G4eLowEnergyLoss::GetNbinEloss()
201 { 203 {
202 return NbinEloss; 204 return NbinEloss;
203 } 205 }
204 // 206 //
205 207
206 void G4eLowEnergyLoss::BuildDEDXTable( 208 void G4eLowEnergyLoss::BuildDEDXTable(
207 const G4ParticleDefin 209 const G4ParticleDefinition& aParticleType)
208 { 210 {
209 ParticleMass = aParticleType.GetPDGMass(); 211 ParticleMass = aParticleType.GetPDGMass();
210 Charge = aParticleType.GetPDGCharge()/eplus; 212 Charge = aParticleType.GetPDGCharge()/eplus;
211 213
212 // calculate data members LOGRTable,RTable 214 // calculate data members LOGRTable,RTable first
213 215
214 G4double lrate = std::log(UpperBoundEloss/Lo 216 G4double lrate = std::log(UpperBoundEloss/LowerBoundEloss);
215 LOGRTable=lrate/NbinEloss; 217 LOGRTable=lrate/NbinEloss;
216 RTable =std::exp(LOGRTable); 218 RTable =std::exp(LOGRTable);
217 // Build energy loss table as a sum of the e 219 // Build energy loss table as a sum of the energy loss due to the
218 // different processes. 220 // different processes.
219 // 221 //
220 222
221 const G4ProductionCutsTable* theCoupleTable= 223 const G4ProductionCutsTable* theCoupleTable=
222 G4ProductionCutsTable::GetProductionCu 224 G4ProductionCutsTable::GetProductionCutsTable();
223 size_t numOfCouples = theCoupleTable->GetTab 225 size_t numOfCouples = theCoupleTable->GetTableSize();
224 226
225 // create table for the total energy loss 227 // create table for the total energy loss
226 228
227 if (&aParticleType==G4Electron::Electron()) 229 if (&aParticleType==G4Electron::Electron())
228 { 230 {
229 RecorderOfProcess=RecorderOfElectronProc 231 RecorderOfProcess=RecorderOfElectronProcess;
230 CounterOfProcess=CounterOfElectronProces 232 CounterOfProcess=CounterOfElectronProcess;
231 if (CounterOfProcess == NbOfProcesses) 233 if (CounterOfProcess == NbOfProcesses)
232 { 234 {
233 if (theDEDXElectronTable) 235 if (theDEDXElectronTable)
234 { 236 {
235 theDEDXElectronTable->clearAndDes 237 theDEDXElectronTable->clearAndDestroy();
236 delete theDEDXElectronTable; 238 delete theDEDXElectronTable;
237 } 239 }
238 theDEDXElectronTable = new G4PhysicsT 240 theDEDXElectronTable = new G4PhysicsTable(numOfCouples);
239 theDEDXTable = theDEDXElectronTable; 241 theDEDXTable = theDEDXElectronTable;
240 } 242 }
241 } 243 }
242 if (&aParticleType==G4Positron::Positron()) 244 if (&aParticleType==G4Positron::Positron())
243 { 245 {
244 RecorderOfProcess=RecorderOfPositronProce 246 RecorderOfProcess=RecorderOfPositronProcess;
245 CounterOfProcess=CounterOfPositronProcess 247 CounterOfProcess=CounterOfPositronProcess;
246 if (CounterOfProcess == NbOfProcesses) 248 if (CounterOfProcess == NbOfProcesses)
247 { 249 {
248 if (theDEDXPositronTable) 250 if (theDEDXPositronTable)
249 { 251 {
250 theDEDXPositronTable->clearAndDest 252 theDEDXPositronTable->clearAndDestroy();
251 delete theDEDXPositronTable; 253 delete theDEDXPositronTable;
252 } 254 }
253 theDEDXPositronTable = new G4PhysicsTa 255 theDEDXPositronTable = new G4PhysicsTable(numOfCouples);
254 theDEDXTable = theDEDXPositronTable; 256 theDEDXTable = theDEDXPositronTable;
255 } 257 }
256 } 258 }
257 259
258 if (CounterOfProcess == NbOfProcesses) 260 if (CounterOfProcess == NbOfProcesses)
259 { 261 {
260 // fill the tables 262 // fill the tables
261 // loop for materials 263 // loop for materials
262 G4double LowEdgeEnergy , Value; 264 G4double LowEdgeEnergy , Value;
263 G4bool isOutRange; 265 G4bool isOutRange;
264 G4PhysicsTable* pointer; 266 G4PhysicsTable* pointer;
265 267
266 for (size_t J=0; J<numOfCouples; J++) 268 for (size_t J=0; J<numOfCouples; J++)
267 { 269 {
268 // create physics vector and fill it 270 // create physics vector and fill it
269 271
270 G4PhysicsLogVector* aVector = new G4P 272 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
271 LowerBoundEloss, UpperBoun 273 LowerBoundEloss, UpperBoundEloss, NbinEloss);
272 274
273 // loop for the kinetic energy 275 // loop for the kinetic energy
274 for (G4int i=0; i<NbinEloss; i++) 276 for (G4int i=0; i<NbinEloss; i++)
275 { 277 {
276 LowEdgeEnergy = aVector->GetLowE 278 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
277 //here comes the sum of the diff 279 //here comes the sum of the different tables created by the
278 //processes (ionisation,bremsstr 280 //processes (ionisation,bremsstrahlung,etc...)
279 Value = 0.; 281 Value = 0.;
280 for (G4int process=0; process < 282 for (G4int process=0; process < NbOfProcesses; process++)
281 { 283 {
282 pointer= RecorderOfProcess[ 284 pointer= RecorderOfProcess[process];
283 Value += (*pointer)[J]->Get 285 Value += (*pointer)[J]->GetValue(LowEdgeEnergy,isOutRange);
284 } 286 }
285 287
286 aVector->PutValue(i,Value) ; 288 aVector->PutValue(i,Value) ;
287 } 289 }
288 290
289 theDEDXTable->insert(aVector) ; 291 theDEDXTable->insert(aVector) ;
290 292
291 } 293 }
292 294
293 295
294 //reset counter to zero 296 //reset counter to zero
295 if (&aParticleType==G4Electron::Electron( 297 if (&aParticleType==G4Electron::Electron()) CounterOfElectronProcess=0;
296 if (&aParticleType==G4Positron::Positron( 298 if (&aParticleType==G4Positron::Positron()) CounterOfPositronProcess=0;
297 299
298 ParticleMass = aParticleType.GetPDGMass() 300 ParticleMass = aParticleType.GetPDGMass();
299 301
300 if (&aParticleType==G4Electron::Electron( 302 if (&aParticleType==G4Electron::Electron())
301 { 303 {
302 // Build range table 304 // Build range table
303 theRangeElectronTable = BuildRangeTable 305 theRangeElectronTable = BuildRangeTable(theDEDXElectronTable,
304 306 theRangeElectronTable,
305 LowerBoundEloss, 307 LowerBoundEloss,UpperBoundEloss,NbinEloss);
306 308
307 // Build lab/proper time tables 309 // Build lab/proper time tables
308 theLabTimeElectronTable = BuildLabTimeT 310 theLabTimeElectronTable = BuildLabTimeTable(theDEDXElectronTable,
309 theLabTimeElectronTab 311 theLabTimeElectronTable,
310 LowerBoundEloss,Upper 312 LowerBoundEloss,UpperBoundEloss,NbinEloss);
311 theProperTimeElectronTable = BuildPrope 313 theProperTimeElectronTable = BuildProperTimeTable(theDEDXElectronTable,
312 theProperTimeElect 314 theProperTimeElectronTable,
313 LowerBoundEloss,Up 315 LowerBoundEloss,UpperBoundEloss,NbinEloss);
314 316
315 // Build coeff tables for the energy lo 317 // Build coeff tables for the energy loss calculation
316 theeRangeCoeffATable = BuildRangeCoeffA 318 theeRangeCoeffATable = BuildRangeCoeffATable(theRangeElectronTable,
317 theeRangeCoeffATa 319 theeRangeCoeffATable,
318 LowerBoundEloss,U 320 LowerBoundEloss,UpperBoundEloss,NbinEloss);
319 321
320 theeRangeCoeffBTable = BuildRangeCoeffB 322 theeRangeCoeffBTable = BuildRangeCoeffBTable(theRangeElectronTable,
321 theeRangeCoeffBTa 323 theeRangeCoeffBTable,
322 LowerBoundEloss,U 324 LowerBoundEloss,UpperBoundEloss,NbinEloss);
323 325
324 theeRangeCoeffCTable = BuildRangeCoeffC 326 theeRangeCoeffCTable = BuildRangeCoeffCTable(theRangeElectronTable,
325 theeRangeCoeffCTa 327 theeRangeCoeffCTable,
326 LowerBoundEloss,U 328 LowerBoundEloss,UpperBoundEloss,NbinEloss);
327 329
328 // invert the range table 330 // invert the range table
329 theInverseRangeElectronTable = BuildInv 331 theInverseRangeElectronTable = BuildInverseRangeTable(theRangeElectronTable,
330 theeRangeCoeffAT 332 theeRangeCoeffATable,
331 theeRangeCoeffBT 333 theeRangeCoeffBTable,
332 theeRangeCoeffCT 334 theeRangeCoeffCTable,
333 theInverseRangeE 335 theInverseRangeElectronTable,
334 LowerBoundEloss, 336 LowerBoundEloss,UpperBoundEloss,NbinEloss);
335 } 337 }
336 if (&aParticleType==G4Positron::Positron( 338 if (&aParticleType==G4Positron::Positron())
337 { 339 {
338 // Build range table 340 // Build range table
339 theRangePositronTable = BuildRangeTable 341 theRangePositronTable = BuildRangeTable(theDEDXPositronTable,
340 342 theRangePositronTable,
341 LowerBoundEloss, 343 LowerBoundEloss,UpperBoundEloss,NbinEloss);
342 344
343 345
344 // Build lab/proper time tables 346 // Build lab/proper time tables
345 theLabTimePositronTable = BuildLabTimeT 347 theLabTimePositronTable = BuildLabTimeTable(theDEDXPositronTable,
346 theLabTimePositronTab 348 theLabTimePositronTable,
347 LowerBoundEloss,Upper 349 LowerBoundEloss,UpperBoundEloss,NbinEloss);
348 theProperTimePositronTable = BuildPrope 350 theProperTimePositronTable = BuildProperTimeTable(theDEDXPositronTable,
349 theProperTimePosit 351 theProperTimePositronTable,
350 LowerBoundEloss,Up 352 LowerBoundEloss,UpperBoundEloss,NbinEloss);
351 353
352 // Build coeff tables for the energy lo 354 // Build coeff tables for the energy loss calculation
353 thepRangeCoeffATable = BuildRangeCoeffA 355 thepRangeCoeffATable = BuildRangeCoeffATable(theRangePositronTable,
354 thepRangeCoeffATa 356 thepRangeCoeffATable,
355 LowerBoundEloss,U 357 LowerBoundEloss,UpperBoundEloss,NbinEloss);
356 358
357 thepRangeCoeffBTable = BuildRangeCoeffB 359 thepRangeCoeffBTable = BuildRangeCoeffBTable(theRangePositronTable,
358 thepRangeCoeffBTa 360 thepRangeCoeffBTable,
359 LowerBoundEloss,U 361 LowerBoundEloss,UpperBoundEloss,NbinEloss);
360 362
361 thepRangeCoeffCTable = BuildRangeCoeffC 363 thepRangeCoeffCTable = BuildRangeCoeffCTable(theRangePositronTable,
362 thepRangeCoeffCTa 364 thepRangeCoeffCTable,
363 LowerBoundEloss,U 365 LowerBoundEloss,UpperBoundEloss,NbinEloss);
364 366
365 // invert the range table 367 // invert the range table
366 theInverseRangePositronTable = BuildInv 368 theInverseRangePositronTable = BuildInverseRangeTable(theRangePositronTable,
367 thepRangeCoeffAT 369 thepRangeCoeffATable,
368 thepRangeCoeffBT 370 thepRangeCoeffBTable,
369 thepRangeCoeffCT 371 thepRangeCoeffCTable,
370 theInverseRangeP 372 theInverseRangePositronTable,
371 LowerBoundEloss, 373 LowerBoundEloss,UpperBoundEloss,NbinEloss);
372 } 374 }
373 375
374 if(verboseLevel > 1) { 376 if(verboseLevel > 1) {
375 G4cout << (*theDEDXElectronTable) << G4 377 G4cout << (*theDEDXElectronTable) << G4endl;
376 } 378 }
377 379
378 380
379 // make the energy loss and the range tab 381 // make the energy loss and the range table available
380 G4EnergyLossTables::Register(&aParticleTy 382 G4EnergyLossTables::Register(&aParticleType,
381 (&aParticleType==G4Electron::Electron() 383 (&aParticleType==G4Electron::Electron())?
382 theDEDXElectronTable: theDEDXPositronTa 384 theDEDXElectronTable: theDEDXPositronTable,
383 (&aParticleType==G4Electron::Electron() 385 (&aParticleType==G4Electron::Electron())?
384 theRangeElectronTable: theRangePositron 386 theRangeElectronTable: theRangePositronTable,
385 (&aParticleType==G4Electron::Electron() 387 (&aParticleType==G4Electron::Electron())?
386 theInverseRangeElectronTable: theInvers 388 theInverseRangeElectronTable: theInverseRangePositronTable,
387 (&aParticleType==G4Electron::Electron() 389 (&aParticleType==G4Electron::Electron())?
388 theLabTimeElectronTable: theLabTimePosi 390 theLabTimeElectronTable: theLabTimePositronTable,
389 (&aParticleType==G4Electron::Electron() 391 (&aParticleType==G4Electron::Electron())?
390 theProperTimeElectronTable: theProperTi 392 theProperTimeElectronTable: theProperTimePositronTable,
391 LowerBoundEloss, UpperBoundEloss, 1.,Nb 393 LowerBoundEloss, UpperBoundEloss, 1.,NbinEloss);
392 } 394 }
393 } 395 }
394 396
395 // 397 //
396 398
397 G4VParticleChange* G4eLowEnergyLoss::AlongStep 399 G4VParticleChange* G4eLowEnergyLoss::AlongStepDoIt( const G4Track& trackData,
398 400 const G4Step& stepData)
399 { 401 {
400 // compute the energy loss after a Step 402 // compute the energy loss after a Step
401 403
402 static const G4double faclow = 1.5 ; 404 static const G4double faclow = 1.5 ;
403 405
404 // get particle and material pointers from t 406 // get particle and material pointers from trackData
405 const G4DynamicParticle* aParticle = trackDa 407 const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
406 G4double E = aParticle->GetKineticEnerg 408 G4double E = aParticle->GetKineticEnergy();
407 409
408 const G4MaterialCutsCouple* couple = trackDa 410 const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
409 411
410 G4double Step = stepData.GetStepLength(); 412 G4double Step = stepData.GetStepLength();
411 413
412 aParticleChange.Initialize(trackData); 414 aParticleChange.Initialize(trackData);
413 //fParticleChange.Initialize(trackData); 415 //fParticleChange.Initialize(trackData);
414 416
415 G4double MeanLoss, finalT; 417 G4double MeanLoss, finalT;
416 418
417 if (E < MinKineticEnergy) finalT = 0.; 419 if (E < MinKineticEnergy) finalT = 0.;
418 420
419 else if ( E< faclow*LowerBoundEloss) 421 else if ( E< faclow*LowerBoundEloss)
420 { 422 {
421 if (Step >= fRangeNow) finalT = 0.; 423 if (Step >= fRangeNow) finalT = 0.;
422 // else finalT = E*(1.-Step/fRangeNow) ; 424 // else finalT = E*(1.-Step/fRangeNow) ;
423 else finalT = E*(1.-std::sqrt(Step/fRangeN 425 else finalT = E*(1.-std::sqrt(Step/fRangeNow)) ;
424 } 426 }
425 427
426 else if (E>=UpperBoundEloss) finalT = E - St 428 else if (E>=UpperBoundEloss) finalT = E - Step*fdEdx;
427 429
428 else if (Step >= fRangeNow) finalT = 0.; 430 else if (Step >= fRangeNow) finalT = 0.;
429 431
430 else 432 else
431 { 433 {
432 if(Step/fRangeNow < linLossLimit) finalT = 434 if(Step/fRangeNow < linLossLimit) finalT = E-Step*fdEdx ;
433 else 435 else
434 { 436 {
435 if (Charge<0.) finalT = G4EnergyLossTabl 437 if (Charge<0.) finalT = G4EnergyLossTables::GetPreciseEnergyFromRange
436 (G4Electron::Elec 438 (G4Electron::Electron(),fRangeNow-Step,couple);
437 else finalT = G4EnergyLossTabl 439 else finalT = G4EnergyLossTables::GetPreciseEnergyFromRange
438 (G4Positron::Posi 440 (G4Positron::Positron(),fRangeNow-Step,couple);
439 } 441 }
440 } 442 }
441 443
442 if(finalT < MinKineticEnergy) finalT = 0. ; 444 if(finalT < MinKineticEnergy) finalT = 0. ;
443 445
444 MeanLoss = E-finalT ; 446 MeanLoss = E-finalT ;
445 447
446 //now the loss with fluctuation 448 //now the loss with fluctuation
447 if ((EnlossFlucFlag) && (finalT > 0.) && (fi 449 if ((EnlossFlucFlag) && (finalT > 0.) && (finalT < E)&&(E > LowerBoundEloss))
448 { 450 {
449 finalT = E-GetLossWithFluct(aParticle,coup 451 finalT = E-GetLossWithFluct(aParticle,couple,MeanLoss,Step);
450 if (finalT < 0.) finalT = 0.; 452 if (finalT < 0.) finalT = 0.;
451 } 453 }
452 454
453 // kill the particle if the kinetic energy < 455 // kill the particle if the kinetic energy <= 0
454 if (finalT <= 0. ) 456 if (finalT <= 0. )
455 { 457 {
456 finalT = 0.; 458 finalT = 0.;
457 if(Charge > 0.0) aParticleChange.ProposeTr 459 if(Charge > 0.0) aParticleChange.ProposeTrackStatus(fStopButAlive);
458 else aParticleChange.ProposeTr 460 else aParticleChange.ProposeTrackStatus(fStopAndKill);
459 } 461 }
460 462
461 G4double edep = E - finalT; 463 G4double edep = E - finalT;
462 464
463 aParticleChange.ProposeEnergy(finalT); 465 aParticleChange.ProposeEnergy(finalT);
464 466
465 // Deexcitation of ionised atoms 467 // Deexcitation of ionised atoms
466 std::vector<G4DynamicParticle*>* deexcitatio 468 std::vector<G4DynamicParticle*>* deexcitationProducts = 0;
467 if (theFluo) deexcitationProducts = Deexcite 469 if (theFluo) deexcitationProducts = DeexciteAtom(couple,E,edep);
468 470
469 size_t nSecondaries = 0; 471 size_t nSecondaries = 0;
470 if (deexcitationProducts != 0) nSecondaries 472 if (deexcitationProducts != 0) nSecondaries = deexcitationProducts->size();
471 aParticleChange.SetNumberOfSecondaries(nSeco 473 aParticleChange.SetNumberOfSecondaries(nSecondaries);
472 474
473 if (nSecondaries > 0) { 475 if (nSecondaries > 0) {
474 476
475 const G4StepPoint* preStep = stepData.GetP 477 const G4StepPoint* preStep = stepData.GetPreStepPoint();
476 const G4StepPoint* postStep = stepData.Get 478 const G4StepPoint* postStep = stepData.GetPostStepPoint();
477 G4ThreeVector r = preStep->GetPosition(); 479 G4ThreeVector r = preStep->GetPosition();
478 G4ThreeVector deltaR = postStep->GetPositi 480 G4ThreeVector deltaR = postStep->GetPosition();
479 deltaR -= r; 481 deltaR -= r;
480 G4double t = preStep->GetGlobalTime(); 482 G4double t = preStep->GetGlobalTime();
481 G4double deltaT = postStep->GetGlobalTime( 483 G4double deltaT = postStep->GetGlobalTime();
482 deltaT -= t; 484 deltaT -= t;
483 G4double time, q; 485 G4double time, q;
484 G4ThreeVector position; 486 G4ThreeVector position;
485 487
486 for (size_t i=0; i<nSecondaries; i++) { 488 for (size_t i=0; i<nSecondaries; i++) {
487 489
488 G4DynamicParticle* part = (*deexcitation 490 G4DynamicParticle* part = (*deexcitationProducts)[i];
489 if (part != 0) { 491 if (part != 0) {
490 G4double eSecondary = part->GetKinetic 492 G4double eSecondary = part->GetKineticEnergy();
491 edep -= eSecondary; 493 edep -= eSecondary;
492 if (edep > 0.) 494 if (edep > 0.)
493 { 495 {
494 q = G4UniformRand(); 496 q = G4UniformRand();
495 time = deltaT*q + t; 497 time = deltaT*q + t;
496 position = deltaR*q; 498 position = deltaR*q;
497 position += r; 499 position += r;
498 G4Track* newTrack = new G4Track(part, ti 500 G4Track* newTrack = new G4Track(part, time, position);
499 aParticleChange.AddSecondary(newTrack); 501 aParticleChange.AddSecondary(newTrack);
500 } 502 }
501 else 503 else
502 { 504 {
503 edep += eSecondary; 505 edep += eSecondary;
504 delete part; 506 delete part;
505 part = 0; 507 part = 0;
506 } 508 }
507 } 509 }
508 } 510 }
509 } 511 }
510 delete deexcitationProducts; 512 delete deexcitationProducts;
511 513
512 aParticleChange.ProposeLocalEnergyDeposit(ed 514 aParticleChange.ProposeLocalEnergyDeposit(edep);
513 515
514 return &aParticleChange; 516 return &aParticleChange;
515 } 517 }
516 518
517 // 519 //
518 520
519 521