58 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),
59 logAtomicCrossSection(0),
60 fEffectiveCharge(0),fMaterialInvScreeningRadius(0),
61 fScreeningFunction(0),isInitialised(false),fLocalTable(false)
95 if (item.second)
delete item.second;
114 G4cout <<
"Calling G4PenelopeGammaConversionModel::Initialise()" <<
G4endl;
157 G4int iZ = (
G4int) theElementVector->at(j)->GetZ();
170 G4cout <<
"Penelope Gamma Conversion model v2008 is initialized " << G4endl
191 G4cout <<
"Calling G4PenelopeGammaConversionModel::InitialiseLocal()" <<
G4endl;
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()",
261 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
275 G4cout <<
"Gamma conversion cross section at " << energy/
MeV <<
" MeV for Z=" << Z <<
302 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" <<
G4endl;
339 ed <<
"Unable to allocate the EffectiveCharge data for " <<
341 ed <<
"This can happen only in Unit Tests" <<
G4endl;
342 G4Exception(
"G4PenelopeGammaConversionModel::SampleSecondaries()",
346 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
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;
390 eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
392 eps = 0.5+xr*std::pow(ru2m1,1./3.);
416 G4double electronTotEnergy = eps*photonEnergy;
417 G4double positronTotEnergy = (1.0-
eps)*photonEnergy;
425 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
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);
434 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*
electron_mass_c2));
435 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
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);
446 if (electronKineEnergy > 0.0)
448 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
449 electronDirection.
rotateUz(photonDirection);
453 fvect->push_back(electron);
457 localEnergyDeposit += electronKineEnergy;
458 electronKineEnergy = 0;
463 if (positronKineEnergy < 0.0)
465 localEnergyDeposit += positronKineEnergy;
466 positronKineEnergy = 0;
469 positronDirection.
rotateUz(photonDirection);
471 positronDirection, positronKineEnergy);
472 fvect->push_back(positron);
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"
486 if (positronKineEnergy)
487 G4cout <<
"Positron (not at rest) " << positronKineEnergy/
keV <<
" keV" <<
G4endl;
489 if (localEnergyDeposit)
490 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
491 G4cout <<
"Total final state: " << (electronKineEnergy+positronKineEnergy+
494 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
498 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
500 if (energyDiff > 0.05*
keV)
501 G4cout <<
"Warning from G4PenelopeGammaConversion: problem with energy conservation: "
502 << (electronKineEnergy+positronKineEnergy+
504 <<
" keV (final) vs. " << photonEnergy/
keV <<
" keV (initial)" << G4endl;
514 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
519 G4cout <<
"G4PenelopeGammaConversionModel::ReadDataFile()" <<
G4endl;
520 G4cout <<
"Going to read Gamma Conversion data files for Z=" << Z <<
G4endl;
523 char* path = getenv(
"G4LEDATA");
527 "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
528 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
536 std::ostringstream ost;
538 ost << path <<
"/penelope/pairproduction/pdgpp" << Z <<
".p08";
540 ost << path <<
"/penelope/pairproduction/pdgpp0" << Z <<
".p08";
541 std::ifstream
file(ost.str().c_str());
544 G4String excep =
"G4PenelopeGammaConversionModel - data file " +
545 G4String(ost.str()) +
" not found!";
546 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
554 while( getline(
file, line) )
561 file.open(ost.str().c_str());
572 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
573 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
579 for (
size_t i=0;i<ndata;i++)
587 theVec->
PutValue(i,std::log(ene),std::log(xs));
594 ed <<
"Problem with allocation of logAtomicCrossSection data table " <<
G4endl;
595 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
610 G4double temp[99] = {1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,6.4696e+01,
611 6.1228e+01,5.7524e+01,5.4033e+01,5.0787e+01,4.7851e+01,4.6373e+01,
612 4.5401e+01,4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,4.1586e+01,
613 4.0953e+01,4.0524e+01,4.0256e+01,3.9756e+01,3.9144e+01,3.8462e+01,
614 3.7778e+01,3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,3.4688e+01,
615 3.4197e+01,3.3786e+01,3.3422e+01,3.3068e+01,3.2740e+01,3.2438e+01,
616 3.2143e+01,3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,3.0950e+01,
617 3.0758e+01,3.0561e+01,3.0285e+01,3.0097e+01,2.9832e+01,2.9581e+01,
618 2.9411e+01,2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,2.8580e+01,
619 2.8442e+01,2.8312e+01,2.8139e+01,2.7973e+01,2.7819e+01,2.7675e+01,
620 2.7496e+01,2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,2.6516e+01,
621 2.6304e+01,2.6108e+01,2.5929e+01,2.5730e+01,2.5577e+01,2.5403e+01,
622 2.5245e+01,2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,2.4506e+01,
623 2.4391e+01,2.4262e+01,2.4145e+01,2.4039e+01,2.3922e+01,2.3813e+01,
624 2.3712e+01,2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,2.3238e+01,
625 2.3139e+01,2.3048e+01,2.2967e+01,2.2833e+01,2.2694e+01,2.2624e+01,
626 2.2545e+01,2.2446e+01,2.2358e+01,2.2264e+01};
629 for (
G4int i=0;i<99;i++)
656 zeff = (*elementVector)[0]->GetZ();
664 for (
G4int i=0;i<nElements;i++)
666 G4double Zelement = (*elementVector)[i]->GetZ();
667 G4double Aelement = (*elementVector)[i]->GetAtomicMassAmu();
668 atot += Aelement*fractionVector[i];
669 zeff += Zelement*Aelement*fractionVector[i];
674 intZ = (
G4int) (zeff+0.25);
689 G4double fc = alzSquared*(0.202059-alzSquared*
691 (0.00835-alzSquared*(0.00201-alzSquared*
693 (0.00012-alzSquared*0.00003)))))
694 +1.0/(alzSquared+1.0));
702 std::pair<G4double,G4double> myPair(0,0);
713 G4cout <<
"Average Z for material " << material->
GetName() <<
" = " <<
715 G4cout <<
"Effective radius for material " << material->
GetName() <<
" = " <<
718 G4cout <<
"Screening parameters F0 for material " << material->
GetName() <<
" = " <<
719 f0a <<
"," << f0b <<
G4endl;
726 std::pair<G4double,G4double>
734 std::pair<G4double,G4double>
result(0.,0.);
744 f2 += 2.0*BSquared*(4.0-a0-3.0*std::log((1.0+BSquared)/BSquared));
virtual ~G4PenelopeGammaConversionModel()
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double fIntrinsicLowEnergyLimit
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
std::ostringstream G4ExceptionDescription
const G4ParticleDefinition * fParticle
static constexpr double MeV
void SetHighEnergyLimit(G4double)
static constexpr double keV
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
const G4ThreeVector & GetMomentumDirection() const
const G4ElementVector * GetElementVector() const
G4ParticleChangeForGamma * fParticleChange
void InitializeScreeningFunctions(const G4Material *)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
std::pair< G4double, G4double > GetScreeningFunctions(G4double)
Hep3Vector & rotateUz(const Hep3Vector &)
const G4double * GetVecNbOfAtomsPerVolume() const
G4double Value(G4double theEnergy, size_t &lastidx) const
std::map< G4int, G4PhysicsFreeVector * > * logAtomicCrossSection
G4double LowEnergyLimit() const
const G4String & GetName() const
#define G4MUTEX_INITIALIZER
std::map< const G4Material *, std::pair< G4double, G4double > > * fScreeningFunction
static constexpr double cm2
void ReadDataFile(const G4int Z)
static constexpr double electron_mass_c2
G4double fIntrinsicHighEnergyLimit
static constexpr double twopi
void SetParticle(const G4ParticleDefinition *)
void InitializeScreeningRadii()
static G4Positron * Positron()
static constexpr double eV
static G4Electron * Electron()
size_t GetTableSize() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double G4ParticleHPJENDLHEData::G4double result
static G4ProductionCutsTable * GetProductionCutsTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::vector< G4Element * > G4ElementVector
std::map< const G4Material *, G4double > * fEffectiveCharge
G4double fAtomicScreeningRadius[99]
static constexpr double barn
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
std::map< const G4Material *, G4double > * fMaterialInvScreeningRadius
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
const G4Material * GetMaterial() const
void PutValue(size_t index, G4double energy, G4double dataValue)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double GetTotNbOfAtomsPerVolume() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static const G4double eps
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *)
static constexpr double fine_structure_const
G4double HighEnergyLimit() const
double B(double temperature)
void ProposeTrackStatus(G4TrackStatus status)
static constexpr double GeV
G4PenelopeGammaConversionModel(const G4ParticleDefinition *p=0, const G4String &processName="PenConversion")
size_t GetNumberOfElements() const