Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4PenelopeGammaConversionModel Class Reference

#include <G4PenelopeGammaConversionModel.hh>

Inheritance diagram for G4PenelopeGammaConversionModel:
Collaboration diagram for G4PenelopeGammaConversionModel:

Public Member Functions

 G4PenelopeGammaConversionModel (const G4ParticleDefinition *p=0, const G4String &processName="PenConversion")
 
virtual ~G4PenelopeGammaConversionModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
void SetVerbosityLevel (G4int lev)
 
G4int GetVerbosityLevel ()
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector
< G4EmElementSelector * > * 
GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
virtual G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy)
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Attributes

G4ParticleChangeForGammafParticleChange
 
const G4ParticleDefinitionfParticle
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 57 of file G4PenelopeGammaConversionModel.hh.

Constructor & Destructor Documentation

G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel ( const G4ParticleDefinition p = 0,
const G4String processName = "PenConversion" 
)

Definition at line 56 of file G4PenelopeGammaConversionModel.cc.

59  logAtomicCrossSection(0),
60  fEffectiveCharge(0),fMaterialInvScreeningRadius(0),
61  fScreeningFunction(0),isInitialised(false),fLocalTable(false)
62 
63 {
64  fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
65  fIntrinsicHighEnergyLimit = 100.0*GeV;
66  fSmallEnergy = 1.1*MeV;
67 
68  InitializeScreeningRadii();
69 
70  if (part)
71  SetParticle(part);
72 
73  // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
74  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
75  //
76  verboseLevel= 0;
77  // Verbosity scale:
78  // 0 = nothing
79  // 1 = warning for energy non-conservation
80  // 2 = details of energy budget
81  // 3 = calculation of cross sections, file openings, sampling of atoms
82  // 4 = entering in methods
83 }
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
static constexpr double electron_mass_c2
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:724
static constexpr double GeV
Definition: G4SIunits.hh:217
static constexpr double MeV
Definition: G4SIunits.hh:214

Here is the call graph for this function:

G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel ( )
virtual

Definition at line 87 of file G4PenelopeGammaConversionModel.cc.

88 {
89  //Delete shared tables, they exist only in the master model
90  if (IsMaster() || fLocalTable)
91  {
92  if (logAtomicCrossSection)
93  {
94  for (auto& item : (*logAtomicCrossSection))
95  if (item.second) delete item.second;
96  delete logAtomicCrossSection;
97  }
98  if (fEffectiveCharge)
99  delete fEffectiveCharge;
100  if (fMaterialInvScreeningRadius)
101  delete fMaterialInvScreeningRadius;
102  if (fScreeningFunction)
103  delete fScreeningFunction;
104  }
105 }
G4bool IsMaster() const
Definition: G4VEmModel.hh:717

Here is the call graph for this function:

Member Function Documentation

G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition ,
G4double  kinEnergy,
G4double  Z,
G4double  A = 0,
G4double  cut = 0,
G4double  emax = DBL_MAX 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 219 of file G4PenelopeGammaConversionModel.cc.

224 {
225  //
226  // Penelope model v2008.
227  // Cross section (including triplet production) read from database and managed
228  // through the G4CrossSectionHandler utility. Cross section data are from
229  // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
230  //
231 
232  if (energy < fIntrinsicLowEnergyLimit)
233  return 0;
234 
235  G4int iZ = (G4int) Z;
236 
237  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
238  //not invoked
239  if (!logAtomicCrossSection)
240  {
241  //create a **thread-local** version of the table. Used only for G4EmCalculator and
242  //Unit Tests
243  fLocalTable = true;
244  logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
245  }
246  //now it should be ok
247  if (!logAtomicCrossSection->count(iZ))
248  {
249  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
250  //not filled up. This can happen in a UnitTest or via G4EmCalculator
251  if (verboseLevel > 0)
252  {
253  //Issue a G4Exception (warning) only in verbose mode
255  ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
256  ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
257  G4Exception("G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom()",
258  "em2018",JustWarning,ed);
259  }
260  //protect file reading via autolock
261  G4AutoLock lock(&PenelopeGammaConversionModelMutex);
262  ReadDataFile(iZ);
263  lock.unlock();
264  }
265 
266  G4double cs = 0;
267  G4double logene = std::log(energy);
268 
269  G4PhysicsFreeVector* theVec = logAtomicCrossSection->find(iZ)->second;
270 
271  G4double logXS = theVec->Value(logene);
272  cs = G4Exp(logXS);
273 
274  if (verboseLevel > 2)
275  G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z <<
276  " = " << cs/barn << " barn" << G4endl;
277  return cs;
278 }
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double energy(const ThreeVector &p, const G4double m)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
static constexpr double barn
Definition: G4SIunits.hh:105

