Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PenelopeGammaConversionModel Class Reference

#include <G4PenelopeGammaConversionModel.hh>

+ Inheritance 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 Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
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 ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
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)
 
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, 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 *)
 
G4int SelectIsotopeNumber (const G4Element *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
G4VEmModelGetTripletModel ()
 
void SetTripletModel (G4VEmModel *)
 
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) const
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
void SetUseBaseMaterials (G4bool val)
 
G4bool UseBaseMaterials () 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)
 
G4VEmModeloperator= (const G4VEmModel &right)=delete
 
 G4VEmModel (const G4VEmModel &)=delete
 

Protected Attributes

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

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 *)
 

Detailed Description

Definition at line 56 of file G4PenelopeGammaConversionModel.hh.

Constructor & Destructor Documentation

◆ G4PenelopeGammaConversionModel()

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

Definition at line 55 of file G4PenelopeGammaConversionModel.cc.

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

◆ ~G4PenelopeGammaConversionModel()

G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel ( )
virtual

Definition at line 86 of file G4PenelopeGammaConversionModel.cc.

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

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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 218 of file G4PenelopeGammaConversionModel.cc.

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

◆ GetVerbosityLevel()

G4int G4PenelopeGammaConversionModel::GetVerbosityLevel ( )
inline

Definition at line 83 of file G4PenelopeGammaConversionModel.hh.

83{return verboseLevel;};

◆ Initialise()

void G4PenelopeGammaConversionModel::Initialise ( const G4ParticleDefinition part,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 109 of file G4PenelopeGammaConversionModel.cc.

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

◆ InitialiseLocal()

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

Reimplemented from G4VEmModel.

Definition at line 186 of file G4PenelopeGammaConversionModel.cc.

188{
189 if (verboseLevel > 3)
190 G4cout << "Calling G4PenelopeGammaConversionModel::InitialiseLocal()" << G4endl;
191
192 //
193 //Check that particle matches: one might have multiple master models (e.g.
194 //for e+ and e-).
195 //
196 if (part == fParticle)
197 {
198 //Get the const table pointers from the master to the workers
199 const G4PenelopeGammaConversionModel* theModel =
200 static_cast<G4PenelopeGammaConversionModel*> (masterModel);
201
202 //Copy pointers to the data tables
203 fEffectiveCharge = theModel->fEffectiveCharge;
204 fMaterialInvScreeningRadius = theModel->fMaterialInvScreeningRadius;
205 fScreeningFunction = theModel->fScreeningFunction;
206 logAtomicCrossSection = theModel->logAtomicCrossSection;
207
208 //Same verbosity for all workers, as the master
209 verboseLevel = theModel->verboseLevel;
210 }
211
212 return;
213}

◆ SampleSecondaries()

void G4PenelopeGammaConversionModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  fvect,
const G4MaterialCutsCouple couple,
const G4DynamicParticle aDynamicGamma,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 282 of file G4PenelopeGammaConversionModel.cc.

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

◆ SetVerbosityLevel()

void G4PenelopeGammaConversionModel::SetVerbosityLevel ( G4int  lev)
inline

Definition at line 82 of file G4PenelopeGammaConversionModel.hh.

82{verboseLevel = lev;};

Member Data Documentation

◆ fParticle

const G4ParticleDefinition* G4PenelopeGammaConversionModel::fParticle
protected

Definition at line 87 of file G4PenelopeGammaConversionModel.hh.

Referenced by Initialise(), and InitialiseLocal().

◆ fParticleChange

G4ParticleChangeForGamma* G4PenelopeGammaConversionModel::fParticleChange
protected

Definition at line 86 of file G4PenelopeGammaConversionModel.hh.

Referenced by Initialise(), and SampleSecondaries().


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