57 1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,
58 6.4696e+01,6.1228e+01,5.7524e+01,5.4033e+01,
59 5.0787e+01,4.7851e+01,4.6373e+01,4.5401e+01,
60 4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,
61 4.1586e+01,4.0953e+01,4.0524e+01,4.0256e+01,
62 3.9756e+01,3.9144e+01,3.8462e+01,3.7778e+01,
63 3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,
64 3.4688e+01,3.4197e+01,3.3786e+01,3.3422e+01,
65 3.3068e+01,3.2740e+01,3.2438e+01,3.2143e+01,
66 3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,
67 3.0950e+01,3.0758e+01,3.0561e+01,3.0285e+01,
68 3.0097e+01,2.9832e+01,2.9581e+01,2.9411e+01,
69 2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,
70 2.8580e+01,2.8442e+01,2.8312e+01,2.8139e+01,
71 2.7973e+01,2.7819e+01,2.7675e+01,2.7496e+01,
72 2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,
73 2.6516e+01,2.6304e+01,2.6108e+01,2.5929e+01,
74 2.5730e+01,2.5577e+01,2.5403e+01,2.5245e+01,
75 2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,
76 2.4506e+01,2.4391e+01,2.4262e+01,2.4145e+01,
77 2.4039e+01,2.3922e+01,2.3813e+01,2.3712e+01,
78 2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,
79 2.3238e+01,2.3139e+01,2.3048e+01,2.2967e+01,
80 2.2833e+01,2.2694e+01,2.2624e+01,2.2545e+01,
81 2.2446e+01,2.2358e+01,2.2264e+01};
87 :
G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr),
88 fEffectiveCharge(nullptr),fMaterialInvScreeningRadius(nullptr),
89 fScreeningFunction(nullptr),fIsInitialised(false),fLocalTable(false)
139 G4cout <<
"Calling G4PenelopeGammaConversionModel::Initialise()" <<
G4endl;
176 for (
size_t j=0;j<
material->GetNumberOfElements();j++)
178 G4int iZ = theElementVector->at(j)->GetZasInt();
189 G4cout <<
"Penelope Gamma Conversion model v2008 is initialized " <<
G4endl
207 G4cout <<
"Calling G4PenelopeGammaConversionModel::InitialiseLocal()" <<
G4endl;
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()",
279 G4cout <<
"Gamma conversion cross section at " <<
energy/
MeV <<
" MeV for Z=" <<
Z <<
306 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" <<
G4endl;
343 ed <<
"Unable to allocate the EffectiveCharge data for " <<
345 ed <<
"This can happen only in Unit Tests" <<
G4endl;
346 G4Exception(
"G4PenelopeGammaConversionModel::SampleSecondaries()",
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;
394 eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
396 eps = 0.5+xr*std::pow(ru2m1,1./3.);
420 G4double positronTotEnergy = (1.0-
eps)*photonEnergy;
428 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
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);
437 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*
electron_mass_c2));
438 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
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);
449 if (electronKineEnergy > 0.0)
451 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
452 electronDirection.
rotateUz(photonDirection);
460 localEnergyDeposit += electronKineEnergy;
461 electronKineEnergy = 0;
466 if (positronKineEnergy < 0.0)
468 localEnergyDeposit += positronKineEnergy;
469 positronKineEnergy = 0;
472 positronDirection.
rotateUz(photonDirection);
474 positronDirection, positronKineEnergy);
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"
489 if (positronKineEnergy)
490 G4cout <<
"Positron (not at rest) " << positronKineEnergy/
keV <<
" keV" <<
G4endl;
492 if (localEnergyDeposit)
493 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
494 G4cout <<
"Total final state: " << (electronKineEnergy+positronKineEnergy+
497 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
501 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
503 if (energyDiff > 0.05*
keV)
504 G4cout <<
"Warning from G4PenelopeGammaConversion: problem with energy conservation: "
505 << (electronKineEnergy+positronKineEnergy+
507 <<
" keV (final) vs. " << photonEnergy/
keV <<
" keV (initial)" <<
G4endl;
517 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
522 G4cout <<
"G4PenelopeGammaConversionModel::ReadDataFile()" <<
G4endl;
523 G4cout <<
"Going to read Gamma Conversion data files for Z=" <<
Z <<
G4endl;
526 char* path = std::getenv(
"G4LEDATA");
530 "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
531 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
539 std::ostringstream ost;
541 ost << path <<
"/penelope/pairproduction/pdgpp" <<
Z <<
".p08";
543 ost << path <<
"/penelope/pairproduction/pdgpp0" <<
Z <<
".p08";
544 std::ifstream
file(ost.str().c_str());
547 G4String excep =
"G4PenelopeGammaConversionModel - data file " +
548 G4String(ost.str()) +
" not found!";
549 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
557 while( getline(
file, line) )
563 file.open(ost.str().c_str());
574 ed <<
"Corrupted data file for Z=" <<
Z <<
G4endl;
575 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
581 for (
size_t i=0;i<ndata;i++)
612 zeff = (*elementVector)[0]->GetZ();
620 for (
G4int i=0;i<nElements;i++)
622 G4double Zelement = (*elementVector)[i]->GetZ();
623 G4double Aelement = (*elementVector)[i]->GetAtomicMassAmu();
624 atot += Aelement*fractionVector[i];
625 zeff += Zelement*Aelement*fractionVector[i];
627 atot /=
material->GetTotNbOfAtomsPerVolume();
628 zeff /= (
material->GetTotNbOfAtomsPerVolume()*atot);
630 intZ = (
G4int) (zeff+0.25);
645 G4double fc = alzSquared*(0.202059-alzSquared*
647 (0.00835-alzSquared*(0.00201-alzSquared*
649 (0.00012-alzSquared*0.00003)))))
650 +1.0/(alzSquared+1.0));
658 std::pair<G4double,G4double> myPair(0,0);
669 G4cout <<
"Average Z for material " <<
material->GetName() <<
" = " <<
671 G4cout <<
"Effective radius for material " <<
material->GetName() <<
" = " <<
674 G4cout <<
"Screening parameters F0 for material " <<
material->GetName() <<
" = " <<
675 f0a <<
"," << f0b <<
G4endl;
682std::pair<G4double,G4double>
690 std::pair<G4double,G4double> result(0.,0.);
700 f2 += 2.0*BSquared*(4.0-
a0-3.0*
G4Log((1.0+BSquared)/BSquared));
G4double B(G4double temperature)
static const G4double eps
std::vector< const G4Element * > G4ElementVector
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
static constexpr double twopi
static constexpr double barn
static constexpr double cm2
static constexpr double keV
static constexpr double eV
static constexpr double GeV
static constexpr double MeV
#define G4MUTEX_INITIALIZER
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
const G4Material * GetMaterial() const
const G4String & GetName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
std::pair< G4double, G4double > GetScreeningFunctions(G4double)
void ReadDataFile(const G4int Z)
const G4ParticleDefinition * fParticle
void InitializeScreeningFunctions(const G4Material *)
static G4double fAtomicScreeningRadius[fMaxZ+1]
G4double fIntrinsicHighEnergyLimit
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
std::map< const G4Material *, G4double > * fEffectiveCharge
G4ParticleChangeForGamma * fParticleChange
G4PenelopeGammaConversionModel(const G4ParticleDefinition *p=nullptr, const G4String &processName="PenConversion")
G4double fIntrinsicLowEnergyLimit
void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *) override
std::map< const G4Material *, G4double > * fMaterialInvScreeningRadius
std::map< const G4Material *, std::pair< G4double, G4double > > * fScreeningFunction
void SetParticle(const G4ParticleDefinition *)
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
virtual ~G4PenelopeGammaConversionModel()
static G4PhysicsFreeVector * fLogAtomicCrossSection[fMaxZ+1]
void PutValue(const std::size_t index, const G4double e, const G4double value)
G4double Value(const G4double energy, std::size_t &lastidx) const
static G4Positron * Positron()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
void SetHighEnergyLimit(G4double)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double energy(const ThreeVector &p, const G4double m)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4Mutex PenelopeGammaConversionModelMutex