66 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),
67 isInitialised(false),fAtomDeexcitation(0),
104 G4cout <<
"Calling G4PenelopeComptonModel::Initialise()" <<
G4endl;
111 G4cout <<
"WARNING from G4PenelopeComptonModel " <<
G4endl;
112 G4cout <<
"Atomic de-excitation module is not instantiated, so there will not be ";
114 G4cout <<
"Please make sure this is intended" <<
G4endl;
124 G4cout <<
"Penelope Compton model v2008 is initialized " << G4endl
134 ed <<
"Using the Penelope Compton model outside its intrinsic validity range. "
139 ed <<
"Result of the simulation have to be taken with care" <<
G4endl;
140 G4Exception(
"G4PenelopeComptonModel::Initialise()",
157 G4cout <<
"Calling G4PenelopeComptonModel::InitialiseLocal()" <<
G4endl;
197 G4cout <<
"Calling CrossSectionPerVolume() of G4PenelopeComptonModel" <<
G4endl;
211 size_t numberOfOscillators = theTable->size();
212 for (
size_t i=0;i<numberOfOscillators;i++)
232 G4cout <<
"Material " << material->
GetName() <<
" has " << atPerMol <<
233 "atoms per molecule" <<
G4endl;
238 moleculeDensity = atomDensity/atPerMol;
240 G4double csvolume = cs*moleculeDensity;
243 G4cout <<
"Compton mean free path at " << energy/
keV <<
" keV for material " <<
244 material->
GetName() <<
" = " << (1./csvolume)/
mm <<
" mm" << G4endl;
260 G4cout <<
"*** G4PenelopeComptonModel -- WARNING ***" <<
G4endl;
261 G4cout <<
"Penelope Compton model v2008 does not calculate cross section _per atom_ " <<
G4endl;
262 G4cout <<
"so the result is always zero. For physics values, please invoke " <<
G4endl;
263 G4cout <<
"GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" <<
G4endl;
297 G4cout <<
"Calling SampleSecondaries() of G4PenelopeComptonModel" <<
G4endl;
320 size_t numberOfOscillators = theTable->size();
321 size_t targetOscillator = 0;
338 if (photonEnergy0 > 5*
MeV)
347 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
350 cosTheta = 1.0 - (1.0-tau)/(ek*tau);
354 targetOscillator = numberOfOscillators-1;
356 G4bool levelFound =
false;
357 for (
size_t j=0;j<numberOfOscillators && !levelFound; j++)
359 S += (*theTable)[j]->GetOscillatorStrength();
362 targetOscillator = j;
367 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
368 }
while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
377 for (
size_t i=0;i<numberOfOscillators;i++)
379 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
380 if (photonEnergy0 > ionEnergy)
382 G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
383 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
384 oscStren = (*theTable)[i]->GetOscillatorStrength();
388 rni = 1.0-0.5*
G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
389 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
391 rni = 0.5*
G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
392 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
404 cdt1 = (1.0-tau)/(ek*tau);
407 for (
size_t i=0;i<numberOfOscillators;i++)
409 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
410 if (photonEnergy0 > ionEnergy)
412 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
413 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
414 oscStren = (*theTable)[i]->GetOscillatorStrength();
418 rn[i] = 1.0-0.5*
G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
419 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
421 rn[i] = 0.5*
G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
422 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
430 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
433 cosTheta = 1.0 - cdt1;
442 targetOscillator = numberOfOscillators-1;
443 G4bool levelFound =
false;
444 for (
size_t i=0;i<numberOfOscillators && !levelFound;i++)
448 targetOscillator = i;
453 hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor();
454 oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength();
456 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
457 (std::sqrt(2.0)*hartreeFunc);
459 pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/
460 (std::sqrt(2.0)*hartreeFunc);
461 }
while (pzomc < -1);
464 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
465 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
478 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
480 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
484 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
486 G4double dirx = sinTheta * std::cos(phi);
487 G4double diry = sinTheta * std::sin(phi);
492 photonDirection1.
rotateUz(photonDirection0);
495 G4double photonEnergy1 = epsilon * photonEnergy0;
497 if (photonEnergy1 > 0.)
507 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
510 photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
514 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
517 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
521 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
522 G4int Z = (
G4int) (*theTable)[targetOscillator]->GetParentZ();
529 if (Z > 0 && shFlag<30)
535 G4double ionEnergyInPenelopeDatabase = ionEnergy;
537 ionEnergy =
std::max(bindingEnergy,ionEnergyInPenelopeDatabase);
541 G4double eKineticEnergy = diffEnergy - ionEnergy;
542 G4double localEnergyDeposit = ionEnergy;
546 if (eKineticEnergy < 0)
552 localEnergyDeposit = diffEnergy;
553 eKineticEnergy = 0.0;
564 size_t nBefore = fvect->size();
566 size_t nAfter = fvect->size();
568 if (nAfter > nBefore)
570 for (
size_t j=nBefore;j<nAfter;j++)
572 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
573 localEnergyDeposit -= itsEnergy;
575 energyInFluorescence += itsEnergy;
577 energyInAuger += itsEnergy;
660 eDirection.
rotateUz(photonDirection0);
662 eDirection,eKineticEnergy) ;
663 fvect->push_back(electron);
666 if (localEnergyDeposit < 0)
669 <<
"G4PenelopeComptonModel::SampleSecondaries - Negative energy deposit"
671 localEnergyDeposit=0.;
677 electronEnergy = eKineticEnergy;
680 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
681 G4cout <<
"Energy balance from G4PenelopeCompton" <<
G4endl;
682 G4cout <<
"Incoming photon energy: " << photonEnergy0/
keV <<
" keV" <<
G4endl;
683 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
684 G4cout <<
"Scattered photon: " << photonEnergy1/
keV <<
" keV" <<
G4endl;
685 G4cout <<
"Scattered electron " << electronEnergy/
keV <<
" keV" <<
G4endl;
686 if (energyInFluorescence)
687 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
689 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
690 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
691 G4cout <<
"Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+
692 localEnergyDeposit+energyInAuger)/
keV <<
694 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
698 G4double energyDiff = std::fabs(photonEnergy1+
699 electronEnergy+energyInFluorescence+
700 localEnergyDeposit+energyInAuger-photonEnergy0);
701 if (energyDiff > 0.05*
keV)
702 G4cout <<
"Warning from G4PenelopeCompton: problem with energy conservation: " <<
703 (photonEnergy1+electronEnergy+energyInFluorescence+energyInAuger+localEnergyDeposit)/
keV <<
704 " keV (final) vs. " <<
705 photonEnergy0/
keV <<
" keV (initial)" << G4endl;
728 static const G4double k2 = std::sqrt(2.);
731 if (energy < ionEnergy)
740 G4double aux = energy*(energy-ionEnergy)*cdt1;
746 sia = 1.0-0.5*
G4Exp(0.5-(k1+k2*x)*(k1+k2*x));
748 sia = 0.5*
G4Exp(0.5-(k1-k2*x)*(k1-k2*x));
753 if (std::fabs(Pzimax) < pf)
755 G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*cosTheta;
758 (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc
759 *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf));
763 G4double XKN = EOEC+ECOE-1.0+cosTheta*cosTheta;
766 G4double diffCS = ECOE*ECOE*XKN*sia;
783 const G4int npoints=10;
784 const G4int ncallsmax=20000;
786 static const G4double Abscissas[10] = {7.652651133497334e-02,2.2778585114164508e-01,3.7370608871541956e-01,
787 5.1086700195082710e-01,6.3605368072651503e-01,7.4633190646015079e-01,
788 8.3911697182221882e-01,9.1223442825132591e-01,9.6397192727791379e-01,
789 9.9312859918509492e-01};
790 static const G4double Weights[10] = {1.5275338713072585e-01,1.4917298647260375e-01,1.4209610931838205e-01,
791 1.3168863844917663e-01,1.1819453196151842e-01,1.0193011981724044e-01,
792 8.3276741576704749e-02,6.2672048334109064e-02,4.0601429800386941e-02,
793 1.7614007139152118e-02};
808 G4double b=0.5*(HighPoint+LowPoint);
812 for (
G4int i=2;i<=npoints;i++)
818 G4int icall = 2*npoints;
833 for (
G4int i=1;i<=LH;i++){
844 for (
G4int j=1;j<npoints;j++)
847 dLocal += Weights[j]*
857 for (
G4int j=1;j<npoints;j++)
860 dLocal += Weights[j]*
864 icall=icall+4*npoints;
866 if (std::abs(s12-si)<
std::max(Ptol*std::abs(s12),1
e-35))
878 if (icall>ncallsmax || LHN>nst)
881 G4cout <<
"LowPoint: " << LowPoint <<
", High Point: " << HighPoint <<
G4endl;
883 G4cout <<
"Calls: " << icall <<
", Integral: " << sumga <<
", Error: " << Err <<
G4endl;
884 G4cout <<
"Number of open subintervals: " << LHN <<
G4endl;
885 G4cout <<
"WARNING: the required accuracy has not been attained" <<
G4endl;
889 Err=std::abs(sumr)/
std::max(std::abs(sumr+sumga),1
e-35);
890 if (Err < Ctol || LHN == 0)
893 for (
G4int i=0;i<LH;i++)
898 }
while(Ctol < 1.0 && loopAgain);
922 G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*std::log(t0)-(1.0/t0);
926 for (
size_t i=0;i<theTable->size();i++)
930 G4double tau=(energy-ionEnergy)/energy;
933 G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1*std::log(tau)-(1.0/tau);
936 cs += stre*(csu-csl);
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4Gamma * Definition()
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
G4ParticleChangeForGamma * fParticleChange
G4VAtomDeexcitation * AtomDeexcitation()
static constexpr double MeV
void SetHighEnergyLimit(G4double)
G4double fIntrinsicLowEnergyLimit
G4double DifferentialCrossSection(G4double cdt, G4double energy, G4PenelopeOscillator *osc)
static constexpr double keV
static constexpr double mm
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4ThreeVector & GetMomentumDirection() const
G4double GetIonisationEnergy()
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Hep3Vector & rotateUz(const Hep3Vector &)
static constexpr double classic_electr_radius
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
G4double LowEnergyLimit() const
const G4String & GetName() const
G4VAtomDeexcitation * fAtomDeexcitation
static G4PenelopeOscillatorManager * GetOscillatorManager()
const G4AtomicTransitionManager * fTransitionManager
void SetDeexcitationFlag(G4bool val)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
static constexpr double electron_mass_c2
G4double GetOscillatorStrength()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
static constexpr double twopi
G4double GetHartreeFactor()
double A(double temperature)
G4double OscillatorTotalCrossSection(G4double energy, G4PenelopeOscillator *osc)
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
G4PenelopeComptonModel(const G4ParticleDefinition *p=0, const G4String &processName="PenCompton")
static constexpr double eV
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
static G4Electron * Electron()
static G4Electron * Definition()
double epsilon(double density, double temperature)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double fIntrinsicHighEnergyLimit
virtual ~G4PenelopeComptonModel()
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
static G4LossTableManager * Instance()
const G4Material * GetMaterial() const
static constexpr double pi
void SetParticle(const G4ParticleDefinition *)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
G4double GetTotNbOfAtomsPerVolume() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double GetTotalZ(const G4Material *)
static G4AtomicTransitionManager * Instance()
G4PenelopeOscillatorManager * oscManager
G4double bindingEnergy(G4int A, G4int Z)
G4double HighEnergyLimit() const
void ProposeTrackStatus(G4TrackStatus status)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static constexpr double GeV
G4double KleinNishinaCrossSection(G4double energy, const G4Material *)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
const G4ParticleDefinition * fParticle
G4double BindingEnergy() const