Geant4 11.1.1
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=nullptr, const G4String &processName="PenConversion")
 
virtual ~G4PenelopeGammaConversionModel ()
 
void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *) override
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
 
void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
 
void SetVerbosityLevel (G4int lev)
 
G4int GetVerbosityLevel ()
 
G4PenelopeGammaConversionModeloperator= (const G4PenelopeGammaConversionModel &right)=delete
 
 G4PenelopeGammaConversionModel (const G4PenelopeGammaConversionModel &)=delete
 
- 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 *, const G4double &length, G4double &eloss)
 
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 FillNumberOfSecondaries (G4int &numberOfTriplets, G4int &numberOfRecoil)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
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)
 
const G4ElementGetCurrentElement (const G4Material *mat=nullptr) const
 
G4int SelectRandomAtomNumber (const G4Material *) const
 
const G4IsotopeGetCurrentIsotope (const G4Element *elm=nullptr) const
 
G4int SelectIsotopeNumber (const G4Element *) const
 
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 *)
 
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 = nullptr
 
G4VParticleChangepParticleChange = nullptr
 
G4PhysicsTablexSectionTable = nullptr
 
const G4MaterialpBaseMaterial = nullptr
 
const std::vector< G4double > * theDensityFactor = nullptr
 
const std::vector< G4int > * theDensityIdx = nullptr
 
G4double inveplus
 
G4double pFactor = 1.0
 
size_t currentCoupleIndex = 0
 
size_t basedCoupleIndex = 0
 
G4bool lossFlucFlag = true
 

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 57 of file G4PenelopeGammaConversionModel.hh.

Constructor & Destructor Documentation

◆ G4PenelopeGammaConversionModel() [1/2]

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

Definition at line 85 of file G4PenelopeGammaConversionModel.cc.

87 :G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr),
88 fEffectiveCharge(nullptr),fMaterialInvScreeningRadius(nullptr),
89 fScreeningFunction(nullptr),fIsInitialised(false),fLocalTable(false)
90{
91 fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
92 fIntrinsicHighEnergyLimit = 100.0*GeV;
93 fSmallEnergy = 1.1*MeV;
94
95 if (part)
96 SetParticle(part);
97
98 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
99 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
100 //
101 fVerboseLevel= 0;
102 // Verbosity scale:
103 // 0 = nothing
104 // 1 = warning for energy non-conservation
105 // 2 = details of energy budget
106 // 3 = calculation of cross sections, file openings, sampling of atoms
107 // 4 = entering in methods
108}
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:746

◆ ~G4PenelopeGammaConversionModel()

G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel ( )
virtual

Definition at line 112 of file G4PenelopeGammaConversionModel.cc.

113{
114 //Delete shared tables, they exist only in the master model
115 if (IsMaster() || fLocalTable)
116 {
117 for(G4int i=0; i<=fMaxZ; ++i)
118 {
119 if(fLogAtomicCrossSection[i]) {
120 delete fLogAtomicCrossSection[i];
121 fLogAtomicCrossSection[i] = nullptr;
122 }
123 }
124 if (fEffectiveCharge)
125 delete fEffectiveCharge;
126 if (fMaterialInvScreeningRadius)
127 delete fMaterialInvScreeningRadius;
128 if (fScreeningFunction)
129 delete fScreeningFunction;
130 }
131}
int G4int
Definition: G4Types.hh:85
G4bool IsMaster() const
Definition: G4VEmModel.hh:725

◆ G4PenelopeGammaConversionModel() [2/2]

G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel ( const G4PenelopeGammaConversionModel )
delete

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

Reimplemented from G4VEmModel.

Definition at line 235 of file G4PenelopeGammaConversionModel.cc.

