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