Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4LivermoreGammaConversionModelRC Class Reference

#include <G4LivermoreGammaConversionModelRC.hh>

+ Inheritance diagram for G4LivermoreGammaConversionModelRC:

Public Member Functions

 G4LivermoreGammaConversionModelRC (const G4ParticleDefinition *p=0, const G4String &nam="LivermoreGammaConversionRC_1")
 
virtual ~G4LivermoreGammaConversionModelRC ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
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)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double)
 
- 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 G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (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)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
G4VEmModel * GetTripletModel ()
 
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 G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4Element * GetCurrentElement () const
 
const G4Isotope * GetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 
G4VEmModel & operator= (const G4VEmModel &right)=delete
 
 G4VEmModel (const G4VEmModel &)=delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const G4Material * pBaseMaterial
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 
G4double inveplus
 
G4double pFactor
 

Detailed Description

Definition at line 39 of file G4LivermoreGammaConversionModelRC.hh.

Constructor & Destructor Documentation

◆ G4LivermoreGammaConversionModelRC()

G4LivermoreGammaConversionModelRC::G4LivermoreGammaConversionModelRC ( const G4ParticleDefinition *  p = 0,
const G4String &  nam = "LivermoreGammaConversionRC_1" 
)

Definition at line 50 of file G4LivermoreGammaConversionModelRC.cc.

52:G4VEmModel(nam),isInitialised(false),smallEnergy(2.*MeV)
53{
54 fParticleChange = nullptr;
55
56 lowEnergyLimit = 2.0*electron_mass_c2;
57
58 verboseLevel= 0;
59 // Verbosity scale for debugging purposes:
60 // 0 = nothing
61 // 1 = calculation of cross sections, file openings...
62 // 2 = entering in methods
63
64 if(verboseLevel > 0)
65 {
66 G4cout << "G4LivermoreGammaConversionModelRC is constructed " << G4endl;
67 }
68}
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout

◆ ~G4LivermoreGammaConversionModelRC()

G4LivermoreGammaConversionModelRC::~G4LivermoreGammaConversionModelRC ( )
virtual

Definition at line 72 of file G4LivermoreGammaConversionModelRC.cc.

73{
74 if(IsMaster()) {
75 for(G4int i=0; i<maxZ; ++i) {
76 if(data[i]) {
77 delete data[i];
78 data[i] = 0;
79 }
80 }
81 }
82}
G4int
int G4int
Definition: G4Types.hh:85
G4VEmModel::IsMaster
G4bool IsMaster() const
Definition: G4VEmModel.hh:736

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

G4double G4LivermoreGammaConversionModelRC::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 G4LivermoreGammaConversionModelRC.cc.

223{
224 if (verboseLevel > 1)
225 {
226 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModelRC"
227 << G4endl;
228 }
229
230 if (GammaEnergy < lowEnergyLimit) { return 0.0; }
231
232 G4double xs = 0.0;
233
234 G4int intZ=G4int(Z);
235
236 if(intZ < 1 || intZ > maxZ) { return xs; }
237
238 G4LPhysicsFreeVector* pv = data[intZ];
239
240 // if element was not initialised
241 // do initialisation safely for MT mode
242 if(!pv)
243 {
244 InitialiseForElement(0, intZ);
245 pv = data[intZ];
246 if(!pv) { return xs; }
247 }
248 // x-section is taken from the table
249 xs = pv->Value(GammaEnergy);
250
251 if(verboseLevel > 0)
252 {
253 G4int n = pv->GetVectorLength() - 1;
254 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
255 << GammaEnergy/MeV << G4endl;
256 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
257 G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
258 G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
259 G4cout << "*********************************************************" << G4endl;
260 }
261
262 return xs;
263
264}
G4double
double G4double
Definition: G4Types.hh:83
G4LivermoreGammaConversionModelRC::InitialiseForElement
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
Definition: G4LivermoreGammaConversionModelRC.cc:604
G4PhysicsVector::Value
G4double Value(G4double theEnergy, std::size_t &lastidx) const
Definition: G4PhysicsVector.cc:519
G4PhysicsVector::GetVectorLength
std::size_t GetVectorLength() const

◆ Initialise()

void G4LivermoreGammaConversionModelRC::Initialise ( const G4ParticleDefinition *  particle,
const G4DataVector &  cuts 
)
virtual

Implements G4VEmModel.

Definition at line 86 of file G4LivermoreGammaConversionModelRC.cc.

89{
90 if (verboseLevel > 1)
91 {
92 G4cout << "Calling Initialise() of G4LivermoreGammaConversionModelRC."
93 << G4endl
94 << "Energy range: "
95 << LowEnergyLimit() / MeV << " MeV - "
96 << HighEnergyLimit() / GeV << " GeV"
97 << G4endl;
98 }
99
100 if(IsMaster())
101 {
102
103 // Initialise element selector
104
105 InitialiseElementSelectors(particle, cuts);
106
107 // Access to elements
108
109 char* path = std::getenv("G4LEDATA");
110
111 G4ProductionCutsTable* theCoupleTable =
112 G4ProductionCutsTable::GetProductionCutsTable();
113
114 G4int numOfCouples = theCoupleTable->GetTableSize();
115
116 for(G4int i=0; i<numOfCouples; ++i)
117 {
118 const G4Material* material =
119 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
120 const G4ElementVector* theElementVector = material->GetElementVector();
121 G4int nelm = material->GetNumberOfElements();
122
123 for (G4int j=0; j<nelm; ++j)
124 {
125 G4int Z = (G4int)(*theElementVector)[j]->GetZ();
126 if(Z < 1) { Z = 1; }
127 else if(Z > maxZ) { Z = maxZ; }
128 if(!data[Z]) { ReadData(Z, path); }
129 }
130 }
131 }
132 if(isInitialised) { return; }
133 fParticleChange = GetParticleChangeForGamma();
134 isInitialised = true;
135}
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:43
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:246
G4ProductionCutsTable::GetTableSize
std::size_t GetTableSize() const
Definition: G4ProductionCutsTable.hh:239
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:58
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:133
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:652
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:645
G4VEmModel::InitialiseElementSelectors
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:148

◆ InitialiseForElement()