240{
241 //
242 // Penelope model v2008.
243 // Cross section (including triplet production) read from database and managed
244 // through the G4CrossSectionHandler utility. Cross section data are from
245 // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
246 //
247
248 if (energy < fIntrinsicLowEnergyLimit)
249 return 0;
250
251 G4int iZ = G4int(Z);
252
253 if (!fLogAtomicCrossSection[iZ])
254 {
255 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
256 //not filled up. This can happen in a UnitTest or via G4EmCalculator
257 if (fVerboseLevel > 0)
258 {
259 //Issue a G4Exception (warning) only in verbose mode
261 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
262 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
263 G4Exception("G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom()",
264 "em2018",JustWarning,ed);
265 }
266 //protect file reading via autolock
267 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
268 ReadDataFile(iZ);
269 lock.unlock();
270 fLocalTable = true;
271 }
272 G4double cs = 0;
273 G4double logene = G4Log(energy);
274 G4PhysicsFreeVector* theVec = fLogAtomicCrossSection[iZ];
275 G4double logXS = theVec->Value(logene);
276 cs = G4Exp(logXS);
277
278 if (fVerboseLevel > 2)
279 G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z <<
280 " = " << cs/barn << " barn" << G4endl;
281 return cs;
282}
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4double G4Log(G4double x)
Definition: G4Log.hh:227
double G4double
Definition: G4Types.hh:83
const G4int Z[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
G4double Value(const G4double energy, std::size_t &lastidx) const
G4double energy(const ThreeVector &p, const G4double m)

◆ GetVerbosityLevel()

G4int G4PenelopeGammaConversionModel::GetVerbosityLevel ( )
inline

Definition at line 82 of file G4PenelopeGammaConversionModel.hh.

82{return fVerboseLevel;};

◆ Initialise()

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

Implements G4VEmModel.

Definition at line 135 of file G4PenelopeGammaConversionModel.cc.

137{
138 if (fVerboseLevel > 3)
139 G4cout << "Calling G4PenelopeGammaConversionModel::Initialise()" << G4endl;
140
141 SetParticle(part);
142
143 //Only the master model creates/fills/destroys the tables
144 if (IsMaster() && part == fParticle)
145 {
146 //delete old material data...
147 if (fEffectiveCharge)
148 {
149 delete fEffectiveCharge;
150 fEffectiveCharge = nullptr;
151 }
152 if (fMaterialInvScreeningRadius)
153 {
154 delete fMaterialInvScreeningRadius;
155 fMaterialInvScreeningRadius = nullptr;
156 }
157 if (fScreeningFunction)
158 {
159 delete fScreeningFunction;
160 fScreeningFunction = nullptr;
161 }
162 //and create new ones
163 fEffectiveCharge = new std::map<const G4Material*,G4double>;
164 fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
165 fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
166
167 G4ProductionCutsTable* theCoupleTable =
169
170 for (G4int i=0;i<(G4int)theCoupleTable->GetTableSize();++i)
171 {
172 const G4Material* material =
173 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
174 const G4ElementVector* theElementVector = material->GetElementVector();
175
176 for (std::size_t j=0;j<material->GetNumberOfElements();++j)
177 {
178 G4int iZ = theElementVector->at(j)->GetZasInt();
179 //read data files only in the master
180 if (iZ <= fMaxZ && !fLogAtomicCrossSection[iZ])
181 ReadDataFile(iZ);
182 }
183
184 //check if material data are available
185 if (!fEffectiveCharge->count(material))
186 InitializeScreeningFunctions(material);
187 }
188 if (fVerboseLevel > 0) {
189 G4cout << "Penelope Gamma Conversion model v2008 is initialized " << G4endl
190 << "Energy range: "
191 << LowEnergyLimit() / MeV << " MeV - "
192 << HighEnergyLimit() / GeV << " GeV"
193 << G4endl;
194 }
195 }
196 if(fIsInitialised) return;
198 fIsInitialised = true;
199}
std::vector< const G4Element * > G4ElementVector
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:185
size_t GetNumberOfElements() const
Definition: G4Material.hh:181
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:124
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634

◆ InitialiseLocal()

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

Reimplemented from G4VEmModel.

Definition at line 203 of file G4PenelopeGammaConversionModel.cc.

205{
206 if (fVerboseLevel > 3)
207 G4cout << "Calling G4PenelopeGammaConversionModel::InitialiseLocal()" << G4endl;
208 //
209 //Check that particle matches: one might have multiple master models (e.g.
210 //for e+ and e-).
211 //
212 if (part == fParticle)
213 {
214 //Get the const table pointers from the master to the workers
215 const G4PenelopeGammaConversionModel* theModel =
216 static_cast<G4PenelopeGammaConversionModel*> (masterModel);
217
218 //Copy pointers to the data tables
219 fEffectiveCharge = theModel->fEffectiveCharge;
220 fMaterialInvScreeningRadius = theModel->fMaterialInvScreeningRadius;
221 fScreeningFunction = theModel->fScreeningFunction;
222 for(G4int i=0; i<=fMaxZ; ++i)
223 fLogAtomicCrossSection[i] = theModel->fLogAtomicCrossSection[i];
224
225 //Same verbosity for all workers, as the master
226 fVerboseLevel = theModel->fVerboseLevel;
227 }
228
229 return;
230}

◆ operator=()

G4PenelopeGammaConversionModel & G4PenelopeGammaConversionModel::operator= ( const G4PenelopeGammaConversionModel right)
delete

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 287 of file G4PenelopeGammaConversionModel.cc.

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

81{fVerboseLevel = lev;};

Member Data Documentation

◆ fParticle

const G4ParticleDefinition* G4PenelopeGammaConversionModel::fParticle
protected

Definition at line 90 of file G4PenelopeGammaConversionModel.hh.

Referenced by Initialise(), and InitialiseLocal().

◆ fParticleChange

G4ParticleChangeForGamma* G4PenelopeGammaConversionModel::fParticleChange
protected

Definition at line 89 of file G4PenelopeGammaConversionModel.hh.

Referenced by Initialise(), and SampleSecondaries().


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