Here is the call graph for this function:

G4int G4PenelopeGammaConversionModel::GetVerbosityLevel ( )
inline

Definition at line 84 of file G4PenelopeGammaConversionModel.hh.

84 {return verboseLevel;};
void G4PenelopeGammaConversionModel::Initialise ( const G4ParticleDefinition part,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 110 of file G4PenelopeGammaConversionModel.cc.

112 {
113  if (verboseLevel > 3)
114  G4cout << "Calling G4PenelopeGammaConversionModel::Initialise()" << G4endl;
115 
116  SetParticle(part);
117 
118  //Only the master model creates/fills/destroys the tables
119  if (IsMaster() && part == fParticle)
120  {
121  // logAtomicCrossSection is created only once, since it is never cleared
122  if (!logAtomicCrossSection)
123  logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
124 
125  //delete old material data...
126  if (fEffectiveCharge)
127  {
128  delete fEffectiveCharge;
129  fEffectiveCharge = nullptr;
130  }
131  if (fMaterialInvScreeningRadius)
132  {
133  delete fMaterialInvScreeningRadius;
134  fMaterialInvScreeningRadius = nullptr;
135  }
136  if (fScreeningFunction)
137  {
138  delete fScreeningFunction;
139  fScreeningFunction = nullptr;
140  }
141  //and create new ones
142  fEffectiveCharge = new std::map<const G4Material*,G4double>;
143  fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
144  fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
145 
146  G4ProductionCutsTable* theCoupleTable =
148 
149  for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
150  {
151  const G4Material* material =
152  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
153  const G4ElementVector* theElementVector = material->GetElementVector();
154 
155  for (size_t j=0;j<material->GetNumberOfElements();j++)
156  {
157  G4int iZ = (G4int) theElementVector->at(j)->GetZ();
158  //read data files only in the master
159  if (!logAtomicCrossSection->count(iZ))
160  ReadDataFile(iZ);
161  }
162 
163  //check if material data are available
164  if (!fEffectiveCharge->count(material))
165  InitializeScreeningFunctions(material);
166  }
167 
168 
169  if (verboseLevel > 0) {
170  G4cout << "Penelope Gamma Conversion model v2008 is initialized " << G4endl
171  << "Energy range: "
172  << LowEnergyLimit() / MeV << " MeV - "
173  << HighEnergyLimit() / GeV << " GeV"
174  << G4endl;
175  }
176 
177  }
178 
179 
180  if(isInitialised) return;
182  isInitialised = true;
183 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
std::vector< G4Element * > G4ElementVector
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
G4bool IsMaster() const
Definition: G4VEmModel.hh:717
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
const G4Material * GetMaterial() const

Here is the call graph for this function:

void G4PenelopeGammaConversionModel::InitialiseLocal ( const G4ParticleDefinition part,
G4VEmModel masterModel 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 187 of file G4PenelopeGammaConversionModel.cc.

189 {
190  if (verboseLevel > 3)
191  G4cout << "Calling G4PenelopeGammaConversionModel::InitialiseLocal()" << G4endl;
192 
193  //
194  //Check that particle matches: one might have multiple master models (e.g.
195  //for e+ and e-).
196  //
197  if (part == fParticle)
198  {
199  //Get the const table pointers from the master to the workers
200  const G4PenelopeGammaConversionModel* theModel =
201  static_cast<G4PenelopeGammaConversionModel*> (masterModel);
202 
203  //Copy pointers to the data tables
204  fEffectiveCharge = theModel->fEffectiveCharge;
205  fMaterialInvScreeningRadius = theModel->fMaterialInvScreeningRadius;
206  fScreeningFunction = theModel->fScreeningFunction;
207  logAtomicCrossSection = theModel->logAtomicCrossSection;
208 
209  //Same verbosity for all workers, as the master
210  verboseLevel = theModel->verboseLevel;
211  }
212 
213  return;
214 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
void G4PenelopeGammaConversionModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  fvect,
const G4MaterialCutsCouple couple,
const G4DynamicParticle aDynamicGamma,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 283 of file G4PenelopeGammaConversionModel.cc.

288 {
289  //
290  // Penelope model v2008.
291  // Final state is sampled according to the Bethe-Heitler model with Coulomb
292  // corrections, according to the semi-empirical model of
293  // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
294  //
295  // The model uses the high energy Coulomb correction from
296  // H. Davies et al., Phys. Rev. 93 (1954) 788
297  // and atomic screening radii tabulated from
298  // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
299  // for Z= 1 to 92.
300  //
301  if (verboseLevel > 3)
302  G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl;
303 
304  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
305 
306  // Always kill primary
309 
310  if (photonEnergy <= fIntrinsicLowEnergyLimit)
311  {
313  return ;
314  }
315 
316  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
317  const G4Material* mat = couple->GetMaterial();
318 
319  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
320  //not invoked
321  if (!fEffectiveCharge)
322  {
323  //create a **thread-local** version of the table. Used only for G4EmCalculator and
324  //Unit Tests
325  fLocalTable = true;
326  fEffectiveCharge = new std::map<const G4Material*,G4double>;
327  fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
328  fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
329  }
330 
331  if (!fEffectiveCharge->count(mat))
332  {
333  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
334  //not filled up. This can happen in a UnitTest or via G4EmCalculator
335  if (verboseLevel > 0)
336  {
337  //Issue a G4Exception (warning) only in verbose mode
339  ed << "Unable to allocate the EffectiveCharge data for " <<
340  mat->GetName() << G4endl;
341  ed << "This can happen only in Unit Tests" << G4endl;
342  G4Exception("G4PenelopeGammaConversionModel::SampleSecondaries()",
343  "em2019",JustWarning,ed);
344  }
345  //protect file reading via autolock
346  G4AutoLock lock(&PenelopeGammaConversionModelMutex);
347  InitializeScreeningFunctions(mat);
348  lock.unlock();
349  }
350 
351  // eps is the fraction of the photon energy assigned to e- (including rest mass)
352  G4double eps = 0;
353  G4double eki = electron_mass_c2/photonEnergy;
354 
355  //Do it fast for photon energy < 1.1 MeV (close to threshold)
356  if (photonEnergy < fSmallEnergy)
357  eps = eki + (1.0-2.0*eki)*G4UniformRand();
358  else
359  {
360  //Complete calculation
361  G4double effC = fEffectiveCharge->find(mat)->second;
362  G4double alz = effC*fine_structure_const;
363  G4double T = std::sqrt(2.0*eki);
364  G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
365  +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
366  -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
367  +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
368 
369  G4double F0b = fScreeningFunction->find(mat)->second.second;
370  G4double g0 = F0b + F00;
371  G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
372  G4double bmin = 4.0*eki/invRad;
373  std::pair<G4double,G4double> scree = GetScreeningFunctions(bmin);
374  G4double g1 = scree.first;
375  G4double g2 = scree.second;
376  G4double g1min = g1+g0;
377  G4double g2min = g2+g0;
378  G4double xr = 0.5-eki;
379  G4double a1 = 2.*g1min*xr*xr/3.;
380  G4double p1 = a1/(a1+g2min);
381 
382  G4bool loopAgain = false;
383  //Random sampling of eps
384  do{
385  loopAgain = false;
386  if (G4UniformRand() <= p1)
387  {
388  G4double ru2m1 = 2.0*G4UniformRand()-1.0;
389  if (ru2m1 < 0)
390  eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
391  else
392  eps = 0.5+xr*std::pow(ru2m1,1./3.);
393  G4double B = eki/(invRad*eps*(1.0-eps));
394  scree = GetScreeningFunctions(B);
395  g1 = scree.first;
396  g1 = std::max(g1+g0,0.);
397  if (G4UniformRand()*g1min > g1)
398  loopAgain = true;
399  }
400  else
401  {
402  eps = eki+2.0*xr*G4UniformRand();
403  G4double B = eki/(invRad*eps*(1.0-eps));
404  scree = GetScreeningFunctions(B);
405  g2 = scree.second;
406  g2 = std::max(g2+g0,0.);
407  if (G4UniformRand()*g2min > g2)
408  loopAgain = true;
409  }
410  }while(loopAgain);
411 
412  }
413  if (verboseLevel > 4)
414  G4cout << "Sampled eps = " << eps << G4endl;
415 
416  G4double electronTotEnergy = eps*photonEnergy;
417  G4double positronTotEnergy = (1.0-eps)*photonEnergy;
418 
419  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
420 
421  //electron kinematics
422  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
423  G4double costheta_el = G4UniformRand()*2.0-1.0;
424  G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
425  costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
426  G4double phi_el = twopi * G4UniformRand() ;
427  G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
428  G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
429  G4double dirZ_el = costheta_el;
430 
431  //positron kinematics
432  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
433  G4double costheta_po = G4UniformRand()*2.0-1.0;
434  kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
435  costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
436  G4double phi_po = twopi * G4UniformRand() ;
437  G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
438  G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
439  G4double dirZ_po = costheta_po;
440 
441  // Kinematics of the created pair:
442  // the electron and positron are assumed to have a symetric angular
443  // distribution with respect to the Z axis along the parent photon
444  G4double localEnergyDeposit = 0. ;
445 
446  if (electronKineEnergy > 0.0)
447  {
448  G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
449  electronDirection.rotateUz(photonDirection);
451  electronDirection,
452  electronKineEnergy);
453  fvect->push_back(electron);
454  }
455  else
456  {
457  localEnergyDeposit += electronKineEnergy;
458  electronKineEnergy = 0;
459  }
460 
461  //Generate the positron. Real particle in any case, because it will annihilate. If below
462  //threshold, produce it at rest
463  if (positronKineEnergy < 0.0)
464  {
465  localEnergyDeposit += positronKineEnergy;
466  positronKineEnergy = 0; //produce it at rest
467  }
468  G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
469  positronDirection.rotateUz(photonDirection);
471  positronDirection, positronKineEnergy);
472  fvect->push_back(positron);
473 
474  //Add rest of energy to the local energy deposit
475  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
476 
477  if (verboseLevel > 1)
478  {
479  G4cout << "-----------------------------------------------------------" << G4endl;
480  G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl;
481  G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
482  G4cout << "-----------------------------------------------------------" << G4endl;
483  if (electronKineEnergy)
484  G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV"
485  << G4endl;
486  if (positronKineEnergy)
487  G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
488  G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
489  if (localEnergyDeposit)
490  G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
491  G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
492  localEnergyDeposit+2.0*electron_mass_c2)/keV <<
493  " keV" << G4endl;
494  G4cout << "-----------------------------------------------------------" << G4endl;
495  }
496  if (verboseLevel > 0)
497  {
498  G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
499  localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
500  if (energyDiff > 0.05*keV)
501  G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: "
502  << (electronKineEnergy+positronKineEnergy+
503  localEnergyDeposit+2.0*electron_mass_c2)/keV
504  << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
505  }
506 }
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
#define F00
const G4String & GetName() const
Definition: G4Material.hh:178
static const G4double eps
double B(double temperature)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
static constexpr double electron_mass_c2
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4Positron * Positron()
Definition: G4Positron.cc:94
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
static constexpr double fine_structure_const
static constexpr double keV
Definition: G4SIunits.hh:216
const G4Material * GetMaterial() const

Here is the call graph for this function:

void G4PenelopeGammaConversionModel::SetVerbosityLevel ( G4int  lev)
inline

Definition at line 83 of file G4PenelopeGammaConversionModel.hh.

83 {verboseLevel = lev;};

Member Data Documentation

const G4ParticleDefinition* G4PenelopeGammaConversionModel::fParticle
protected

Definition at line 88 of file G4PenelopeGammaConversionModel.hh.

G4ParticleChangeForGamma* G4PenelopeGammaConversionModel::fParticleChange
protected

Definition at line 84 of file G4PenelopeGammaConversionModel.hh.


The documentation for this class was generated from the following files: