57 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),
58 logAtomicCrossSection(0),
59 fEffectiveCharge(0),fMaterialInvScreeningRadius(0),
60 fScreeningFunction(0),isInitialised(false),fLocalTable(false)
64 fIntrinsicHighEnergyLimit = 100.0*
GeV;
65 fSmallEnergy = 1.1*
MeV;
67 InitializeScreeningRadii();
91 if (logAtomicCrossSection)
98 delete logAtomicCrossSection;
100 if (fEffectiveCharge)
101 delete fEffectiveCharge;
102 if (fMaterialInvScreeningRadius)
103 delete fMaterialInvScreeningRadius;
104 if (fScreeningFunction)
105 delete fScreeningFunction;
115 if (verboseLevel > 3)
116 G4cout <<
"Calling G4PenelopeGammaConversionModel::Initialise()" <<
G4endl;
124 if (!logAtomicCrossSection)
125 logAtomicCrossSection =
new std::map<G4int,G4PhysicsFreeVector*>;
128 if (fEffectiveCharge)
130 delete fEffectiveCharge;
131 fEffectiveCharge = 0;
133 if (fMaterialInvScreeningRadius)
135 delete fMaterialInvScreeningRadius;
136 fMaterialInvScreeningRadius = 0;
138 if (fScreeningFunction)
140 delete fScreeningFunction;
141 fScreeningFunction = 0;
144 fEffectiveCharge =
new std::map<const G4Material*,G4double>;
145 fMaterialInvScreeningRadius =
new std::map<const G4Material*,G4double>;
146 fScreeningFunction =
new std::map<const G4Material*,std::pair<G4double,G4double> >;
159 G4int iZ = (
G4int) theElementVector->at(j)->GetZ();
161 if (!logAtomicCrossSection->count(iZ))
166 if (!fEffectiveCharge->count(material))
167 InitializeScreeningFunctions(material);
171 if (verboseLevel > 0) {
172 G4cout <<
"Penelope Gamma Conversion model v2008 is initialized " << G4endl
182 if(isInitialised)
return;
184 isInitialised =
true;
192 if (verboseLevel > 3)
193 G4cout <<
"Calling G4PenelopeGammaConversionModel::InitialiseLocal()" <<
G4endl;
206 fEffectiveCharge = theModel->fEffectiveCharge;
207 fMaterialInvScreeningRadius = theModel->fMaterialInvScreeningRadius;
208 fScreeningFunction = theModel->fScreeningFunction;
209 logAtomicCrossSection = theModel->logAtomicCrossSection;
212 verboseLevel = theModel->verboseLevel;
234 if (energy < fIntrinsicLowEnergyLimit)
241 if (!logAtomicCrossSection)
246 logAtomicCrossSection =
new std::map<G4int,G4PhysicsFreeVector*>;
249 if (!logAtomicCrossSection->count(iZ))
253 if (verboseLevel > 0)
257 ed <<
"Unable to retrieve the cross section table for Z=" << iZ <<
G4endl;
258 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
259 G4Exception(
"G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom()",
263 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
274 cs = std::exp(logXS);
276 if (verboseLevel > 2)
277 G4cout <<
"Gamma conversion cross section at " << energy/
MeV <<
" MeV for Z=" << Z <<
303 if (verboseLevel > 3)
304 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" <<
G4endl;
312 if (photonEnergy <= fIntrinsicLowEnergyLimit)
323 if (!fEffectiveCharge)
328 fEffectiveCharge =
new std::map<const G4Material*,G4double>;
329 fMaterialInvScreeningRadius =
new std::map<const G4Material*,G4double>;
330 fScreeningFunction =
new std::map<const G4Material*,std::pair<G4double,G4double> >;
333 if (!fEffectiveCharge->count(mat))
337 if (verboseLevel > 0)
341 ed <<
"Unable to allocate the EffectiveCharge data for " <<
343 ed <<
"This can happen only in Unit Tests" <<
G4endl;
344 G4Exception(
"G4PenelopeGammaConversionModel::SampleSecondaries()",
348 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
349 InitializeScreeningFunctions(mat);
358 if (photonEnergy < fSmallEnergy)
363 G4double effC = fEffectiveCharge->find(mat)->second;
366 G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
367 +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
368 -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
369 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
371 G4double F0b = fScreeningFunction->find(mat)->second.second;
373 G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
375 std::pair<G4double,G4double> scree = GetScreeningFunctions(bmin);
392 eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
394 eps = 0.5+xr*std::pow(ru2m1,1./3.);
395 G4double B = eki/(invRad*eps*(1.0-eps));
396 scree = GetScreeningFunctions(B);
405 G4double B = eki/(invRad*eps*(1.0-eps));
406 scree = GetScreeningFunctions(B);
415 if (verboseLevel > 4)
418 G4double electronTotEnergy = eps*photonEnergy;
419 G4double positronTotEnergy = (1.0-eps)*photonEnergy;
427 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
429 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
430 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
436 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*
electron_mass_c2));
437 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
439 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
440 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
448 if (electronKineEnergy > 0.0)
450 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
451 electronDirection.
rotateUz(photonDirection);
455 fvect->push_back(electron);
459 localEnergyDeposit += electronKineEnergy;
460 electronKineEnergy = 0;
465 if (positronKineEnergy < 0.0)
467 localEnergyDeposit += positronKineEnergy;
468 positronKineEnergy = 0;
471 positronDirection.
rotateUz(photonDirection);
473 positronDirection, positronKineEnergy);
474 fvect->push_back(positron);
479 if (verboseLevel > 1)
481 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
482 G4cout <<
"Energy balance from G4PenelopeGammaConversion" <<
G4endl;
483 G4cout <<
"Incoming photon energy: " << photonEnergy/
keV <<
" keV" <<
G4endl;
484 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
485 if (electronKineEnergy)
486 G4cout <<
"Electron (explicitely produced) " << electronKineEnergy/
keV <<
" keV"
488 if (positronKineEnergy)
489 G4cout <<
"Positron (not at rest) " << positronKineEnergy/
keV <<
" keV" <<
G4endl;
491 if (localEnergyDeposit)
492 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
493 G4cout <<
"Total final state: " << (electronKineEnergy+positronKineEnergy+
496 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
498 if (verboseLevel > 0)
500 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
502 if (energyDiff > 0.05*
keV)
503 G4cout <<
"Warning from G4PenelopeGammaConversion: problem with energy conservation: "
504 << (electronKineEnergy+positronKineEnergy+
506 <<
" keV (final) vs. " << photonEnergy/
keV <<
" keV (initial)" << G4endl;
512 void G4PenelopeGammaConversionModel::ReadDataFile(
const G4int Z)
516 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
519 if (verboseLevel > 2)
521 G4cout <<
"G4PenelopeGammaConversionModel::ReadDataFile()" <<
G4endl;
522 G4cout <<
"Going to read Gamma Conversion data files for Z=" << Z <<
G4endl;
525 char* path = getenv(
"G4LEDATA");
529 "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
530 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
538 std::ostringstream ost;
540 ost << path <<
"/penelope/pairproduction/pdgpp" << Z <<
".p08";
542 ost << path <<
"/penelope/pairproduction/pdgpp0" << Z <<
".p08";
543 std::ifstream file(ost.str().c_str());
546 G4String excep =
"G4PenelopeGammaConversionModel - data file " +
547 G4String(ost.str()) +
" not found!";
548 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
556 while( getline(file, line) )
563 file.open(ost.str().c_str());
567 if (verboseLevel > 3)
574 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
575 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
581 for (
size_t i=0;i<ndata;i++)
589 theVec->
PutValue(i,std::log(ene),std::log(xs));
593 if (!logAtomicCrossSection)
596 ed <<
"Problem with allocation of logAtomicCrossSection data table " <<
G4endl;
597 G4Exception(
"G4PenelopeGammaConversionModel::ReadDataFile()",
602 logAtomicCrossSection->insert(std::make_pair(Z,theVec));
610 void G4PenelopeGammaConversionModel::InitializeScreeningRadii()
612 G4double temp[99] = {1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,6.4696e+01,
613 6.1228e+01,5.7524e+01,5.4033e+01,5.0787e+01,4.7851e+01,4.6373e+01,
614 4.5401e+01,4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,4.1586e+01,
615 4.0953e+01,4.0524e+01,4.0256e+01,3.9756e+01,3.9144e+01,3.8462e+01,
616 3.7778e+01,3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,3.4688e+01,
617 3.4197e+01,3.3786e+01,3.3422e+01,3.3068e+01,3.2740e+01,3.2438e+01,
618 3.2143e+01,3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,3.0950e+01,
619 3.0758e+01,3.0561e+01,3.0285e+01,3.0097e+01,2.9832e+01,2.9581e+01,
620 2.9411e+01,2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,2.8580e+01,
621 2.8442e+01,2.8312e+01,2.8139e+01,2.7973e+01,2.7819e+01,2.7675e+01,
622 2.7496e+01,2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,2.6516e+01,
623 2.6304e+01,2.6108e+01,2.5929e+01,2.5730e+01,2.5577e+01,2.5403e+01,
624 2.5245e+01,2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,2.4506e+01,
625 2.4391e+01,2.4262e+01,2.4145e+01,2.4039e+01,2.3922e+01,2.3813e+01,
626 2.3712e+01,2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,2.3238e+01,
627 2.3139e+01,2.3048e+01,2.2967e+01,2.2833e+01,2.2694e+01,2.2624e+01,
628 2.2545e+01,2.2446e+01,2.2358e+01,2.2264e+01};
631 for (
G4int i=0;i<99;i++)
632 fAtomicScreeningRadius[i] = temp[i];
637 void G4PenelopeGammaConversionModel::InitializeScreeningFunctions(
const G4Material*
material)
658 zeff = (*elementVector)[0]->GetZ();
666 for (
G4int i=0;i<nElements;i++)
668 G4double Zelement = (*elementVector)[i]->GetZ();
669 G4double Aelement = (*elementVector)[i]->GetAtomicMassAmu();
670 atot += Aelement*fractionVector[i];
671 zeff += Zelement*Aelement*fractionVector[i];
676 intZ = (
G4int) (zeff+0.25);
683 if (fEffectiveCharge)
684 fEffectiveCharge->insert(std::make_pair(material,zeff));
691 G4double fc = alzSquared*(0.202059-alzSquared*
693 (0.00835-alzSquared*(0.00201-alzSquared*
695 (0.00012-alzSquared*0.00003)))))
696 +1.0/(alzSquared+1.0));
700 G4double matRadius = 2.0/ fAtomicScreeningRadius[intZ-1];
701 if (fMaterialInvScreeningRadius)
702 fMaterialInvScreeningRadius->insert(std::make_pair(material,matRadius));
704 std::pair<G4double,G4double> myPair(0,0);
705 G4double f0a = 4.0*std::log(fAtomicScreeningRadius[intZ-1]);
710 if (fScreeningFunction)
711 fScreeningFunction->insert(std::make_pair(material,myPair));
713 if (verboseLevel > 2)
715 G4cout <<
"Average Z for material " << material->
GetName() <<
" = " <<
717 G4cout <<
"Effective radius for material " << material->
GetName() <<
" = " <<
718 fAtomicScreeningRadius[intZ-1] <<
" m_e*c/hbar --> BCB = " <<
720 G4cout <<
"Screening parameters F0 for material " << material->
GetName() <<
" = " <<
721 f0a <<
"," << f0b <<
G4endl;
728 std::pair<G4double,G4double>
729 G4PenelopeGammaConversionModel::GetScreeningFunctions(
G4double B)
736 std::pair<G4double,G4double> result(0.,0.);
738 G4double f1 = 2.0-2.0*std::log(1.0+BSquared);
744 G4double a0 = 4.0*B*std::atan(1./B);
746 f2 += 2.0*BSquared*(4.0-a0-3.0*std::log((1.0+BSquared)/BSquared));
G4double LowEnergyLimit() const
std::vector< G4Element * > G4ElementVector
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
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)
const G4String & GetName() const
G4PenelopeGammaConversionModel(const G4ParticleDefinition *p=0, const G4String &processName="PenConversion")
const G4ElementVector * GetElementVector() const
#define G4MUTEX_INITIALIZER
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
const G4double * GetVecNbOfAtomsPerVolume() const
double precision function energy(A, Z)
G4GLOB_DLL std::ostream G4cout
G4ParticleChangeForGamma * fParticleChange
size_t GetTableSize() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
const G4ParticleDefinition * fParticle
G4double GetTotNbOfAtomsPerVolume() const
static G4ProductionCutsTable * GetProductionCutsTable()
static G4Positron * Positron()
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
size_t GetNumberOfElements() const
void ProposeTrackStatus(G4TrackStatus status)
virtual ~G4PenelopeGammaConversionModel()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const