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