source: trunk/source/processes/electromagnetic/utils/include/G4VEmProcess.hh @ 1230

Last change on this file since 1230 was 1228, checked in by garnier, 15 years ago

update geant4.9.3 tag
File size: 20.2 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// $Id: G4VEmProcess.hh,v 1.55 2009/09/23 14:42:47 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class header file
32//
33//
34// File name:     G4VEmProcess
35//
36// Author:        Vladimir Ivanchenko
37//
38// Creation date: 01.10.2003
39//
40// Modifications:
41// 30-06-04 make destructor virtual (V.Ivanchenko)
42// 09-08-04 optimise integral option (V.Ivanchenko)
43// 11-08-04 add protected methods to access cuts (V.Ivanchenko)
44// 09-09-04 Bug fix for the integral mode with 2 peaks (V.Ivanchneko)
45// 16-09-04 Add flag for LambdaTable and method RecalculateLambda (VI)
46// 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko)
47// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
48// 18-04-05 Use G4ParticleChangeForGamma (V.Ivantchenko)
49// 09-05-05 Fix problem in logic when path boundary between materials (VI)
50// 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI)
51// 01-02-06 put default value A=0. to keep compatibility with v5.2 (mma)
52// 13-05-06 Add method to access model by index (V.Ivanchenko)
53// 12-09-06 add SetModel() (mma)
54// 25-09-07 More accurate handling zero xsect in
55//          PostStepGetPhysicalInteractionLength (V.Ivanchenko)
56// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
57// 15-07-08 Reorder class members for further multi-thread development (VI)
58//
59// Class Description:
60//
61// It is the unified Discrete process
62
63// -------------------------------------------------------------------
64//
65
66#ifndef G4VEmProcess_h
67#define G4VEmProcess_h 1
68
69#include "G4VDiscreteProcess.hh"
70#include "globals.hh"
71#include "G4Material.hh"
72#include "G4MaterialCutsCouple.hh"
73#include "G4Track.hh"
74#include "G4EmModelManager.hh"
75#include "G4UnitsTable.hh"
76#include "G4ParticleDefinition.hh"
77#include "G4ParticleChangeForGamma.hh"
78
79class G4Step;
80class G4VEmModel;
81class G4DataVector;
82class G4VParticleChange;
83class G4PhysicsTable;
84class G4PhysicsVector;
85
86//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87
88class G4VEmProcess : public G4VDiscreteProcess
89{
90public:
91
92  G4VEmProcess(const G4String& name,
93               G4ProcessType type = fElectromagnetic);
94
95  virtual ~G4VEmProcess();
96
97  //------------------------------------------------------------------------
98  // Virtual methods to be implemented in concrete processes
99  //------------------------------------------------------------------------
100
101  virtual G4bool IsApplicable(const G4ParticleDefinition& p) = 0;
102
103  virtual void PrintInfo() = 0;
104
105protected:
106
107  virtual void InitialiseProcess(const G4ParticleDefinition*) = 0;
108
109  //------------------------------------------------------------------------
110  // Implementation of virtual methods common to all Discrete processes
111  //------------------------------------------------------------------------
112
113public:
114
115  // Initialise for build of tables
116  void PreparePhysicsTable(const G4ParticleDefinition&);
117
118  // Build physics table during initialisation
119  void BuildPhysicsTable(const G4ParticleDefinition&);
120
121  void PrintInfoDefinition();
122
123  // implementation of virtual method, specific for G4VEmProcess
124  G4double PostStepGetPhysicalInteractionLength(
125                             const G4Track& track,
126                             G4double   previousStepSize,
127                             G4ForceCondition* condition
128                            );
129
130  // implementation of virtual method, specific for G4VEmProcess
131  G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&);
132
133  // Store PhysicsTable in a file.
134  // Return false in case of failure at I/O
135  G4bool StorePhysicsTable(const G4ParticleDefinition*,
136                           const G4String& directory,
137                           G4bool ascii = false);
138
139  // Retrieve Physics from a file.
140  // (return true if the Physics Table can be build by using file)
141  // (return false if the process has no functionality or in case of failure)
142  // File name should is constructed as processName+particleName and the
143  // should be placed under the directory specifed by the argument.
144  G4bool RetrievePhysicsTable(const G4ParticleDefinition*,
145                              const G4String& directory,
146                              G4bool ascii);
147
148  // deexcitation activated per G4Region
149  void ActivateDeexcitation(G4bool, const G4Region* r = 0);
150
151  //------------------------------------------------------------------------
152  // Specific methods for Discrete EM post step simulation
153  //------------------------------------------------------------------------
154
155  // It returns the cross section per volume for energy/ material
156  G4double CrossSectionPerVolume(G4double kineticEnergy,
157                                 const G4MaterialCutsCouple* couple);
158
159  // It returns the cross section of the process per atom
160  inline G4double ComputeCrossSectionPerAtom(G4double kineticEnergy, 
161                                             G4double Z, G4double A=0., 
162                                             G4double cut=0.0);
163
164  inline G4double MeanFreePath(const G4Track& track);
165
166  // It returns cross section per volume
167  inline G4double GetLambda(G4double& kinEnergy, 
168                            const G4MaterialCutsCouple* couple);
169
170  //------------------------------------------------------------------------
171  // Specific methods to build and access Physics Tables
172  //------------------------------------------------------------------------
173
174  // Binning for lambda table
175  inline void SetLambdaBinning(G4int nbins);
176  inline G4int LambdaBinning() const;
177
178  // Min kinetic energy for tables
179  inline void SetMinKinEnergy(G4double e);
180  inline G4double MinKinEnergy() const;
181
182  // Max kinetic energy for tables
183  inline void SetMaxKinEnergy(G4double e);
184  inline G4double MaxKinEnergy() const;
185
186  inline void SetPolarAngleLimit(G4double a);
187  inline G4double PolarAngleLimit() const;
188
189  inline const G4PhysicsTable* LambdaTable() const;
190
191  //------------------------------------------------------------------------
192  // Define and access particle type
193  //------------------------------------------------------------------------
194
195  inline const G4ParticleDefinition* Particle() const;
196  inline const G4ParticleDefinition* SecondaryParticle() const;
197
198  //------------------------------------------------------------------------
199  // Specific methods to set, access, modify models and basic parameters
200  //------------------------------------------------------------------------
201
202protected:
203  // Select model in run time
204  inline G4VEmModel* SelectModel(G4double& kinEnergy, size_t index); 
205
206public:
207  // Select model by energy and region index
208  inline G4VEmModel* SelectModelForMaterial(G4double kinEnergy, 
209                                            size_t& idxRegion) const;
210   
211  // Add model for region, smaller value of order defines which
212  // model will be selected for a given energy interval 
213  void AddEmModel(G4int, G4VEmModel*, const G4Region* region = 0);
214
215  // Assign a model to a process
216  void SetModel(G4VEmModel*, G4int index = 1);
217 
218  // return the assigned model
219  G4VEmModel* Model(G4int index = 1);
220   
221  // Define new energy range for the model identified by the name
222  void UpdateEmModel(const G4String&, G4double, G4double);
223
224  // Access to models
225  G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
226
227  inline void SetLambdaFactor(G4double val);
228
229  inline void SetIntegral(G4bool val);
230  inline G4bool IsIntegral() const;
231
232  inline void SetApplyCuts(G4bool val);
233
234  //------------------------------------------------------------------------
235  // Other generic methods
236  //------------------------------------------------------------------------
237 
238protected:
239
240  G4double GetMeanFreePath(const G4Track& track,
241                           G4double previousStepSize,
242                           G4ForceCondition* condition);
243
244  G4PhysicsVector* LambdaPhysicsVector(const G4MaterialCutsCouple*);
245
246  inline G4double RecalculateLambda(G4double kinEnergy,
247                                    const G4MaterialCutsCouple* couple);
248
249  inline G4ParticleChangeForGamma* GetParticleChange();
250
251  inline void SetParticle(const G4ParticleDefinition* p);
252 
253  inline void SetSecondaryParticle(const G4ParticleDefinition* p);
254
255  inline size_t CurrentMaterialCutsCoupleIndex() const;
256
257  inline G4double GetGammaEnergyCut();
258
259  inline G4double GetElectronEnergyCut();
260
261  inline void SetBuildTableFlag(G4bool val);
262
263  inline void SetStartFromNullFlag(G4bool val);
264
265private:
266
267  void Clear();
268
269  void BuildLambdaTable();
270
271  void FindLambdaMax();
272
273  inline void InitialiseStep(const G4Track&);
274
275  inline void DefineMaterial(const G4MaterialCutsCouple* couple);
276
277  inline void ComputeIntegralLambda(G4double kinEnergy);
278
279  inline G4double GetLambdaFromTable(G4double kinEnergy);
280
281  inline G4double GetCurrentLambda(G4double kinEnergy);
282
283  inline G4double ComputeCurrentLambda(G4double kinEnergy);
284
285  // copy constructor and hide assignment operator
286  G4VEmProcess(G4VEmProcess &);
287  G4VEmProcess & operator=(const G4VEmProcess &right);
288
289  // ======== Parameters of the class fixed at construction =========
290
291  G4EmModelManager*            modelManager;
292  const G4ParticleDefinition*  theGamma;
293  const G4ParticleDefinition*  theElectron;
294  const G4ParticleDefinition*  thePositron;
295  const G4ParticleDefinition*  secondaryParticle;
296
297  G4bool                       buildLambdaTable;
298
299  // ======== Parameters of the class fixed at initialisation =======
300
301  std::vector<G4VEmModel*>     emModels;
302
303  // tables and vectors
304  G4PhysicsTable*              theLambdaTable;
305  G4double*                    theEnergyOfCrossSectionMax;
306  G4double*                    theCrossSectionMax;
307
308  const std::vector<G4double>* theCuts;
309  const std::vector<G4double>* theCutsGamma;
310  const std::vector<G4double>* theCutsElectron;
311  const std::vector<G4double>* theCutsPositron;
312
313  G4int                        nLambdaBins;
314
315  G4double                     minKinEnergy;
316  G4double                     maxKinEnergy;
317  G4double                     lambdaFactor;
318  G4double                     polarAngleLimit;
319
320  G4bool                       integral;
321  G4bool                       applyCuts;
322  G4bool                       startFromNull;
323  G4bool                       useDeexcitation;
324
325  G4int                        nDERegions;
326  std::vector<const G4Region*> deRegions;
327  G4bool*                      idxDERegions;
328
329  // ======== Cashed values - may be state dependent ================
330
331protected:
332
333  G4ParticleChangeForGamma     fParticleChange;
334
335private:
336
337  std::vector<G4DynamicParticle*> secParticles;
338
339  G4VEmModel*                  currentModel; 
340
341  const G4ParticleDefinition*  particle;
342
343  // cash
344  const G4Material*            currentMaterial;
345  const G4MaterialCutsCouple*  currentCouple;
346  size_t                       currentCoupleIndex;
347
348  G4double                     mfpKinEnergy;
349  G4double                     preStepKinEnergy;
350  G4double                     preStepLambda;
351
352};
353
354//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
356
357inline G4double G4VEmProcess::ComputeCrossSectionPerAtom(
358                 G4double kineticEnergy, G4double Z, G4double A, G4double cut)
359{
360  SelectModel(kineticEnergy, currentCoupleIndex);
361  G4double x = 0.0;
362  if(currentModel) {
363   x = currentModel->ComputeCrossSectionPerAtom(particle,kineticEnergy,
364                                                 Z,A,cut);
365  }
366  return x;
367}
368
369//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
370
371inline G4double G4VEmProcess::MeanFreePath(const G4Track& track)
372{
373  DefineMaterial(track.GetMaterialCutsCouple());
374  preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
375  G4double x = DBL_MAX;
376  if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
377  return x;
378}
379
380//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381
382inline G4double G4VEmProcess::GetLambda(G4double& kineticEnergy,
383                                        const G4MaterialCutsCouple* couple)
384{
385  DefineMaterial(couple);
386  return GetCurrentLambda(kineticEnergy);
387}
388
389//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390
391inline void G4VEmProcess::SetLambdaBinning(G4int nbins)
392{
393  nLambdaBins = nbins;
394}
395
396//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
397
398inline G4int G4VEmProcess::LambdaBinning() const
399{
400  return nLambdaBins;
401}
402
403//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
404
405inline void G4VEmProcess::SetMinKinEnergy(G4double e)
406{
407  minKinEnergy = e;
408}
409
410//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
411
412inline G4double G4VEmProcess::MinKinEnergy() const
413{
414  return minKinEnergy;
415}
416
417//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
418
419inline void G4VEmProcess::SetMaxKinEnergy(G4double e)
420{
421  maxKinEnergy = e;
422}
423
424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
425
426inline G4double G4VEmProcess::MaxKinEnergy() const
427{
428  return maxKinEnergy;
429}
430
431//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
432
433inline void G4VEmProcess::SetPolarAngleLimit(G4double val)
434{
435  if(val < 0.0)     polarAngleLimit = 0.0;
436  else if(val > pi) polarAngleLimit = pi;
437  else              polarAngleLimit = val;
438}
439
440//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
441
442inline G4double G4VEmProcess::PolarAngleLimit() const
443{
444  return polarAngleLimit;
445}
446
447//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
448
449inline const G4PhysicsTable* G4VEmProcess::LambdaTable() const
450{
451  return theLambdaTable;
452}
453
454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
455
456inline const G4ParticleDefinition* G4VEmProcess::Particle() const
457{
458  return particle;
459}
460
461//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
462
463inline const G4ParticleDefinition* G4VEmProcess::SecondaryParticle() const
464{
465  return secondaryParticle;
466}
467
468//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
469
470inline 
471G4VEmModel* G4VEmProcess::SelectModel(G4double& kinEnergy, size_t index)
472{
473  currentModel = modelManager->SelectModel(kinEnergy, index);
474  currentModel->SetCurrentCouple(currentCouple);
475  return currentModel;
476}
477
478//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
479
480inline 
481G4VEmModel* G4VEmProcess::SelectModelForMaterial(G4double kinEnergy, 
482                                                 size_t& idxRegion) const
483{
484  return modelManager->SelectModel(kinEnergy, idxRegion);
485}
486
487//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
488
489inline void G4VEmProcess::SetLambdaFactor(G4double val)
490{
491  if(val > 0.0 && val <= 1.0) lambdaFactor = val;
492}
493
494//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
495
496inline void G4VEmProcess::SetIntegral(G4bool val)
497{
498  if(particle && particle != theGamma) integral = val;
499  if(integral) buildLambdaTable = true;
500}
501
502//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
503
504inline G4bool G4VEmProcess::IsIntegral() const
505{
506  return integral;
507}
508
509//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
510
511inline void G4VEmProcess::SetApplyCuts(G4bool val)
512{
513  applyCuts = val;
514}
515
516//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
517
518inline G4double G4VEmProcess::RecalculateLambda(G4double e, 
519                                                const G4MaterialCutsCouple* couple)
520{
521  DefineMaterial(couple);
522  return ComputeCurrentLambda(e);
523}
524
525//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
526
527inline G4ParticleChangeForGamma* G4VEmProcess::GetParticleChange()
528{
529  return &fParticleChange;
530}
531
532//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
533
534inline void G4VEmProcess::SetParticle(const G4ParticleDefinition* p)
535{
536  particle = p;
537}
538
539//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
540
541inline void G4VEmProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
542{
543  secondaryParticle = p;
544}
545
546//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
547
548inline size_t G4VEmProcess::CurrentMaterialCutsCoupleIndex() const 
549{
550  return currentCoupleIndex;
551}
552
553//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
554
555inline G4double G4VEmProcess::GetGammaEnergyCut()
556{
557  return (*theCutsGamma)[currentCoupleIndex];
558}
559
560//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
561
562inline G4double G4VEmProcess::GetElectronEnergyCut()
563{
564  return (*theCutsElectron)[currentCoupleIndex];
565}
566
567//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
568
569inline void G4VEmProcess::SetBuildTableFlag(G4bool val)
570{
571  buildLambdaTable = val;
572  if(!val) integral = false;
573}
574
575//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
576
577inline void G4VEmProcess::SetStartFromNullFlag(G4bool val)
578{
579  startFromNull = val;
580}
581
582//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
583
584inline void G4VEmProcess::InitialiseStep(const G4Track& track)
585{
586  preStepKinEnergy = track.GetKineticEnergy();
587  DefineMaterial(track.GetMaterialCutsCouple());
588  SelectModel(preStepKinEnergy, currentCoupleIndex);
589  if (theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
590}
591
592//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
593
594inline void G4VEmProcess::DefineMaterial(const G4MaterialCutsCouple* couple)
595{
596  if(couple != currentCouple) {
597    currentCouple   = couple;
598    currentMaterial = couple->GetMaterial();
599    currentCoupleIndex = couple->GetIndex();
600    mfpKinEnergy = DBL_MAX;
601  }
602}
603
604//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
605
606inline void G4VEmProcess::ComputeIntegralLambda(G4double e)
607{
608  mfpKinEnergy  = theEnergyOfCrossSectionMax[currentCoupleIndex];
609  if (e <= mfpKinEnergy) {
610    preStepLambda = GetLambdaFromTable(e);
611
612  } else {
613    G4double e1 = e*lambdaFactor;
614    if(e1 > mfpKinEnergy) {
615      preStepLambda  = GetLambdaFromTable(e);
616      G4double preStepLambda1 = GetLambdaFromTable(e1);
617      if(preStepLambda1 > preStepLambda) {
618        mfpKinEnergy = e1;
619        preStepLambda = preStepLambda1;
620      }
621    } else {
622      preStepLambda = theCrossSectionMax[currentCoupleIndex];
623    }
624  }
625}
626
627//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
628
629inline G4double G4VEmProcess::GetLambdaFromTable(G4double e)
630{
631  return (((*theLambdaTable)[currentCoupleIndex])->Value(e));
632}
633
634//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
635
636inline G4double G4VEmProcess::GetCurrentLambda(G4double e)
637{
638  G4double x = 0.0;
639  if(theLambdaTable) { x = GetLambdaFromTable(e); }
640  else               { x = ComputeCurrentLambda(e); }
641  return x;
642}
643
644//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
645
646inline G4double G4VEmProcess::ComputeCurrentLambda(G4double e)
647{
648  SelectModel(e, currentCoupleIndex);
649  return currentModel->CrossSectionPerVolume(currentMaterial,particle,
650                                             e,(*theCuts)[currentCoupleIndex]);
651}
652
653//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
654
655#endif
Note: See TracBrowser for help on using the repository browser.