ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eLowEnergyLoss.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eLowEnergyLoss.cc
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 "G4EnergyLossMessenger.hh"
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 
85 
89  new G4PhysicsTable*[10];
91  new G4PhysicsTable*[10];
92 
93 
104 
111 
117 
118 
119 G4EnergyLossMessenger* G4eLowEnergyLoss::eLossMessenger = 0;
120 
121 //
122 
123 // constructor and destructor
124 
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 
147 {
148  if (theLossTable)
149  {
151  delete theLossTable;
152  }
153 }
154 
156 {
157  NbOfProcesses=nb;
158 }
159 
161 {
162  NbOfProcesses++;
163 }
164 
166 {
167  NbOfProcesses--;
168 }
169 
171 {
172  return NbOfProcesses;
173 }
174 
176 {
177  LowerBoundEloss=val;
178 }
179 
181 {
182  UpperBoundEloss=val;
183 }
184 
186 {
187  NbinEloss=nb;
188 }
189 
191 {
192  return LowerBoundEloss;
193 }
194 
196 {
197  return UpperBoundEloss;
198 }
199 
201 {
202  return NbinEloss;
203 }
204 //
205 
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=
223  size_t numOfCouples = theCoupleTable->GetTableSize();
224 
225  // create table for the total energy loss
226 
227  if (&aParticleType==G4Electron::Electron())
228  {
232  {
234  {
236  delete theDEDXElectronTable;
237  }
238  theDEDXElectronTable = new G4PhysicsTable(numOfCouples);
240  }
241  }
242  if (&aParticleType==G4Positron::Positron())
243  {
247  {
249  {
251  delete theDEDXPositronTable;
252  }
253  theDEDXPositronTable = new G4PhysicsTable(numOfCouples);
255  }
256  }
257 
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(
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
306 
307  // Build lab/proper time tables
314 
315  // Build coeff tables for the energy loss calculation
319 
323 
327 
328  // invert the range table
335  }
336  if (&aParticleType==G4Positron::Positron())
337  {
338  // Build range table
342 
343 
344  // Build lab/proper time tables
351 
352  // Build coeff tables for the energy loss calculation
356 
360 
364 
365  // invert the range table
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())?
383  (&aParticleType==G4Electron::Electron())?
385  (&aParticleType==G4Electron::Electron())?
387  (&aParticleType==G4Electron::Electron())?
389  (&aParticleType==G4Electron::Electron())?
392  }
393 }
394 
395 //
396 
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  {
436  (G4Electron::Electron(),fRangeNow-Step,couple);
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.;
459  }
460 
461  G4double edep = E - finalT;
462 
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();
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;
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 
513 
514  return &aParticleChange;
515 }
516 
517 //
518