void G4LivermoreGammaConversionModelRC::InitialiseForElement ( const G4ParticleDefinition *  ,
G4int  Z 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 604 of file G4LivermoreGammaConversionModelRC.cc.

607{
608 G4AutoLock l(&LivermoreGammaConversionModelRCMutex);
609 // G4cout << "G4LivermoreGammaConversionModelRC::InitialiseForElement Z= "
610 // << Z << G4endl;
611 if(!data[Z]) { ReadData(Z); }
612 l.unlock();
613}

Referenced by ComputeCrossSectionPerAtom().

◆ InitialiseLocal()

void G4LivermoreGammaConversionModelRC::InitialiseLocal ( const G4ParticleDefinition *  ,
G4VEmModel *  masterModel 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 139 of file G4LivermoreGammaConversionModelRC.cc.

141{
142 SetElementSelectors(masterModel->GetElementSelectors());
143}
G4VEmModel::SetElementSelectors
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:842
G4VEmModel::GetElementSelectors
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:834

◆ MinPrimaryEnergy()

G4double G4LivermoreGammaConversionModelRC::MinPrimaryEnergy ( const G4Material *  ,
const G4ParticleDefinition *  ,
G4double   
)
virtual

Reimplemented from G4VEmModel.

Definition at line 148 of file G4LivermoreGammaConversionModelRC.cc.

151{
152 return lowEnergyLimit;
153}

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 268 of file G4LivermoreGammaConversionModelRC.cc.

273{
274
275// The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
276// cross sections with Coulomb correction. A modified version of the random
277// number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
278
279// Note 1 : Effects due to the breakdown of the Born approximation at low
280// energy are ignored.
281// Note 2 : The differential cross section implicitly takes account of
282// pair creation in both nuclear and atomic electron fields. However triplet
283// prodution is not generated.
284
285 if (verboseLevel > 1) {
286 G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModelRC"
287 << G4endl;
288 }
289
290 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
291 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
292
293 G4double epsilon ;
294 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
295
296 G4double electronTotEnergy = 0.0;
297 G4double positronTotEnergy = 0.0;
298 G4double HardPhotonEnergy = 0.0;
299
300
301 // Do it fast if photon energy < 2. MeV
302 if (photonEnergy < smallEnergy )
303 {
304 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
305 if (G4UniformRand() > 0.5)
306 {
307 electronTotEnergy = (1. - epsilon) * photonEnergy;
308 positronTotEnergy = epsilon * photonEnergy;
309 }
310 else
311 {
312 positronTotEnergy = (1. - epsilon) * photonEnergy;
313 electronTotEnergy = epsilon * photonEnergy;
314 }
315 }
316 else
317 {
318 // Select randomly one element in the current material
319
320 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
321 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
322
323 if (element == 0)
324 {
325 G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries - element = 0"
326 << G4endl;
327 return;
328 }
329 G4IonisParamElm* ionisation = element->GetIonisation();
330 if (ionisation == 0)
331 {
332 G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries - ionisation = 0"
333 << G4endl;
334 return;
335 }
336
337 // Extract Coulomb factor for this Element
338 G4double fZ = 8. * (ionisation->GetlogZ3());
339 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
340
341 // Limits of the screening variable
342 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
343 G4double screenMax = G4Exp ((42.24 - fZ)/8.368) - 0.952 ;
344 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
345
346 // Limits of the energy sampling
347 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
348 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
349 G4double epsilonRange = 0.5 - epsilonMin ;
350
351 // Sample the energy rate of the created electron (or positron)
352 G4double screen;
353 G4double gReject ;
354
355 G4double f10 = ScreenFunction1(screenMin) - fZ;
356 G4double f20 = ScreenFunction2(screenMin) - fZ;
357 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
358 G4double normF2 = std::max(1.5 * f20,0.);
359
360 // Method for Radiative corrections
361
362 G4double a=393.3750918, b=115.3070201, c=810.6428451, d=19.96497475, e=1016.874592, f=1.936685510,
363 gLocal=751.2140962, h=0.099751048, i=299.9466339, j=0.002057250, k=49.81034926;
364 G4double aa=-18.6371131, bb=-1729.95248, cc=9450.971186, dd=106336.0145, ee=55143.09287, ff=-117602.840,
365 gg=-721455.467, hh=693957.8635, ii=156266.1085, jj=533209.9347;
366 G4double Rechazo = 0.;
367 G4double logepsMin = log(epsilonMin);
368 G4double NormaRC = a + b*logepsMin + c/logepsMin + d*pow(logepsMin,2.) + e/pow(logepsMin,2.) + f*pow(logepsMin,3.) +
369 gLocal/pow(logepsMin,3.) + h*pow(logepsMin,4.) + i/pow(logepsMin,4.) + j*pow(logepsMin,5.) +
370 k/pow(logepsMin,5.);
371
372 G4double HardPhotonThreshold = 0.08;
373 G4double r1, r2, r3, beta=0, gbeta, sigt = 582.068, sigh, rejet;
374 // , Pi = 2.*acos(0.);
375 G4double cg = (11./2.)/(G4Exp(-11.*HardPhotonThreshold/2.)-G4Exp(-11./2.));
376
377 r1 = G4UniformRand();
378 sigh = 1028.58*G4Exp(-HardPhotonThreshold/0.09033) + 136.63; // sigma hard
379
380
381 if (r1 > 1.- sigh/sigt) {
382 r2 = G4UniformRand();
383 rejet = 0.;
384 while (r2 > rejet) {
385 r3 = G4UniformRand();
386 beta = (-2./11.)*log(G4Exp(-0.08*11./2.)-r3*11./(2.*cg));
387 gbeta = G4Exp(-11.*beta/2.);
388 rejet = fbeta(beta)/(8000.*gbeta);
389 }
390 HardPhotonEnergy = beta * photonEnergy;
391 }
392 else{
393 HardPhotonEnergy = 0.;
394 }
395
396 photonEnergy -= HardPhotonEnergy;
397
398 do
399 {
400 do
401 {
402 if (normF1 / (normF1 + normF2) > G4UniformRand() )
403 {
404 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.333333) ;
405 screen = screenFactor / (epsilon * (1. - epsilon));
406 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
407 }
408 else
409 {
410 epsilon = epsilonMin + epsilonRange * G4UniformRand();
411 screen = screenFactor / (epsilon * (1 - epsilon));
412 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
413 }
414 } while ( gReject < G4UniformRand() );
415
416
417
418 if (G4UniformRand()>0.5) epsilon = (1. - epsilon); // Extención de Epsilon hasta 1.
419
420 G4double logepsilon = log(epsilon);
421 G4double deltaP_R1 = 1. + (a + b*logepsilon + c/logepsilon + d*pow(logepsilon,2.) + e/pow(logepsilon,2.) +
422 f*pow(logepsilon,3.) + gLocal/pow(logepsilon,3.) + h*pow(logepsilon,4.) + i/pow(logepsilon,4.) +
423 j*pow(logepsilon,5.) + k/pow(logepsilon,5.))/100.;
424 G4double deltaP_R2 = 1.+((aa + cc*logepsilon + ee*pow(logepsilon,2.) + gg*pow(logepsilon,3.) + ii*pow(logepsilon,4.))
425 / (1. + bb*logepsilon + dd*pow(logepsilon,2.) + ff*pow(logepsilon,3.) + hh*pow(logepsilon,4.)
426 + jj*pow(logepsilon,5.) ))/100.;
427
428 if (epsilon <= 0.5)
429 {
430 Rechazo = deltaP_R1/NormaRC;
431 }
432 else
433 {
434 Rechazo = deltaP_R2/NormaRC;
435 }
436 //G4cout << Rechazo << " " << NormaRC << " " << epsilon << G4endl;
437
438 } while (Rechazo < G4UniformRand() );
439
440 electronTotEnergy = (1. - epsilon) * photonEnergy;
441 positronTotEnergy = epsilon * photonEnergy;
442
443 } // End of epsilon sampling
444
445
446 // Fix charges randomly
447
448 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
449 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
450 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
451
452 G4double u;
453 const G4double a1 = 0.625;
454 G4double a2 = 3. * a1;
455 // G4double d = 27. ;
456
457 // if (9. / (9. + d) > G4UniformRand())
458 if (0.25 > G4UniformRand())
459 {
460 u = - G4Log(G4UniformRand() * G4UniformRand()) / a1 ;
461 }
462 else
463 {
464 u = - G4Log(G4UniformRand() * G4UniformRand()) / a2 ;
465 }
466
467 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
468 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
469 G4double phi = twopi * G4UniformRand();
470
471 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
472 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
473
474
475 // Kinematics of the created pair:
476 // the electron and positron are assumed to have a symetric angular
477 // distribution with respect to the Z axis along the parent photon
478
479 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
480
481 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
482 electronDirection.rotateUz(photonDirection);
483
484 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
485 electronDirection,
486 electronKineEnergy);
487
488 // The e+ is always created
489 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
490
491 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
492 positronDirection.rotateUz(photonDirection);
493
494 // Create G4DynamicParticle object for the particle2
495 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
496 positronDirection,
497 positronKineEnergy);
498 // Fill output vector
499 fvect->push_back(particle1);
500 fvect->push_back(particle2);
501
502 if (HardPhotonEnergy > 0.)
503 {
504 G4double thetaHardPhoton = u*electron_mass_c2/HardPhotonEnergy;
505 phi = twopi * G4UniformRand();
506 G4double dxHardP= std::sin(thetaHardPhoton)*std::cos(phi);
507 G4double dyHardP= std::sin(thetaHardPhoton)*std::sin(phi);
508 G4double dzHardP =std::cos(thetaHardPhoton);
509
510 G4ThreeVector hardPhotonDirection (dxHardP, dyHardP, dzHardP);
511 hardPhotonDirection.rotateUz(photonDirection);
512 G4DynamicParticle* particle3 = new G4DynamicParticle (G4Gamma::Gamma(),
513 hardPhotonDirection,
514 HardPhotonEnergy);
515 fvect->push_back(particle3);
516 }
517
518 // kill incident photon
519 fParticleChange->SetProposedKineticEnergy(0.);
520 fParticleChange->ProposeTrackStatus(fStopAndKill);
521
522}
epsilon
double epsilon(double density, double temperature)
Definition: G4DNAElectronHoleRecombination.cc:96
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:226
fStopAndKill
@ fStopAndKill
Definition: G4TrackStatus.hh:46
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4Element::GetfCoulomb
G4double GetfCoulomb() const
Definition: G4Element.hh:190
G4Element::GetIonisation
G4IonisParamElm * GetIonisation() const
Definition: G4Element.hh:198
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:85
G4IonisParamElm::GetlogZ3
G4double GetlogZ3() const
Definition: G4IonisParamElm.hh:62
G4IonisParamElm::GetZ3
G4double GetZ3() const
Definition: G4IonisParamElm.hh:60
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:119
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:93
G4VEmModel::SelectRandomAtom
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:570
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
The documentation for this class was generated from the following files: