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