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