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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 18.0 KB
RevLine 
[819]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//
[1192]27// $Id: G4eLowEnergyLoss.cc,v 1.37 2009/07/23 09:15:37 vnivanch Exp $
[1196]28// GEANT4 tag $Name: geant4-09-03-cand-01 $
[819]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
[1192]275 for (G4int i=0; i<=NbinEloss; i++)
[819]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.