Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VEmModel.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: G4VEmModel
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 03.01.2002
38 //
39 // Modifications:
40 //
41 // 23-12-02 V.Ivanchenko change interface before move to cut per region
42 // 24-01-03 Cut per region (V.Ivanchenko)
43 // 13-02-03 Add name (V.Ivanchenko)
44 // 25-02-03 Add sample theta and displacement (V.Ivanchenko)
45 // 23-07-03 Replace G4Material by G4MaterialCutCouple in dE/dx and CrossSection
46 // calculation (V.Ivanchenko)
47 // 01-03-04 L.Urban signature changed in SampleCosineTheta
48 // 23-04-04 L.urban signature of SampleCosineTheta changed back
49 // 17-11-04 Add method CrossSectionPerAtom (V.Ivanchenko)
50 // 14-03-05 Reduce number of pure virtual methods and make inline part
51 // separate (V.Ivanchenko)
52 // 24-03-05 Remove IsInCharge and add G4VParticleChange in the constructor (VI)
53 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
54 // 15-04-05 optimize internal interface for msc (V.Ivanchenko)
55 // 08-05-05 A -> N (V.Ivanchenko)
56 // 25-07-05 Move constructor and destructor to the body (V.Ivanchenko)
57 // 02-02-06 ComputeCrossSectionPerAtom: default value A=0. (mma)
58 // 06-02-06 add method ComputeMeanFreePath() (mma)
59 // 07-03-06 Optimize msc methods (V.Ivanchenko)
60 // 29-06-06 Add member currentElement and Get/Set methods (V.Ivanchenko)
61 // 29-10-07 Added SampleScattering (V.Ivanchenko)
62 // 15-07-08 Reorder class members and improve comments (VI)
63 // 21-07-08 Added vector of G4ElementSelector and methods to use it (VI)
64 // 12-09-08 Added methods GetParticleCharge, GetChargeSquareRatio,
65 // CorrectionsAlongStep, ActivateNuclearStopping (VI)
66 // 16-02-09 Moved implementations of virtual methods to source (VI)
67 // 07-04-09 Moved msc methods from G4VEmModel to G4VMscModel (VI)
68 // 13-10-10 Added G4VEmAngularDistribution (VI)
69 //
70 // Class Description:
71 //
72 // Abstract interface to energy loss models
73 
74 // -------------------------------------------------------------------
75 //
76 
77 #ifndef G4VEmModel_h
78 #define G4VEmModel_h 1
79 
80 #include "globals.hh"
81 #include "G4DynamicParticle.hh"
82 #include "G4ParticleDefinition.hh"
83 #include "G4MaterialCutsCouple.hh"
84 #include "G4Material.hh"
85 #include "G4Element.hh"
86 #include "G4ElementVector.hh"
87 #include "G4DataVector.hh"
88 #include "G4VEmFluctuationModel.hh"
90 #include "G4EmElementSelector.hh"
91 #include "Randomize.hh"
92 #include <vector>
93 
94 class G4PhysicsTable;
95 class G4Region;
96 class G4VParticleChange;
99 class G4Track;
100 
102 {
103 
104 public:
105 
106  G4VEmModel(const G4String& nam);
107 
108  virtual ~G4VEmModel();
109 
110  //------------------------------------------------------------------------
111  // Virtual methods to be implemented for any concrete model
112  //------------------------------------------------------------------------
113 
114  virtual void Initialise(const G4ParticleDefinition*,
115  const G4DataVector&) = 0;
116 
117  virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
118  const G4MaterialCutsCouple*,
119  const G4DynamicParticle*,
120  G4double tmin = 0.0,
121  G4double tmax = DBL_MAX) = 0;
122 
123  //------------------------------------------------------------------------
124  // Methods with standard implementation; may be overwritten if needed
125  //------------------------------------------------------------------------
126 
127  // main method to compute dEdx
129  const G4ParticleDefinition*,
130  G4double kineticEnergy,
131  G4double cutEnergy = DBL_MAX);
132 
133  // main method to compute cross section per Volume
135  const G4ParticleDefinition*,
136  G4double kineticEnergy,
137  G4double cutEnergy = 0.0,
138  G4double maxEnergy = DBL_MAX);
139 
140  // main method to compute cross section per atom
142  G4double kinEnergy,
143  G4double Z,
144  G4double A = 0., /* amu */
145  G4double cutEnergy = 0.0,
146  G4double maxEnergy = DBL_MAX);
147 
148  // Compute effective ion charge square
149  virtual G4double ChargeSquareRatio(const G4Track&);
150 
151  // Compute effective ion charge square
153  const G4Material*,
154  G4double kineticEnergy);
155 
156  // Compute ion charge
158  const G4Material*,
159  G4double kineticEnergy);
160 
161  // Initialisation for a new track
162  virtual void StartTracking(G4Track*);
163 
164  // add correction to energy loss and compute non-ionizing energy loss
165  virtual void CorrectionsAlongStep(const G4MaterialCutsCouple*,
166  const G4DynamicParticle*,
167  G4double& eloss,
168  G4double& niel,
169  G4double length);
170 
171  // value which may be tabulated (by default cross section)
172  virtual G4double Value(const G4MaterialCutsCouple*,
173  const G4ParticleDefinition*,
174  G4double kineticEnergy);
175 
176  // threshold for zero value
177  virtual G4double MinPrimaryEnergy(const G4Material*,
178  const G4ParticleDefinition*);
179 
180  // initilisation at run time for a given material
181  virtual void SetupForMaterial(const G4ParticleDefinition*,
182  const G4Material*,
183  G4double kineticEnergy);
184 
185  // add a region for the model
186  virtual void DefineForRegion(const G4Region*);
187 
188 protected:
189 
190  // initialisation of the ParticleChange for the model
192 
193  // initialisation of the ParticleChange for the model
195 
196  // kinematically allowed max kinetic energy of a secondary
198  G4double kineticEnergy);
199 
200 public:
201 
202  //------------------------------------------------------------------------
203  // Generic methods common to all models
204  //------------------------------------------------------------------------
205 
206  // should be called at initialisation to build element selectors
208  const G4DataVector&);
209 
210  // dEdx per unit length
212  const G4ParticleDefinition*,
213  G4double kineticEnergy,
214  G4double cutEnergy = DBL_MAX);
215 
216  // cross section per volume
218  const G4ParticleDefinition*,
219  G4double kineticEnergy,
220  G4double cutEnergy = 0.0,
221  G4double maxEnergy = DBL_MAX);
222 
223  // compute mean free path via cross section per volume
225  G4double kineticEnergy,
226  const G4Material*,
227  G4double cutEnergy = 0.0,
228  G4double maxEnergy = DBL_MAX);
229 
230  // generic cross section per element
232  const G4Element*,
233  G4double kinEnergy,
234  G4double cutEnergy = 0.0,
235  G4double maxEnergy = DBL_MAX);
236 
237  // select isotope in order to have precise mass of the nucleus
238  inline G4int SelectIsotopeNumber(const G4Element*);
239 
240  // atom can be selected effitiantly if element selectors are initialised
241  inline const G4Element* SelectRandomAtom(const G4MaterialCutsCouple*,
242  const G4ParticleDefinition*,
243  G4double kineticEnergy,
244  G4double cutEnergy = 0.0,
245  G4double maxEnergy = DBL_MAX);
246 
247  // to select atom cross section per volume is recomputed for each element
248  const G4Element* SelectRandomAtom(const G4Material*,
249  const G4ParticleDefinition*,
250  G4double kineticEnergy,
251  G4double cutEnergy = 0.0,
252  G4double maxEnergy = DBL_MAX);
253 
254  //------------------------------------------------------------------------
255  // Get/Set methods
256  //------------------------------------------------------------------------
257 
259 
261 
263 
265 
267 
269 
270  inline G4double HighEnergyLimit() const;
271 
272  inline G4double LowEnergyLimit() const;
273 
274  inline G4double HighEnergyActivationLimit() const;
275 
276  inline G4double LowEnergyActivationLimit() const;
277 
278  inline G4double PolarAngleLimit() const;
279 
280  inline G4double SecondaryThreshold() const;
281 
282  inline G4bool LPMFlag() const;
283 
284  inline G4bool DeexcitationFlag() const;
285 
286  inline G4bool ForceBuildTableFlag() const;
287 
288  inline void SetHighEnergyLimit(G4double);
289 
290  inline void SetLowEnergyLimit(G4double);
291 
293 
295 
296  inline G4bool IsActive(G4double kinEnergy);
297 
298  inline void SetPolarAngleLimit(G4double);
299 
300  inline void SetSecondaryThreshold(G4double);
301 
302  inline void SetLPMFlag(G4bool val);
303 
304  inline void SetDeexcitationFlag(G4bool val);
305 
306  inline void ForceBuildTable(G4bool val);
307 
308  inline G4double MaxSecondaryKinEnergy(const G4DynamicParticle* dynParticle);
309 
310  inline const G4String& GetName() const;
311 
312  inline void SetCurrentCouple(const G4MaterialCutsCouple*);
313 
314  inline const G4Element* GetCurrentElement() const;
315 
316 protected:
317 
318  inline const G4MaterialCutsCouple* CurrentCouple() const;
319 
320  inline void SetCurrentElement(const G4Element*);
321 
322 private:
323 
324  // hide assignment operator
325  G4VEmModel & operator=(const G4VEmModel &right);
326  G4VEmModel(const G4VEmModel&);
327 
328  // ======== Parameters of the class fixed at construction =========
329 
330  G4VEmFluctuationModel* flucModel;
331  G4VEmAngularDistribution* anglModel;
332  const G4String name;
333 
334  // ======== Parameters of the class fixed at initialisation =======
335 
336  G4double lowLimit;
337  G4double highLimit;
338  G4double eMinActive;
339  G4double eMaxActive;
340  G4double polarAngleLimit;
341  G4double secondaryThreshold;
342  G4bool theLPMflag;
343  G4bool flagDeexcitation;
344  G4bool flagForceBuildTable;
345 
346  G4int nSelectors;
347  std::vector<G4EmElementSelector*> elmSelectors;
348 
349 protected:
350 
353  const std::vector<G4double>* theDensityFactor;
354  const std::vector<G4int>* theDensityIdx;
355 
356  // ======== Cashed values - may be state dependent ================
357 
358 private:
359 
360  const G4MaterialCutsCouple* fCurrentCouple;
361  const G4Element* fCurrentElement;
362 
363  G4int nsec;
364  std::vector<G4double> xsec;
365 
366 };
367 
368 // ======== Run time inline methods ================
369 
371 {
372  fCurrentCouple = p;
373 }
374 
375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
376 
378 {
379  return fCurrentCouple;
380 }
381 
382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
383 
385 {
386  fCurrentElement = elm;
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390 
392 {
393  return fCurrentElement;
394 }
395 
396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
397 
398 inline
400 {
401  return MaxSecondaryEnergy(dynPart->GetParticleDefinition(),
402  dynPart->GetKineticEnergy());
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
406 
408  const G4ParticleDefinition* p,
409  G4double kinEnergy,
410  G4double cutEnergy)
411 {
412  fCurrentCouple = c;
413  return ComputeDEDXPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy);
414 }
415 
416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
417 
419  const G4ParticleDefinition* p,
420  G4double kinEnergy,
421  G4double cutEnergy,
422  G4double maxEnergy)
423 {
424  fCurrentCouple = c;
425  return CrossSectionPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy,maxEnergy);
426 }
427 
428 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
429 
431  G4double ekin,
432  const G4Material* material,
433  G4double emin,
434  G4double emax)
435 {
436  G4double mfp = DBL_MAX;
437  G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
438  if (cross > DBL_MIN) { mfp = 1./cross; }
439  return mfp;
440 }
441 
442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
443 
445  const G4ParticleDefinition* part,
446  const G4Element* elm,
447  G4double kinEnergy,
448  G4double cutEnergy,
449  G4double maxEnergy)
450 {
451  fCurrentElement = elm;
452  return ComputeCrossSectionPerAtom(part,kinEnergy,elm->GetZ(),elm->GetN(),
453  cutEnergy,maxEnergy);
454 }
455 
456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
457 
458 inline const G4Element*
460  const G4ParticleDefinition* p,
461  G4double kinEnergy,
462  G4double cutEnergy,
463  G4double maxEnergy)
464 {
465  fCurrentCouple = couple;
466  if(nSelectors > 0) {
467  fCurrentElement =
468  elmSelectors[couple->GetIndex()]->SelectRandomAtom(kinEnergy);
469  } else {
470  fCurrentElement = SelectRandomAtom(couple->GetMaterial(),p,kinEnergy,
471  cutEnergy,maxEnergy);
472  }
473  return fCurrentElement;
474 }
475 
476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
477 
479 {
480  fCurrentElement = elm;
481  G4int N = G4int(elm->GetN() + 0.5);
482  G4int ni = elm->GetNumberOfIsotopes();
483  if(ni > 0) {
484  G4int idx = 0;
485  if(ni > 1) {
488  for(; idx<ni; ++idx) {
489  x -= ab[idx];
490  if (x <= 0.0) { break; }
491  }
492  if(idx >= ni) { idx = ni - 1; }
493  }
494  N = elm->GetIsotope(idx)->GetN();
495  }
496  return N;
497 }
498 
499 // ======== Get/Set inline methods used at initialisation ================
500 
502 {
503  return flucModel;
504 }
505 
506 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
507 
509 {
510  return anglModel;
511 }
512 
513 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
514 
516 {
517  anglModel = p;
518 }
519 
520 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
521 
523 {
524  return highLimit;
525 }
526 
527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
528 
530 {
531  return lowLimit;
532 }
533 
534 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
535 
537 {
538  return eMaxActive;
539 }
540 
541 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
542 
544 {
545  return eMinActive;
546 }
547 
548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
549 
551 {
552  return polarAngleLimit;
553 }
554 
555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
556 
558 {
559  return secondaryThreshold;
560 }
561 
562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
563 
565 {
566  return theLPMflag;
567 }
568 
569 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
570 
572 {
573  return flagDeexcitation;
574 }
575 
576 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
577 
579 {
580  return flagForceBuildTable;
581 }
582 
583 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
584 
586 {
587  highLimit = val;
588 }
589 
590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
591 
593 {
594  lowLimit = val;
595 }
596 
597 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
598 
600 {
601  eMaxActive = val;
602 }
603 
604 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
605 
607 {
608  eMinActive = val;
609 }
610 
611 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
612 
614 {
615  return (kinEnergy >= eMinActive && kinEnergy <= eMaxActive);
616 }
617 
618 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
619 
621 {
622  polarAngleLimit = val;
623 }
624 
625 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
626 
628 {
629  secondaryThreshold = val;
630 }
631 
632 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
633 
635 {
636  theLPMflag = val;
637 }
638 
639 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
640 
642 {
643  flagDeexcitation = val;
644 }
645 
646 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
647 
649 {
650  flagForceBuildTable = val;
651 }
652 
653 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
654 
655 inline const G4String& G4VEmModel::GetName() const
656 {
657  return name;
658 }
659 
661 {
662  return xSectionTable;
663 }
664 
665 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
666 
667 #endif
668