| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 pr
34 // by Laszlo Urban, 20 March 19
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
41 // new table (the inverse of the ran
42 // AlongStepDoit uses now this new t
43 // 08-09-98: cleanup
44 // 24-09-98: rndmStepFlag false by default (no
45 // 14-10-98: messenger file added.
46 // 16-10-98: public method SetStepFunction()
47 // 20-01-99: important correction in AlongStep
48 // 10/02/00 modifications , new e.m. structur
49 // 11/04/00: Bug fix in dE/dx fluctuation simu
50 // 19-09-00 change of fluctuation sampling V.
51 // 20/09/00 update fluctuations V.Ivanchenko
52 // 18/10/01 add fluorescence AlongStepDoIt V.
53 // 18/10/01 Revision to improve code quality
54 // 19/10/01 update according to new design, V
55 // 24/10/01 MGP - Protection against negative
56 // 26/10/01 VI Clean up access to deexcitatio
57 // 23/11/01 VI Move static member-functions f
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
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 br
77 // to 2 . YOU DO NOT HAVE TO CHANGE this varia
78 //
79 // You have to change NbOfProcesses if you inv
80 // to the continuous energy loss.
81 // The NbOfProcesses data member can be change
82 // functions Get/Set/Plus/MinusNbOfProcesses (
83
84 G4int G4eLowEnergyLoss::NbOfProcess
85
86 G4int G4eLowEnergyLoss::CounterOfEl
87 G4int G4eLowEnergyLoss::CounterOfPo
88 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfE
89 new
90 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfP
91 new
92
93
94 G4PhysicsTable* G4eLowEnergyLoss::theDEDXElec
95 G4PhysicsTable* G4eLowEnergyLoss::theDEDXPosi
96 G4PhysicsTable* G4eLowEnergyLoss::theRangeEle
97 G4PhysicsTable* G4eLowEnergyLoss::theRangePos
98 G4PhysicsTable* G4eLowEnergyLoss::theInverseR
99 G4PhysicsTable* G4eLowEnergyLoss::theInverseR
100 G4PhysicsTable* G4eLowEnergyLoss::theLabTimeE
101 G4PhysicsTable* G4eLowEnergyLoss::theLabTimeP
102 G4PhysicsTable* G4eLowEnergyLoss::theProperTi
103 G4PhysicsTable* G4eLowEnergyLoss::theProperTi
104
105 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo
106 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo
107 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo
108 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo
109 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo
110 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo
111
112 G4double G4eLowEnergyLoss::LowerBoundE
113 G4double G4eLowEnergyLoss::UpperBoundE
114 G4int G4eLowEnergyLoss::NbinEloss =
115 G4double G4eLowEnergyLoss::RTable ;
116 G4double G4eLowEnergyLoss::LOGRTable ;
117
118
119 G4EnergyLossMessenger* G4eLowEnergyLoss::eLoss
120
121 //
122
123 // constructor and destructor
124
125 G4eLowEnergyLoss::G4eLowEnergyLoss(const G4Str
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 G4En
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
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(G4do
176 {
177 LowerBoundEloss=val;
178 }
179
180 void G4eLowEnergyLoss::SetUpperBoundEloss(G4do
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 G4ParticleDefin
208 {
209 ParticleMass = aParticleType.GetPDGMass();
210 Charge = aParticleType.GetPDGCharge()/eplus;
211
212 // calculate data members LOGRTable,RTable
213
214 G4double lrate = std::log(UpperBoundEloss/Lo
215 LOGRTable=lrate/NbinEloss;
216 RTable =std::exp(LOGRTable);
217 // Build energy loss table as a sum of the e
218 // different processes.
219 //
220
221 const G4ProductionCutsTable* theCoupleTable=
222 G4ProductionCutsTable::GetProductionCu
223 size_t numOfCouples = theCoupleTable->GetTab
224
225 // create table for the total energy loss
226
227 if (&aParticleType==G4Electron::Electron())
228 {
229 RecorderOfProcess=RecorderOfElectronProc
230 CounterOfProcess=CounterOfElectronProces
231 if (CounterOfProcess == NbOfProcesses)
232 {
233 if (theDEDXElectronTable)
234 {
235 theDEDXElectronTable->clearAndDes
236 delete theDEDXElectronTable;
237 }
238 theDEDXElectronTable = new G4PhysicsT
239 theDEDXTable = theDEDXElectronTable;
240 }
241 }
242 if (&aParticleType==G4Positron::Positron())
243 {
244 RecorderOfProcess=RecorderOfPositronProce
245 CounterOfProcess=CounterOfPositronProcess
246 if (CounterOfProcess == NbOfProcesses)
247 {
248 if (theDEDXPositronTable)
249 {
250 theDEDXPositronTable->clearAndDest
251 delete theDEDXPositronTable;
252 }
253 theDEDXPositronTable = new G4PhysicsTa
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 G4P
271 LowerBoundEloss, UpperBoun
272
273 // loop for the kinetic energy
274 for (G4int i=0; i<NbinEloss; i++)
275 {
276 LowEdgeEnergy = aVector->GetLowE
277 //here comes the sum of the diff
278 //processes (ionisation,bremsstr
279 Value = 0.;
280 for (G4int process=0; process <
281 {
282 pointer= RecorderOfProcess[
283 Value += (*pointer)[J]->Get
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(
296 if (&aParticleType==G4Positron::Positron(
297
298 ParticleMass = aParticleType.GetPDGMass()
299
300 if (&aParticleType==G4Electron::Electron(
301 {
302 // Build range table
303 theRangeElectronTable = BuildRangeTable
304
305 LowerBoundEloss,
306
307 // Build lab/proper time tables
308 theLabTimeElectronTable = BuildLabTimeT
309 theLabTimeElectronTab
310 LowerBoundEloss,Upper
311 theProperTimeElectronTable = BuildPrope
312 theProperTimeElect
313 LowerBoundEloss,Up
314
315 // Build coeff tables for the energy lo
316 theeRangeCoeffATable = BuildRangeCoeffA
317 theeRangeCoeffATa
318 LowerBoundEloss,U
319
320 theeRangeCoeffBTable = BuildRangeCoeffB
321 theeRangeCoeffBTa
322 LowerBoundEloss,U
323
324 theeRangeCoeffCTable = BuildRangeCoeffC
325 theeRangeCoeffCTa
326 LowerBoundEloss,U
327
328 // invert the range table
329 theInverseRangeElectronTable = BuildInv
330 theeRangeCoeffAT
331 theeRangeCoeffBT
332 theeRangeCoeffCT
333 theInverseRangeE
334 LowerBoundEloss,
335 }
336 if (&aParticleType==G4Positron::Positron(
337 {
338 // Build range table
339 theRangePositronTable = BuildRangeTable
340
341 LowerBoundEloss,
342
343
344 // Build lab/proper time tables
345 theLabTimePositronTable = BuildLabTimeT
346 theLabTimePositronTab
347 LowerBoundEloss,Upper
348 theProperTimePositronTable = BuildPrope
349 theProperTimePosit
350 LowerBoundEloss,Up
351
352 // Build coeff tables for the energy lo
353 thepRangeCoeffATable = BuildRangeCoeffA
354 thepRangeCoeffATa
355 LowerBoundEloss,U
356
357 thepRangeCoeffBTable = BuildRangeCoeffB
358 thepRangeCoeffBTa
359 LowerBoundEloss,U
360
361 thepRangeCoeffCTable = BuildRangeCoeffC
362 thepRangeCoeffCTa
363 LowerBoundEloss,U
364
365 // invert the range table
366 theInverseRangePositronTable = BuildInv
367 thepRangeCoeffAT
368 thepRangeCoeffBT
369 thepRangeCoeffCT
370 theInverseRangeP
371 LowerBoundEloss,
372 }
373
374 if(verboseLevel > 1) {
375 G4cout << (*theDEDXElectronTable) << G4
376 }
377
378
379 // make the energy loss and the range tab
380 G4EnergyLossTables::Register(&aParticleTy
381 (&aParticleType==G4Electron::Electron()
382 theDEDXElectronTable: theDEDXPositronTa
383 (&aParticleType==G4Electron::Electron()
384 theRangeElectronTable: theRangePositron
385 (&aParticleType==G4Electron::Electron()
386 theInverseRangeElectronTable: theInvers
387 (&aParticleType==G4Electron::Electron()
388 theLabTimeElectronTable: theLabTimePosi
389 (&aParticleType==G4Electron::Electron()
390 theProperTimeElectronTable: theProperTi
391 LowerBoundEloss, UpperBoundEloss, 1.,Nb
392 }
393 }
394
395 //
396
397 G4VParticleChange* G4eLowEnergyLoss::AlongStep
398
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 t
405 const G4DynamicParticle* aParticle = trackDa
406 G4double E = aParticle->GetKineticEnerg
407
408 const G4MaterialCutsCouple* couple = trackDa
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/fRangeN
424 }
425
426 else if (E>=UpperBoundEloss) finalT = E - St
427
428 else if (Step >= fRangeNow) finalT = 0.;
429
430 else
431 {
432 if(Step/fRangeNow < linLossLimit) finalT =
433 else
434 {
435 if (Charge<0.) finalT = G4EnergyLossTabl
436 (G4Electron::Elec
437 else finalT = G4EnergyLossTabl
438 (G4Positron::Posi
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.) && (fi
448 {
449 finalT = E-GetLossWithFluct(aParticle,coup
450 if (finalT < 0.) finalT = 0.;
451 }
452
453 // kill the particle if the kinetic energy <
454 if (finalT <= 0. )
455 {
456 finalT = 0.;
457 if(Charge > 0.0) aParticleChange.ProposeTr
458 else aParticleChange.ProposeTr
459 }
460
461 G4double edep = E - finalT;
462
463 aParticleChange.ProposeEnergy(finalT);
464
465 // Deexcitation of ionised atoms
466 std::vector<G4DynamicParticle*>* deexcitatio
467 if (theFluo) deexcitationProducts = Deexcite
468
469 size_t nSecondaries = 0;
470 if (deexcitationProducts != 0) nSecondaries
471 aParticleChange.SetNumberOfSecondaries(nSeco
472
473 if (nSecondaries > 0) {
474
475 const G4StepPoint* preStep = stepData.GetP
476 const G4StepPoint* postStep = stepData.Get
477 G4ThreeVector r = preStep->GetPosition();
478 G4ThreeVector deltaR = postStep->GetPositi
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 = (*deexcitation
489 if (part != 0) {
490 G4double eSecondary = part->GetKinetic
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, ti
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(ed
513
514 return &aParticleChange;
515 }
516
517 //
518
519