source: trunk/source/processes/electromagnetic/lowenergy/src/G4eLowEnergyLoss.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 18.0 KB
Line 
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: G4eLowEnergyLoss.cc,v 1.37 2009/07/23 09:15:37 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-04-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 "G4EnergyLossMessenger.hh"
70#include "G4Poisson.hh"
71#include "G4ProductionCutsTable.hh"
72
73//
74
75// Initialisation of static data members
76// -------------------------------------
77// Contributing processes : ion.loss + soft brems->NbOfProcesses is initialized
78// to 2 . YOU DO NOT HAVE TO CHANGE this variable for a 'normal' run.
79//
80// You have to change NbOfProcesses if you invent a new process contributing
81// to the continuous energy loss.
82// The NbOfProcesses data member can be changed using the (public static)
83// functions Get/Set/Plus/MinusNbOfProcesses (see G4eLowEnergyLoss.hh)
84
85G4int            G4eLowEnergyLoss::NbOfProcesses = 2;
86
87G4int            G4eLowEnergyLoss::CounterOfElectronProcess = 0;
88G4int            G4eLowEnergyLoss::CounterOfPositronProcess = 0;
89G4PhysicsTable** G4eLowEnergyLoss::RecorderOfElectronProcess =
90                                           new G4PhysicsTable*[10];
91G4PhysicsTable** G4eLowEnergyLoss::RecorderOfPositronProcess =
92                                           new G4PhysicsTable*[10];
93
94
95G4PhysicsTable*  G4eLowEnergyLoss::theDEDXElectronTable         = 0;
96G4PhysicsTable*  G4eLowEnergyLoss::theDEDXPositronTable         = 0;
97G4PhysicsTable*  G4eLowEnergyLoss::theRangeElectronTable        = 0;
98G4PhysicsTable*  G4eLowEnergyLoss::theRangePositronTable        = 0;
99G4PhysicsTable*  G4eLowEnergyLoss::theInverseRangeElectronTable = 0;
100G4PhysicsTable*  G4eLowEnergyLoss::theInverseRangePositronTable = 0;
101G4PhysicsTable*  G4eLowEnergyLoss::theLabTimeElectronTable      = 0;
102G4PhysicsTable*  G4eLowEnergyLoss::theLabTimePositronTable      = 0;
103G4PhysicsTable*  G4eLowEnergyLoss::theProperTimeElectronTable   = 0;
104G4PhysicsTable*  G4eLowEnergyLoss::theProperTimePositronTable   = 0;
105
106G4PhysicsTable*  G4eLowEnergyLoss::theeRangeCoeffATable         = 0;
107G4PhysicsTable*  G4eLowEnergyLoss::theeRangeCoeffBTable         = 0;
108G4PhysicsTable*  G4eLowEnergyLoss::theeRangeCoeffCTable         = 0;
109G4PhysicsTable*  G4eLowEnergyLoss::thepRangeCoeffATable         = 0;
110G4PhysicsTable*  G4eLowEnergyLoss::thepRangeCoeffBTable         = 0;
111G4PhysicsTable*  G4eLowEnergyLoss::thepRangeCoeffCTable         = 0;
112
113G4double         G4eLowEnergyLoss::LowerBoundEloss = 10.*eV ;
114G4double         G4eLowEnergyLoss::UpperBoundEloss = 100.*GeV ;
115G4int            G4eLowEnergyLoss::NbinEloss = 360 ;
116G4double         G4eLowEnergyLoss::RTable ;
117G4double         G4eLowEnergyLoss::LOGRTable ;
118
119
120G4EnergyLossMessenger* G4eLowEnergyLoss::eLossMessenger         = 0;
121
122//   
123 
124// constructor and destructor
125 
126G4eLowEnergyLoss::G4eLowEnergyLoss(const G4String& processName)
127   : G4VeLowEnergyLoss (processName),
128     theLossTable(0),
129     MinKineticEnergy(1.*eV),
130     Charge(-1.),
131     lastCharge(0.),
132     theDEDXTable(0),
133     CounterOfProcess(0),
134     RecorderOfProcess(0),
135     fdEdx(0),
136     fRangeNow(0),
137     linLossLimit(0.05),
138     theFluo(false)
139{
140 
141 //create (only once) EnergyLoss messenger
142 if(!eLossMessenger) eLossMessenger = new G4EnergyLossMessenger();
143}
144
145//
146
147G4eLowEnergyLoss::~G4eLowEnergyLoss() 
148{
149     if (theLossTable) 
150       {
151         theLossTable->clearAndDestroy();
152         delete theLossTable;
153       }
154}
155
156void G4eLowEnergyLoss::SetNbOfProcesses(G4int nb) 
157{
158    NbOfProcesses=nb;
159}
160
161void G4eLowEnergyLoss::PlusNbOfProcesses()       
162{
163    NbOfProcesses++;
164}
165
166void G4eLowEnergyLoss::MinusNbOfProcesses() 
167{
168    NbOfProcesses--;
169}                                     
170
171G4int G4eLowEnergyLoss::GetNbOfProcesses() 
172{
173    return NbOfProcesses;
174}
175   
176void G4eLowEnergyLoss::SetLowerBoundEloss(G4double val) 
177{
178    LowerBoundEloss=val;
179}
180   
181void G4eLowEnergyLoss::SetUpperBoundEloss(G4double val) 
182{
183    UpperBoundEloss=val;
184} 
185
186void G4eLowEnergyLoss::SetNbinEloss(G4int nb)
187{
188    NbinEloss=nb;
189}
190 
191G4double G4eLowEnergyLoss::GetLowerBoundEloss()
192{
193    return LowerBoundEloss;
194} 
195   
196G4double G4eLowEnergyLoss::GetUpperBoundEloss() 
197{
198    return UpperBoundEloss;
199} 
200
201G4int G4eLowEnergyLoss::GetNbinEloss() 
202{
203    return NbinEloss;
204} 
205//     
206
207void G4eLowEnergyLoss::BuildDEDXTable(
208                         const G4ParticleDefinition& aParticleType)
209{
210  ParticleMass = aParticleType.GetPDGMass(); 
211  Charge = aParticleType.GetPDGCharge()/eplus;
212
213  //  calculate data members LOGRTable,RTable first
214
215  G4double lrate = std::log(UpperBoundEloss/LowerBoundEloss);
216  LOGRTable=lrate/NbinEloss;
217  RTable   =std::exp(LOGRTable);
218  // Build energy loss table as a sum of the energy loss due to the
219  // different processes.
220  //
221
222  const G4ProductionCutsTable* theCoupleTable=
223        G4ProductionCutsTable::GetProductionCutsTable();
224  size_t numOfCouples = theCoupleTable->GetTableSize();
225
226  // create table for the total energy loss
227
228  if (&aParticleType==G4Electron::Electron())
229    {
230      RecorderOfProcess=RecorderOfElectronProcess;
231      CounterOfProcess=CounterOfElectronProcess;
232      if (CounterOfProcess == NbOfProcesses)
233        {
234         if (theDEDXElectronTable)
235           {
236             theDEDXElectronTable->clearAndDestroy();
237             delete theDEDXElectronTable;
238           }
239         theDEDXElectronTable = new G4PhysicsTable(numOfCouples);
240         theDEDXTable = theDEDXElectronTable;
241        }
242    }
243  if (&aParticleType==G4Positron::Positron())
244    {
245     RecorderOfProcess=RecorderOfPositronProcess;
246     CounterOfProcess=CounterOfPositronProcess;
247     if (CounterOfProcess == NbOfProcesses)
248       {
249        if (theDEDXPositronTable)
250          {
251            theDEDXPositronTable->clearAndDestroy();
252            delete theDEDXPositronTable;
253          }
254        theDEDXPositronTable = new G4PhysicsTable(numOfCouples);
255        theDEDXTable = theDEDXPositronTable;
256       }
257    }
258
259  if (CounterOfProcess == NbOfProcesses)
260    {
261     // fill the tables
262     // loop for materials
263     G4double LowEdgeEnergy , Value;
264     G4bool isOutRange;
265     G4PhysicsTable* pointer;
266
267     for (size_t J=0; J<numOfCouples; J++)
268        {
269         // create physics vector and fill it
270
271         G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
272                    LowerBoundEloss, UpperBoundEloss, NbinEloss);
273
274         // loop for the kinetic energy
275         for (G4int i=0; i<=NbinEloss; i++)
276            {
277              LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
278              //here comes the sum of the different tables created by the
279              //processes (ionisation,bremsstrahlung,etc...)
280              Value = 0.;
281              for (G4int process=0; process < NbOfProcesses; process++)
282                 {
283                   pointer= RecorderOfProcess[process];
284                   Value += (*pointer)[J]->GetValue(LowEdgeEnergy,isOutRange);
285                 }
286
287              aVector->PutValue(i,Value) ;
288            }
289
290         theDEDXTable->insert(aVector) ;
291
292        }
293
294
295     //reset counter to zero
296     if (&aParticleType==G4Electron::Electron()) CounterOfElectronProcess=0;
297     if (&aParticleType==G4Positron::Positron()) CounterOfPositronProcess=0;
298
299     ParticleMass = aParticleType.GetPDGMass();
300
301     if (&aParticleType==G4Electron::Electron())
302     {
303       // Build range table
304       theRangeElectronTable = BuildRangeTable(theDEDXElectronTable,
305                                               theRangeElectronTable,
306                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
307
308       // Build lab/proper time tables
309       theLabTimeElectronTable = BuildLabTimeTable(theDEDXElectronTable,
310                         theLabTimeElectronTable,
311                         LowerBoundEloss,UpperBoundEloss,NbinEloss);
312       theProperTimeElectronTable = BuildProperTimeTable(theDEDXElectronTable,
313                            theProperTimeElectronTable,
314                            LowerBoundEloss,UpperBoundEloss,NbinEloss);
315
316       // Build coeff tables for the energy loss calculation
317       theeRangeCoeffATable = BuildRangeCoeffATable(theRangeElectronTable,
318                             theeRangeCoeffATable,
319                             LowerBoundEloss,UpperBoundEloss,NbinEloss);
320
321       theeRangeCoeffBTable = BuildRangeCoeffBTable(theRangeElectronTable,
322                             theeRangeCoeffBTable,
323                             LowerBoundEloss,UpperBoundEloss,NbinEloss);
324
325       theeRangeCoeffCTable = BuildRangeCoeffCTable(theRangeElectronTable,
326                             theeRangeCoeffCTable,
327                             LowerBoundEloss,UpperBoundEloss,NbinEloss);
328
329       // invert the range table
330       theInverseRangeElectronTable = BuildInverseRangeTable(theRangeElectronTable,
331                              theeRangeCoeffATable,
332                              theeRangeCoeffBTable,
333                              theeRangeCoeffCTable,
334                              theInverseRangeElectronTable,
335                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
336     }
337     if (&aParticleType==G4Positron::Positron())
338     {
339       // Build range table
340       theRangePositronTable = BuildRangeTable(theDEDXPositronTable,
341                                               theRangePositronTable,
342                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
343
344
345       // Build lab/proper time tables
346       theLabTimePositronTable = BuildLabTimeTable(theDEDXPositronTable,
347                         theLabTimePositronTable,
348                         LowerBoundEloss,UpperBoundEloss,NbinEloss);
349       theProperTimePositronTable = BuildProperTimeTable(theDEDXPositronTable,
350                            theProperTimePositronTable,
351                            LowerBoundEloss,UpperBoundEloss,NbinEloss);
352
353       // Build coeff tables for the energy loss calculation
354       thepRangeCoeffATable = BuildRangeCoeffATable(theRangePositronTable,
355                             thepRangeCoeffATable,
356                             LowerBoundEloss,UpperBoundEloss,NbinEloss);
357
358       thepRangeCoeffBTable = BuildRangeCoeffBTable(theRangePositronTable,
359                             thepRangeCoeffBTable,
360                             LowerBoundEloss,UpperBoundEloss,NbinEloss);
361
362       thepRangeCoeffCTable = BuildRangeCoeffCTable(theRangePositronTable,
363                             thepRangeCoeffCTable,
364                             LowerBoundEloss,UpperBoundEloss,NbinEloss);
365
366       // invert the range table
367       theInverseRangePositronTable = BuildInverseRangeTable(theRangePositronTable,
368                              thepRangeCoeffATable,
369                              thepRangeCoeffBTable,
370                              thepRangeCoeffCTable,
371                              theInverseRangePositronTable,
372                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
373     }
374
375     if(verboseLevel > 1) {
376       G4cout << (*theDEDXElectronTable) << G4endl;
377     }
378
379
380     // make the energy loss and the range table available
381     G4EnergyLossTables::Register(&aParticleType,
382       (&aParticleType==G4Electron::Electron())?
383       theDEDXElectronTable: theDEDXPositronTable,
384       (&aParticleType==G4Electron::Electron())?
385       theRangeElectronTable: theRangePositronTable,
386       (&aParticleType==G4Electron::Electron())?
387       theInverseRangeElectronTable: theInverseRangePositronTable,
388       (&aParticleType==G4Electron::Electron())?
389       theLabTimeElectronTable: theLabTimePositronTable,
390       (&aParticleType==G4Electron::Electron())?
391       theProperTimeElectronTable: theProperTimePositronTable,
392       LowerBoundEloss, UpperBoundEloss, 1.,NbinEloss);
393    }
394}
395
396//
397
398G4VParticleChange* G4eLowEnergyLoss::AlongStepDoIt( const G4Track& trackData,
399                                                 const G4Step&  stepData)
400{
401 // compute the energy loss after a Step
402
403  static const G4double faclow = 1.5 ;
404
405  // get particle and material pointers from trackData
406  const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
407  G4double E      = aParticle->GetKineticEnergy();
408
409  const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
410
411  G4double Step = stepData.GetStepLength();
412
413  aParticleChange.Initialize(trackData);
414  //fParticleChange.Initialize(trackData);
415
416  G4double MeanLoss, finalT;
417
418  if (E < MinKineticEnergy)   finalT = 0.;
419
420  else if ( E< faclow*LowerBoundEloss)
421  {
422    if (Step >= fRangeNow)  finalT = 0.;
423   //  else finalT = E*(1.-Step/fRangeNow) ;
424    else finalT = E*(1.-std::sqrt(Step/fRangeNow)) ;
425  }
426
427  else if (E>=UpperBoundEloss) finalT = E - Step*fdEdx;
428
429  else if (Step >= fRangeNow)  finalT = 0.;
430
431  else
432  {
433    if(Step/fRangeNow < linLossLimit) finalT = E-Step*fdEdx ;
434    else
435    {
436      if (Charge<0.) finalT = G4EnergyLossTables::GetPreciseEnergyFromRange
437                             (G4Electron::Electron(),fRangeNow-Step,couple);
438      else           finalT = G4EnergyLossTables::GetPreciseEnergyFromRange
439                             (G4Positron::Positron(),fRangeNow-Step,couple);
440     }
441  }
442
443  if(finalT < MinKineticEnergy) finalT = 0. ;
444
445  MeanLoss = E-finalT ;
446
447  //now the loss with fluctuation
448  if ((EnlossFlucFlag) && (finalT > 0.) && (finalT < E)&&(E > LowerBoundEloss))
449  {
450    finalT = E-GetLossWithFluct(aParticle,couple,MeanLoss,Step);
451    if (finalT < 0.) finalT = 0.;
452  }
453
454  // kill the particle if the kinetic energy <= 0
455  if (finalT <= 0. )
456  {
457    finalT = 0.;
458    if(Charge > 0.0) aParticleChange.ProposeTrackStatus(fStopButAlive);
459    else             aParticleChange.ProposeTrackStatus(fStopAndKill);
460  }
461
462  G4double edep = E - finalT;
463
464  aParticleChange.ProposeEnergy(finalT);
465
466  // Deexcitation of ionised atoms
467  std::vector<G4DynamicParticle*>* deexcitationProducts = 0;
468  if (theFluo) deexcitationProducts = DeexciteAtom(couple,E,edep);
469
470  size_t nSecondaries = 0;
471  if (deexcitationProducts != 0) nSecondaries = deexcitationProducts->size();
472  aParticleChange.SetNumberOfSecondaries(nSecondaries);
473
474  if (nSecondaries > 0) {
475
476    const G4StepPoint* preStep = stepData.GetPreStepPoint();
477    const G4StepPoint* postStep = stepData.GetPostStepPoint();
478    G4ThreeVector r = preStep->GetPosition();
479    G4ThreeVector deltaR = postStep->GetPosition();
480    deltaR -= r;
481    G4double t = preStep->GetGlobalTime();
482    G4double deltaT = postStep->GetGlobalTime();
483    deltaT -= t;
484    G4double time, q;
485    G4ThreeVector position;
486
487    for (size_t i=0; i<nSecondaries; i++) {
488
489      G4DynamicParticle* part = (*deexcitationProducts)[i];
490      if (part != 0) {
491        G4double eSecondary = part->GetKineticEnergy();
492        edep -= eSecondary;
493        if (edep > 0.)
494          {
495            q = G4UniformRand();
496            time = deltaT*q + t;
497            position  = deltaR*q;
498            position += r;
499            G4Track* newTrack = new G4Track(part, time, position);
500            aParticleChange.AddSecondary(newTrack);
501          }
502        else
503          {
504            edep += eSecondary;
505            delete part;
506            part = 0;
507          }
508      }
509    }
510  }
511  delete deexcitationProducts;
512
513  aParticleChange.ProposeLocalEnergyDeposit(edep);
514
515  return &aParticleChange;
516}
517
518//
519
Note: See TracBrowser for help on using the repository browser.