65 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),
66 isInitialised(false),fAtomDeexcitation(0),
103 G4cout <<
"Calling G4PenelopeComptonModel::Initialise()" <<
G4endl;
110 G4cout <<
"WARNING from G4PenelopeComptonModel " <<
G4endl;
111 G4cout <<
"Atomic de-excitation module is not instantiated, so there will not be ";
113 G4cout <<
"Please make sure this is intended" <<
G4endl;
123 G4cout <<
"Penelope Compton model v2008 is initialized " << G4endl
133 ed <<
"Using the Penelope Compton model outside its intrinsic validity range. "
138 ed <<
"Result of the simulation have to be taken with care" <<
G4endl;
139 G4Exception(
"G4PenelopeComptonModel::Initialise()",
156 G4cout <<
"Calling G4PenelopeComptonModel::InitialiseLocal()" <<
G4endl;
196 G4cout <<
"Calling CrossSectionPerVolume() of G4PenelopeComptonModel" <<
G4endl;
210 size_t numberOfOscillators = theTable->size();
211 for (
size_t i=0;i<numberOfOscillators;i++)
222 cs *=
pi*classic_electr_radius*classic_electr_radius;
231 G4cout <<
"Material " << material->
GetName() <<
" has " << atPerMol <<
232 "atoms per molecule" <<
G4endl;
237 moleculeDensity = atomDensity/atPerMol;
239 G4double csvolume = cs*moleculeDensity;
242 G4cout <<
"Compton mean free path at " << energy/
keV <<
" keV for material " <<
243 material->
GetName() <<
" = " << (1./csvolume)/
mm <<
" mm" << G4endl;
259 G4cout <<
"*** G4PenelopeComptonModel -- WARNING ***" <<
G4endl;
260 G4cout <<
"Penelope Compton model v2008 does not calculate cross section _per atom_ " <<
G4endl;
261 G4cout <<
"so the result is always zero. For physics values, please invoke " <<
G4endl;
262 G4cout <<
"GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" <<
G4endl;
296 G4cout <<
"Calling SampleSecondaries() of G4PenelopeComptonModel" <<
G4endl;
319 size_t numberOfOscillators = theTable->size();
320 size_t targetOscillator = 0;
323 G4double ek = photonEnergy0/electron_mass_c2;
337 if (photonEnergy0 > 5*
MeV)
346 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
349 cosTheta = 1.0 - (1.0-tau)/(ek*tau);
353 targetOscillator = numberOfOscillators-1;
355 G4bool levelFound =
false;
356 for (
size_t j=0;j<numberOfOscillators && !levelFound; j++)
358 S += (*theTable)[j]->GetOscillatorStrength();
361 targetOscillator = j;
366 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
367 }
while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
376 for (
size_t i=0;i<numberOfOscillators;i++)
378 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
379 if (photonEnergy0 > ionEnergy)
381 G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
382 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
383 oscStren = (*theTable)[i]->GetOscillatorStrength();
384 pzomc = hartreeFunc*(aux2-electron_mass_c2*ionEnergy)/
385 (electron_mass_c2*std::sqrt(2.0*aux2+ionEnergy*ionEnergy));
387 rni = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
388 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
390 rni = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
391 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
403 cdt1 = (1.0-tau)/(ek*tau);
406 for (
size_t i=0;i<numberOfOscillators;i++)
408 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
409 if (photonEnergy0 > ionEnergy)
411 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
412 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
413 oscStren = (*theTable)[i]->GetOscillatorStrength();
414 pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/
415 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
417 rn[i] = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
418 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
420 rn[i] = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
421 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
429 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
432 cosTheta = 1.0 - cdt1;
441 targetOscillator = numberOfOscillators-1;
442 G4bool levelFound =
false;
443 for (
size_t i=0;i<numberOfOscillators && !levelFound;i++)
447 targetOscillator = i;
452 hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor();
453 oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength();
455 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
456 (std::sqrt(2.0)*hartreeFunc);
458 pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/
459 (std::sqrt(2.0)*hartreeFunc);
460 }
while (pzomc < -1);
463 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
464 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
477 epsilon = (tau/
b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
479 epsilon = (tau/
b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
483 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
485 G4double dirx = sinTheta * std::cos(phi);
486 G4double diry = sinTheta * std::sin(phi);
491 photonDirection1.rotateUz(photonDirection0);
494 G4double photonEnergy1 = epsilon * photonEnergy0;
496 if (photonEnergy1 > 0.)
506 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
509 photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
513 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
516 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
520 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
521 G4int Z = (
G4int) (*theTable)[targetOscillator]->GetParentZ();
528 if (Z > 0 && shFlag<30)
534 G4double ionEnergyInPenelopeDatabase = ionEnergy;
536 ionEnergy =
std::max(bindingEnergy,ionEnergyInPenelopeDatabase);
540 G4double eKineticEnergy = diffEnergy - ionEnergy;
541 G4double localEnergyDeposit = ionEnergy;
545 if (eKineticEnergy < 0)
551 localEnergyDeposit = diffEnergy;
552 eKineticEnergy = 0.0;
563 size_t nBefore = fvect->size();
565 size_t nAfter = fvect->size();
567 if (nAfter > nBefore)
569 for (
size_t j=nBefore;j<nAfter;j++)
571 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
572 localEnergyDeposit -= itsEnergy;
574 energyInFluorescence += itsEnergy;
576 energyInAuger += itsEnergy;
659 eDirection.rotateUz(photonDirection0);
661 eDirection,eKineticEnergy) ;
662 fvect->push_back(electron);
665 if (localEnergyDeposit < 0)
668 <<
"G4PenelopeComptonModel::SampleSecondaries - Negative energy deposit"
670 localEnergyDeposit=0.;
676 electronEnergy = eKineticEnergy;
679 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
680 G4cout <<
"Energy balance from G4PenelopeCompton" <<
G4endl;
681 G4cout <<
"Incoming photon energy: " << photonEnergy0/
keV <<
" keV" <<
G4endl;
682 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
683 G4cout <<
"Scattered photon: " << photonEnergy1/
keV <<
" keV" <<
G4endl;
684 G4cout <<
"Scattered electron " << electronEnergy/
keV <<
" keV" <<
G4endl;
685 if (energyInFluorescence)
686 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
688 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
689 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
690 G4cout <<
"Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+
691 localEnergyDeposit+energyInAuger)/
keV <<
693 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
697 G4double energyDiff = std::fabs(photonEnergy1+
698 electronEnergy+energyInFluorescence+
699 localEnergyDeposit+energyInAuger-photonEnergy0);
700 if (energyDiff > 0.05*
keV)
701 G4cout <<
"Warning from G4PenelopeCompton: problem with energy conservation: " <<
702 (photonEnergy1+electronEnergy+energyInFluorescence+energyInAuger+localEnergyDeposit)/
keV <<
703 " keV (final) vs. " <<
704 photonEnergy0/
keV <<
" keV (initial)" << G4endl;
727 static const G4double k2 = std::sqrt(2.);
730 if (energy < ionEnergy)
735 G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1;
739 G4double aux = energy*(energy-ionEnergy)*cdt1;
741 (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
745 sia = 1.0-0.5*std::exp(0.5-(k1+k2*x)*(k1+k2*x));
747 sia = 0.5*std::exp(0.5-(k1-k2*x)*(k1-k2*x));
752 if (std::fabs(Pzimax) < pf)
754 G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*cosTheta;
757 (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc
758 *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf));
762 G4double XKN = EOEC+ECOE-1.0+cosTheta*cosTheta;
765 G4double diffCS = ECOE*ECOE*XKN*sia;
782 const G4int npoints=10;
783 const G4int ncallsmax=20000;
785 static const G4double Abscissas[10] = {7.652651133497334e-02,2.2778585114164508e-01,3.7370608871541956e-01,
786 5.1086700195082710e-01,6.3605368072651503e-01,7.4633190646015079e-01,
787 8.3911697182221882e-01,9.1223442825132591e-01,9.6397192727791379e-01,
788 9.9312859918509492e-01};
789 static const G4double Weights[10] = {1.5275338713072585e-01,1.4917298647260375e-01,1.4209610931838205e-01,
790 1.3168863844917663e-01,1.1819453196151842e-01,1.0193011981724044e-01,
791 8.3276741576704749e-02,6.2672048334109064e-02,4.0601429800386941e-02,
792 1.7614007139152118e-02};
807 G4double b=0.5*(HighPoint+LowPoint);
811 for (
G4int i=2;i<=npoints;i++)
817 G4int icall = 2*npoints;
832 for (
G4int i=1;i<=LH;i++){
843 for (
G4int j=1;j<npoints;j++)
846 dLocal += Weights[j]*
856 for (
G4int j=1;j<npoints;j++)
859 dLocal += Weights[j]*
863 icall=icall+4*npoints;
865 if (std::abs(s12-si)<
std::max(Ptol*std::abs(s12),1e-35))
877 if (icall>ncallsmax || LHN>nst)
880 G4cout <<
"LowPoint: " << LowPoint <<
", High Point: " << HighPoint <<
G4endl;
882 G4cout <<
"Calls: " << icall <<
", Integral: " << sumga <<
", Error: " << Err <<
G4endl;
883 G4cout <<
"Number of open subintervals: " << LHN <<
G4endl;
884 G4cout <<
"WARNING: the required accuracy has not been attained" <<
G4endl;
888 Err=std::abs(sumr)/
std::max(std::abs(sumr+sumga),1e-35);
889 if (Err < Ctol || LHN == 0)
892 for (
G4int i=0;i<LH;i++)
897 }
while(Ctol < 1.0 && loopAgain);
915 G4double ek =energy/electron_mass_c2;
921 G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*std::log(t0)-(1.0/t0);
925 for (
size_t i=0;i<theTable->size();i++)
929 G4double tau=(energy-ionEnergy)/energy;
932 G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1*std::log(tau)-(1.0/tau);
935 cs += stre*(csu-csl);
G4double KleinNishinaCrossSection(G4double energy, const G4Material *)
G4double GetIonisationEnergy()
G4double LowEnergyLimit() const
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
const G4String & GetName() const
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
static G4Electron * Definition()
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double BindingEnergy() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
G4double GetOscillatorStrength()
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
const G4AtomicTransitionManager * fTransitionManager
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
G4double fIntrinsicLowEnergyLimit
static G4PenelopeOscillatorManager * GetOscillatorManager()
void SetParticle(const G4ParticleDefinition *)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double GetTotNbOfAtomsPerVolume() const
G4double fIntrinsicHighEnergyLimit
G4ParticleChangeForGamma * fParticleChange
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4ParticleDefinition * fParticle
G4double energy(const ThreeVector &p, const G4double m)
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
const G4double x[NPOINTSGL]
G4double OscillatorTotalCrossSection(G4double energy, G4PenelopeOscillator *osc)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double DifferentialCrossSection(G4double cdt, G4double energy, G4PenelopeOscillator *osc)
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4PenelopeOscillatorManager * oscManager
virtual ~G4PenelopeComptonModel()
G4PenelopeComptonModel(const G4ParticleDefinition *p=0, const G4String &processName="PenCompton")
G4VAtomDeexcitation * AtomDeexcitation()
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
G4ThreeVector G4ParticleMomentum
G4double bindingEnergy(G4int A, G4int Z)
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
G4double GetTotalZ(const G4Material *)
These are cumulative for the molecule Returns the total Z for the molecule.
double epsilon(double density, double temperature)
G4VAtomDeexcitation * fAtomDeexcitation
G4ParticleChangeForGamma * GetParticleChangeForGamma()
static G4Gamma * Definition()
const G4Material * GetMaterial() const
G4double GetHartreeFactor()
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)