Geant4  10.02.p03
Public Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
G4PenelopeComptonModel Class Reference

#include <G4PenelopeComptonModel.hh>

Inheritance diagram for G4PenelopeComptonModel:
Collaboration diagram for G4PenelopeComptonModel:

Public Member Functions

 G4PenelopeComptonModel (const G4ParticleDefinition *p=0, const G4String &processName="PenCompton")
 
virtual ~G4PenelopeComptonModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle *> *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
void SetVerbosityLevel (G4int lev)
 
G4int GetVerbosityLevel ()
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int, const G4ParticleDefinition *, G4double)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector *> *)
 
G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy)
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4Element * GetCurrentElement () const
 
const G4Isotope * GetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Attributes

G4ParticleChangeForGamma * fParticleChange
 
const G4ParticleDefinition * fParticle
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 

Private Member Functions

void SetParticle (const G4ParticleDefinition *)
 
G4double DifferentialCrossSection (G4double cdt, G4double energy, G4PenelopeOscillator *osc)
 
G4double OscillatorTotalCrossSection (G4double energy, G4PenelopeOscillator *osc)
 
G4double KleinNishinaCrossSection (G4double energy, const G4Material *)
 
G4PenelopeComptonModel & operator= (const G4PenelopeComptonModel &right)
 
 G4PenelopeComptonModel (const G4PenelopeComptonModel &)
 

Private Attributes

G4double fIntrinsicLowEnergyLimit
 
G4double fIntrinsicHighEnergyLimit
 
G4int verboseLevel
 
G4bool isInitialised
 
G4VAtomDeexcitation * fAtomDeexcitation
 
const G4AtomicTransitionManager * fTransitionManager
 
G4PenelopeOscillatorManager * oscManager
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 63 of file G4PenelopeComptonModel.hh.

Constructor & Destructor Documentation

◆ G4PenelopeComptonModel() [1/2]

G4PenelopeComptonModel::G4PenelopeComptonModel ( const G4ParticleDefinition *  p = 0,
const G4String &  processName = "PenCompton" 
)

Definition at line 64 of file G4PenelopeComptonModel.cc.

66  :G4VEmModel(nam),fParticleChange(0),fParticle(0),
67  isInitialised(false),fAtomDeexcitation(0),
68  oscManager(0)
69 {
70  fIntrinsicLowEnergyLimit = 100.0*eV;
71  fIntrinsicHighEnergyLimit = 100.0*GeV;
72  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
73  //
74  oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
75 
76  if (part)
77  SetParticle(part);
78 
79  verboseLevel= 0;
80  // Verbosity scale:
81  // 0 = nothing
82  // 1 = warning for energy non-conservation
83  // 2 = details of energy budget
84  // 3 = calculation of cross sections, file openings, sampling of atoms
85  // 4 = entering in methods
86 
87  //Mark this model as "applicable" for atomic deexcitation
88  SetDeexcitationFlag(true);
89 
90  fTransitionManager = G4AtomicTransitionManager::Instance();
91 }
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:69
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
G4PenelopeComptonModel::fTransitionManager
const G4AtomicTransitionManager * fTransitionManager
Definition: G4PenelopeComptonModel.hh:132
part
TString part[npart]
Definition: Hadr00/scripts/Style.C:32
GeV
static const double GeV
Definition: G4SIunits.hh:214
G4PenelopeOscillatorManager::GetOscillatorManager
static G4PenelopeOscillatorManager * GetOscillatorManager()
Definition: G4PenelopeOscillatorManager.cc:90
G4PenelopeComptonModel::SetParticle
void SetParticle(const G4ParticleDefinition *)
Definition: G4PenelopeComptonModel.cc:948
eV
static const double eV
Definition: G4SIunits.hh:212
G4PenelopeComptonModel::fParticleChange
G4ParticleChangeForGamma * fParticleChange
Definition: G4PenelopeComptonModel.hh:100
G4PenelopeComptonModel::fParticle
const G4ParticleDefinition * fParticle
Definition: G4PenelopeComptonModel.hh:105
G4PenelopeComptonModel::oscManager
G4PenelopeOscillatorManager * oscManager
Definition: G4PenelopeComptonModel.hh:134
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:781
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:45
G4PenelopeComptonModel::fAtomDeexcitation
G4VAtomDeexcitation * fAtomDeexcitation
Definition: G4PenelopeComptonModel.hh:131
Here is the call graph for this function:

◆ ~G4PenelopeComptonModel()

G4PenelopeComptonModel::~G4PenelopeComptonModel ( )
virtual

Definition at line 95 of file G4PenelopeComptonModel.cc.

96 {;}

◆ G4PenelopeComptonModel() [2/2]

G4PenelopeComptonModel::G4PenelopeComptonModel ( const G4PenelopeComptonModel &  )
private

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

G4double G4PenelopeComptonModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  ,
G4double  ,
G4double  ,
G4double  ,
G4double  ,
G4double   
)
virtual

Reimplemented from G4VEmModel.

Definition at line 253 of file G4PenelopeComptonModel.cc.

259 {
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;
264  return 0;
265 }
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61

◆ CrossSectionPerVolume()

G4double G4PenelopeComptonModel::CrossSectionPerVolume ( const G4Material *  material,
const G4ParticleDefinition *  p,
G4double  kineticEnergy,
G4double  cutEnergy = 0.0,
G4double  maxEnergy = DBL_MAX 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 179 of file G4PenelopeComptonModel.cc.

184 {
185  // Penelope model v2008 to calculate the Compton scattering cross section:
186  // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
187  //
188  // The cross section for Compton scattering is calculated according to the Klein-Nishina
189  // formula for energy > 5 MeV.
190  // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega,
191  // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection().
192  // The parametrization includes the J(p)
193  // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations
194  // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
195  //
196  if (verboseLevel > 3)
197  G4cout << "Calling CrossSectionPerVolume() of G4PenelopeComptonModel" << G4endl;
198  SetupForMaterial(p, material, energy);
199 
200 
201  G4double cs = 0;
202  //Force null cross-section if below the low-energy edge of the table
203  if (energy < LowEnergyLimit())
204  return cs;
205 
206  //Retrieve the oscillator table for this material
207  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material);
208 
209  if (energy < 5*MeV) //explicit calculation for E < 5 MeV
210  {
211  size_t numberOfOscillators = theTable->size();
212  for (size_t i=0;i<numberOfOscillators;i++)
213  {
214  G4PenelopeOscillator* theOsc = (*theTable)[i];
215  //sum contributions coming from each oscillator
216  cs += OscillatorTotalCrossSection(energy,theOsc);
217  }
218  }
219  else //use Klein-Nishina for E>5 MeV
220  cs = KleinNishinaCrossSection(energy,material);
221 
222  //cross sections are in units of pi*classic_electr_radius^2
223  cs *= pi*classic_electr_radius*classic_electr_radius;
224 
225  //Now, cs is the cross section *per molecule*, let's calculate the
226  //cross section per volume
227 
228  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
229  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
230 
231  if (verboseLevel > 3)
232  G4cout << "Material " << material->GetName() << " has " << atPerMol <<
233  "atoms per molecule" << G4endl;
234 
235  G4double moleculeDensity = 0.;
236 
237  if (atPerMol)
238  moleculeDensity = atomDensity/atPerMol;
239 
240  G4double csvolume = cs*moleculeDensity;
241 
242  if (verboseLevel > 2)
243  G4cout << "Compton mean free path at " << energy/keV << " keV for material " <<
244  material->GetName() << " = " << (1./csvolume)/mm << " mm" << G4endl;
245  return csvolume;
246 }
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4PenelopeComptonModel::KleinNishinaCrossSection
G4double KleinNishinaCrossSection(G4double energy, const G4Material *)
Definition: G4PenelopeComptonModel.cc:908
MeV
static const double MeV
Definition: G4SIunits.hh:211
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
G4VEmModel::SetupForMaterial
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:403
G4PenelopeOscillatorManager::GetOscillatorTableCompton
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:364
G4cout
G4GLOB_DLL std::ostream G4cout
energy
double energy
Definition: plottest35.C:25
python.hepunit.classic_electr_radius
int classic_electr_radius
Definition: hepunit.py:288
pi
static const double pi
Definition: G4SIunits.hh:74
G4PenelopeOscillatorTable
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
Definition: G4PenelopeOscillatorManager.hh:63
G4PenelopeComptonModel::OscillatorTotalCrossSection
G4double OscillatorTotalCrossSection(G4double energy, G4PenelopeOscillator *osc)
Definition: G4PenelopeComptonModel.cc:773
G4endl
#define G4endl
Definition: G4ios.hh:61
G4PenelopeComptonModel::oscManager
G4PenelopeOscillatorManager * oscManager
Definition: G4PenelopeComptonModel.hh:134
keV
static const double keV
Definition: G4SIunits.hh:213
G4double
double G4double
Definition: G4Types.hh:76
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:178
G4PenelopeOscillatorManager::GetAtomsPerMolecule
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
Definition: G4PenelopeOscillatorManager.cc:1196
mm
static const double mm
Definition: G4SIunits.hh:114
Here is the call graph for this function:

◆ DifferentialCrossSection()

G4double G4PenelopeComptonModel::DifferentialCrossSection ( G4double  cdt,
G4double  energy,
G4PenelopeOscillator *  osc 
)
private

Definition at line 711 of file G4PenelopeComptonModel.cc.

713 {
714  //
715  // Penelope model v2008. Single differential cross section *per electron*
716  // for photon Compton scattering by
717  // electrons in the given atomic oscillator, differential in the direction of the
718  // scattering photon. This is in units of pi*classic_electr_radius**2
719  //
720  // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
721  // The parametrization includes the J(p) distribution profiles for the atomic shells,
722  // that are tabulated from Hartree-Fock calculations
723  // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
724  //
725  G4double ionEnergy = osc->GetIonisationEnergy();
726  G4double harFunc = osc->GetHartreeFactor();
727 
728  static const G4double k2 = std::sqrt(2.);
729  static const G4double k1 = 1./k2;
730 
731  if (energy < ionEnergy)
732  return 0;
733 
734  //energy of the Compton line
735  G4double cdt1 = 1.0-cosTheta;
736  G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1;
737  G4double ECOE = 1.0/EOEC;
738 
739  //Incoherent scattering function (analytical profile)
740  G4double aux = energy*(energy-ionEnergy)*cdt1;
741  G4double Pzimax =
742  (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
743  G4double sia = 0.0;
744  G4double x = harFunc*Pzimax;
745  if (x > 0)
746  sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x));
747  else
748  sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x));
749 
750  //1st order correction, integral of Pz times the Compton profile.
751  //Calculated approximately using a free-electron gas profile
752  G4double pf = 3.0/(4.0*harFunc);
753  if (std::fabs(Pzimax) < pf)
754  {
755  G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*cosTheta;
756  G4double p2 = Pzimax*Pzimax;
757  G4double dspz = std::sqrt(QCOE2)*
758  (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc
759  *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf));
760  sia += std::max(dspz,-1.0*sia);
761  }
762 
763  G4double XKN = EOEC+ECOE-1.0+cosTheta*cosTheta;
764 
765  //Differential cross section (per electron, in units of pi*classic_electr_radius**2)
766  G4double diffCS = ECOE*ECOE*XKN*sia;
767 
768  return diffCS;
769 }
energy
double energy
Definition: plottest35.C:25
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetVerbosityLevel()

G4int G4PenelopeComptonModel::GetVerbosityLevel ( )
inline

Definition at line 100 of file G4PenelopeComptonModel.hh.

100 {return verboseLevel;};

◆ Initialise()

void G4PenelopeComptonModel::Initialise ( const G4ParticleDefinition *  part,
const G4DataVector &   
)
virtual

Implements G4VEmModel.

Definition at line 100 of file G4PenelopeComptonModel.cc.

102 {
103  if (verboseLevel > 3)
104  G4cout << "Calling G4PenelopeComptonModel::Initialise()" << G4endl;
105 
106  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
107  //Issue warning if the AtomicDeexcitation has not been declared
108  if (!fAtomDeexcitation)
109  {
110  G4cout << G4endl;
111  G4cout << "WARNING from G4PenelopeComptonModel " << G4endl;
112  G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
113  G4cout << "any fluorescence/Auger emission." << G4endl;
114  G4cout << "Please make sure this is intended" << G4endl;
115  }
116 
117  SetParticle(part);
118 
119  if (IsMaster() && part == fParticle)
120  {
121 
122  if (verboseLevel > 0)
123  {
124  G4cout << "Penelope Compton model v2008 is initialized " << G4endl
125  << "Energy range: "
126  << LowEnergyLimit() / keV << " keV - "
127  << HighEnergyLimit() / GeV << " GeV";
128  }
129  //Issue a warning, if the model is going to be used down to a
130  //energy which is outside the validity of the model itself
131  if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
132  {
133  G4ExceptionDescription ed;
134  ed << "Using the Penelope Compton model outside its intrinsic validity range. "
135  << G4endl;
136  ed << "-> LowEnergyLimit() in process = " << LowEnergyLimit()/keV << "keV " << G4endl;
137  ed << "-> Instrinsic low-energy limit = " << fIntrinsicLowEnergyLimit/keV << "keV "
138  << G4endl;
139  ed << "Result of the simulation have to be taken with care" << G4endl;
140  G4Exception("G4PenelopeComptonModel::Initialise()",
141  "em2100",JustWarning,ed);
142  }
143  }
144 
145  if(isInitialised) return;
146  fParticleChange = GetParticleChangeForGamma();
147  isInitialised = true;
148 
149 }
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:112
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4VEmModel::IsMaster
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4cout
G4GLOB_DLL std::ostream G4cout
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
GeV
static const double GeV
Definition: G4SIunits.hh:214
G4PenelopeComptonModel::SetParticle
void SetParticle(const G4ParticleDefinition *)
Definition: G4PenelopeComptonModel.cc:948
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PenelopeComptonModel::fParticleChange
G4ParticleChangeForGamma * fParticleChange
Definition: G4PenelopeComptonModel.hh:100
G4PenelopeComptonModel::fParticle
const G4ParticleDefinition * fParticle
Definition: G4PenelopeComptonModel.hh:105
G4endl
#define G4endl
Definition: G4ios.hh:61
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:528
keV
static const double keV
Definition: G4SIunits.hh:213
G4PenelopeComptonModel::fAtomDeexcitation
G4VAtomDeexcitation * fAtomDeexcitation
Definition: G4PenelopeComptonModel.hh:131
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
Here is the call graph for this function:

◆ InitialiseLocal()

void G4PenelopeComptonModel::InitialiseLocal ( const G4ParticleDefinition *  part,
G4VEmModel *  masterModel 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 153 of file G4PenelopeComptonModel.cc.

155 {
156  if (verboseLevel > 3)
157  G4cout << "Calling G4PenelopeComptonModel::InitialiseLocal()" << G4endl;
158 
159  //
160  //Check that particle matches: one might have multiple master models (e.g.
161  //for e+ and e-).
162  //
163  if (part == fParticle)
164  {
165  //Get the const table pointers from the master to the workers
166  const G4PenelopeComptonModel* theModel =
167  static_cast<G4PenelopeComptonModel*> (masterModel);
168 
169  //Same verbosity for all workers, as the master
170  verboseLevel = theModel->verboseLevel;
171  }
172 
173  return;
174 }
G4cout
G4GLOB_DLL std::ostream G4cout
G4PenelopeComptonModel::fParticle
const G4ParticleDefinition * fParticle
Definition: G4PenelopeComptonModel.hh:105
G4endl
#define G4endl
Definition: G4ios.hh:61

◆ KleinNishinaCrossSection()

G4double G4PenelopeComptonModel::KleinNishinaCrossSection ( G4double  energy,
const G4Material *  material 
)
private

Definition at line 908 of file G4PenelopeComptonModel.cc.

910 {
911  // use Klein-Nishina formula
912  // total cross section in units of pi*classic_electr_radius^2
913 
914  G4double cs = 0;
915 
916  G4double ek =energy/electron_mass_c2;
917  G4double eks = ek*ek;
918  G4double ek2 = 1.0+ek+ek;
919  G4double ek1 = eks-ek2-1.0;
920 
921  G4double t0 = 1.0/ek2;
922  G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*std::log(t0)-(1.0/t0);
923 
924  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material);
925 
926  for (size_t i=0;i<theTable->size();i++)
927  {
928  G4PenelopeOscillator* theOsc = (*theTable)[i];
929  G4double ionEnergy = theOsc->GetIonisationEnergy();
930  G4double tau=(energy-ionEnergy)/energy;
931  if (tau > t0)
932  {
933  G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1*std::log(tau)-(1.0/tau);
934  G4double stre = theOsc->GetOscillatorStrength();
935 
936  cs += stre*(csu-csl);
937  }
938  }
939 
940  cs /= (ek*eks);
941 
942  return cs;
943 
944 }
G4PenelopeOscillatorManager::GetOscillatorTableCompton
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:364
energy
double energy
Definition: plottest35.C:25
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4PenelopeOscillatorTable
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
Definition: G4PenelopeOscillatorManager.hh:63
G4PenelopeComptonModel::oscManager
G4PenelopeOscillatorManager * oscManager
Definition: G4PenelopeComptonModel.hh:134
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ operator=()

G4PenelopeComptonModel& G4PenelopeComptonModel::operator= ( const G4PenelopeComptonModel &  right)
private

◆ OscillatorTotalCrossSection()

G4double G4PenelopeComptonModel::OscillatorTotalCrossSection ( G4double  energy,
G4PenelopeOscillator *  osc 
)
private

Definition at line 773 of file G4PenelopeComptonModel.cc.

774 {
775  //Total cross section (integrated) for the given oscillator in units of
776  //pi*classic_electr_radius^2
777 
778  //Integrate differential cross section for each oscillator
779  G4double stre = osc->GetOscillatorStrength();
780 
781  // here one uses the using the 20-point
782  // Gauss quadrature method with an adaptive bipartition scheme
783  const G4int npoints=10;
784  const G4int ncallsmax=20000;
785  const G4int nst=256;
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};
794 
795  G4double MaxError = 1e-5;
796  //Error control
797  G4double Ctol = std::min(std::max(MaxError,1e-13),1e-02);
798  G4double Ptol = 0.01*Ctol;
799  G4double Err=1e35;
800 
801  //Gauss integration from -1 to 1
802  G4double LowPoint = -1.0;
803  G4double HighPoint = 1.0;
804 
805  G4double h=HighPoint-LowPoint;
806  G4double sumga=0.0;
807  G4double a=0.5*(HighPoint-LowPoint);
808  G4double b=0.5*(HighPoint+LowPoint);
809  G4double c=a*Abscissas[0];
810  G4double d= Weights[0]*
811  (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
812  for (G4int i=2;i<=npoints;i++)
813  {
814  c=a*Abscissas[i-1];
815  d += Weights[i-1]*
816  (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
817  }
818  G4int icall = 2*npoints;
819  G4int LH=1;
820  G4double S[nst],x[nst],sn[nst],xrn[nst];
821  S[0]=d*a;
822  x[0]=LowPoint;
823 
824  G4bool loopAgain = true;
825 
826  //Adaptive bipartition scheme
827  do{
828  G4double h0=h;
829  h=0.5*h; //bipartition
830  G4double sumr=0;
831  G4int LHN=0;
832  G4double si,xa,xb,xc;
833  for (G4int i=1;i<=LH;i++){
834  si=S[i-1];
835  xa=x[i-1];
836  xb=xa+h;
837  xc=xa+h0;
838  a=0.5*(xb-xa);
839  b=0.5*(xb+xa);
840  c=a*Abscissas[0];
841  G4double dLocal = Weights[0]*
842  (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
843 
844  for (G4int j=1;j<npoints;j++)
845  {
846  c=a*Abscissas[j];
847  dLocal += Weights[j]*
848  (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
849  }
850  G4double s1=dLocal*a;
851  a=0.5*(xc-xb);
852  b=0.5*(xc+xb);
853  c=a*Abscissas[0];
854  dLocal=Weights[0]*
855  (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
856 
857  for (G4int j=1;j<npoints;j++)
858  {
859  c=a*Abscissas[j];
860  dLocal += Weights[j]*
861  (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
862  }
863  G4double s2=dLocal*a;
864  icall=icall+4*npoints;
865  G4double s12=s1+s2;
866  if (std::abs(s12-si)<std::max(Ptol*std::abs(s12),1e-35))
867  sumga += s12;
868  else
869  {
870  sumr += s12;
871  LHN += 2;
872  sn[LHN-1]=s2;
873  xrn[LHN-1]=xb;
874  sn[LHN-2]=s1;
875  xrn[LHN-2]=xa;
876  }
877 
878  if (icall>ncallsmax || LHN>nst)
879  {
880  G4cout << "G4PenelopeComptonModel: " << G4endl;
881  G4cout << "LowPoint: " << LowPoint << ", High Point: " << HighPoint << G4endl;
882  G4cout << "Tolerance: " << MaxError << 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;
886  loopAgain = false;
887  }
888  }
889  Err=std::abs(sumr)/std::max(std::abs(sumr+sumga),1e-35);
890  if (Err < Ctol || LHN == 0)
891  loopAgain = false; //end of cycle
892  LH=LHN;
893  for (G4int i=0;i<LH;i++)
894  {
895  S[i]=sn[i];
896  x[i]=xrn[i];
897  }
898  }while(Ctol < 1.0 && loopAgain);
899 
900 
901  G4double xs = stre*sumga;
902 
903  return xs;
904 }
S
double S(double temp)
Definition: G4DNAElectronHoleRecombination.cc:76
d
Float_t d
Definition: advanced/microbeam/plot.C:237
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
energy
double energy
Definition: plottest35.C:25
G4bool
bool G4bool
Definition: G4Types.hh:79
G4PenelopeComptonModel::DifferentialCrossSection
G4double DifferentialCrossSection(G4double cdt, G4double energy, G4PenelopeOscillator *osc)
Definition: G4PenelopeComptonModel.cc:711
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleSecondaries()

void G4PenelopeComptonModel::SampleSecondaries ( std::vector< G4DynamicParticle *> *  fvect,
const G4MaterialCutsCouple *  couple,
const G4DynamicParticle *  aDynamicGamma,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 269 of file G4PenelopeComptonModel.cc.

274 {
275 
276  // Penelope model v2008 to sample the Compton scattering final state.
277  // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
278  // The model determines also the original shell from which the electron is expelled,
279  // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
280  //
281  // The final state for Compton scattering is calculated according to the Klein-Nishina
282  // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and
283  // one can assume that the target electron is at rest.
284  // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega,
285  // to sample the scattering angle and the energy of the emerging electron, which is
286  // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is
287  // used to sample cos(theta). The efficiency increases monotonically with photon energy and is
288  // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV,
289  // respectively.
290  // The parametrization includes the J(p) distribution profiles for the atomic shells, that are
291  // tabulated
292  // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201.
293  // Doppler broadening is included.
294  //
295 
296  if (verboseLevel > 3)
297  G4cout << "Calling SampleSecondaries() of G4PenelopeComptonModel" << G4endl;
298 
299  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
300 
301  // do nothing below the threshold
302  // should never get here because the XS is zero below the limit
303  if(photonEnergy0 < LowEnergyLimit())
304  return;
305 
306  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
307  const G4Material* material = couple->GetMaterial();
308 
309  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material);
310 
311  const G4int nmax = 64;
312  G4double rn[nmax]={0.0};
313  G4double pac[nmax]={0.0};
314 
315  G4double S=0.0;
316  G4double epsilon = 0.0;
317  G4double cosTheta = 1.0;
318  G4double hartreeFunc = 0.0;
319  G4double oscStren = 0.0;
320  size_t numberOfOscillators = theTable->size();
321  size_t targetOscillator = 0;
322  G4double ionEnergy = 0.0*eV;
323 
324  G4double ek = photonEnergy0/electron_mass_c2;
325  G4double ek2 = 2.*ek+1.0;
326  G4double eks = ek*ek;
327  G4double ek1 = eks-ek2-1.0;
328 
329  G4double taumin = 1.0/ek2;
330  G4double a1 = std::log(ek2);
331  G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2);
332 
333  G4double TST = 0;
334  G4double tau = 0.;
335 
336  //If the incoming photon is above 5 MeV, the quicker approach based on the
337  //pure Klein-Nishina formula is used
338  if (photonEnergy0 > 5*MeV)
339  {
340  do{
341  do{
342  if ((a2*G4UniformRand()) < a1)
343  tau = std::pow(taumin,G4UniformRand());
344  else
345  tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
346  //rejection function
347  TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
348  }while (G4UniformRand()> TST);
349  epsilon=tau;
350  cosTheta = 1.0 - (1.0-tau)/(ek*tau);
351 
352  //Target shell electrons
353  TST = oscManager->GetTotalZ(material)*G4UniformRand();
354  targetOscillator = numberOfOscillators-1; //last level
355  S=0.0;
356  G4bool levelFound = false;
357  for (size_t j=0;j<numberOfOscillators && !levelFound; j++)
358  {
359  S += (*theTable)[j]->GetOscillatorStrength();
360  if (S > TST)
361  {
362  targetOscillator = j;
363  levelFound = true;
364  }
365  }
366  //check whether the level is valid
367  ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
368  }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
369  }
370  else //photonEnergy0 < 5 MeV
371  {
372  //Incoherent scattering function for theta=PI
373  G4double s0=0.0;
374  G4double pzomc=0.0;
375  G4double rni=0.0;
376  G4double aux=0.0;
377  for (size_t i=0;i<numberOfOscillators;i++)
378  {
379  ionEnergy = (*theTable)[i]->GetIonisationEnergy();
380  if (photonEnergy0 > ionEnergy)
381  {
382  G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
383  hartreeFunc = (*theTable)[i]->GetHartreeFactor();
384  oscStren = (*theTable)[i]->GetOscillatorStrength();
385  pzomc = hartreeFunc*(aux2-electron_mass_c2*ionEnergy)/
386  (electron_mass_c2*std::sqrt(2.0*aux2+ionEnergy*ionEnergy));
387  if (pzomc > 0)
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));
390  else
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));
393  s0 += oscStren*rni;
394  }
395  }
396  //Sampling tau
397  G4double cdt1 = 0.;
398  do
399  {
400  if ((G4UniformRand()*a2) < a1)
401  tau = std::pow(taumin,G4UniformRand());
402  else
403  tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
404  cdt1 = (1.0-tau)/(ek*tau);
405  //Incoherent scattering function
406  S = 0.;
407  for (size_t i=0;i<numberOfOscillators;i++)
408  {
409  ionEnergy = (*theTable)[i]->GetIonisationEnergy();
410  if (photonEnergy0 > ionEnergy) //sum only on excitable levels
411  {
412  aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
413  hartreeFunc = (*theTable)[i]->GetHartreeFactor();
414  oscStren = (*theTable)[i]->GetOscillatorStrength();
415  pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/
416  (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
417  if (pzomc > 0)
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));
420  else
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));
423  S += oscStren*rn[i];
424  pac[i] = S;
425  }
426  else
427  pac[i] = S-1e-6;
428  }
429  //Rejection function
430  TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
431  }while ((G4UniformRand()*s0) > TST);
432 
433  cosTheta = 1.0 - cdt1;
434  G4double fpzmax=0.0,fpz=0.0;
435  G4double A=0.0;
436  //Target electron shell
437  do
438  {
439  do
440  {
441  TST = S*G4UniformRand();
442  targetOscillator = numberOfOscillators-1; //last level
443  G4bool levelFound = false;
444  for (size_t i=0;i<numberOfOscillators && !levelFound;i++)
445  {
446  if (pac[i]>TST)
447  {
448  targetOscillator = i;
449  levelFound = true;
450  }
451  }
452  A = G4UniformRand()*rn[targetOscillator];
453  hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor();
454  oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength();
455  if (A < 0.5)
456  pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
457  (std::sqrt(2.0)*hartreeFunc);
458  else
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);
462 
463  // F(EP) rejection
464  G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
465  G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
466  if (AF > 0)
467  fpzmax = 1.0+AF*0.2;
468  else
469  fpzmax = 1.0-AF*0.2;
470  fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
471  }while ((fpzmax*G4UniformRand())>fpz);
472 
473  //Energy of the scattered photon
474  G4double T = pzomc*pzomc;
475  G4double b1 = 1.0-T*tau*tau;
476  G4double b2 = 1.0-T*tau*cosTheta;
477  if (pzomc > 0.0)
478  epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
479  else
480  epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
481  } //energy < 5 MeV
482 
483  //Ok, the kinematics has been calculated.
484  G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
485  G4double phi = twopi * G4UniformRand() ;
486  G4double dirx = sinTheta * std::cos(phi);
487  G4double diry = sinTheta * std::sin(phi);
488  G4double dirz = cosTheta ;
489 
490  // Update G4VParticleChange for the scattered photon
491  G4ThreeVector photonDirection1(dirx,diry,dirz);
492  photonDirection1.rotateUz(photonDirection0);
493  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
494 
495  G4double photonEnergy1 = epsilon * photonEnergy0;
496 
497  if (photonEnergy1 > 0.)
498  fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
499  else
500  {
501  fParticleChange->SetProposedKineticEnergy(0.) ;
502  fParticleChange->ProposeTrackStatus(fStopAndKill);
503  }
504 
505  //Create scattered electron
506  G4double diffEnergy = photonEnergy0*(1-epsilon);
507  ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
508 
509  G4double Q2 =
510  photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
511  G4double cosThetaE = 0.; //scattering angle for the electron
512 
513  if (Q2 > 1.0e-12)
514  cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
515  else
516  cosThetaE = 1.0;
517  G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
518 
519  //Now, try to handle fluorescence
520  //Notice: merged levels are indicated with Z=0 and flag=30
521  G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
522  G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ();
523 
524  //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
525  G4double bindingEnergy = 0.*eV;
526  const G4AtomicShell* shell = 0;
527 
528  //Real level
529  if (Z > 0 && shFlag<30)
530  {
531  shell = fTransitionManager->Shell(Z,shFlag-1);
532  bindingEnergy = shell->BindingEnergy();
533  }
534 
535  G4double ionEnergyInPenelopeDatabase = ionEnergy;
536  //protection against energy non-conservation
537  ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase);
538 
539  //subtract the excitation energy. If not emitted by fluorescence
540  //the ionization energy is deposited as local energy deposition
541  G4double eKineticEnergy = diffEnergy - ionEnergy;
542  G4double localEnergyDeposit = ionEnergy;
543  G4double energyInFluorescence = 0.; //testing purposes only
544  G4double energyInAuger = 0; //testing purposes
545 
546  if (eKineticEnergy < 0)
547  {
548  //It means that there was some problem/mismatch between the two databases.
549  //Try to make it work
550  //In this case available Energy (diffEnergy) < ionEnergy
551  //Full residual energy is deposited locally
552  localEnergyDeposit = diffEnergy;
553  eKineticEnergy = 0.0;
554  }
555 
556  //the local energy deposit is what remains: part of this may be spent for fluorescence.
557  //Notice: shell might be NULL (invalid!) if shFlag=30. Must be protected
558  //Now, take care of fluorescence, if required
559  if (fAtomDeexcitation && shell)
560  {
561  G4int index = couple->GetIndex();
562  if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
563  {
564  size_t nBefore = fvect->size();
565  fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
566  size_t nAfter = fvect->size();
567 
568  if (nAfter > nBefore) //actual production of fluorescence
569  {
570  for (size_t j=nBefore;j<nAfter;j++) //loop on products
571  {
572  G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
573  localEnergyDeposit -= itsEnergy;
574  if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
575  energyInFluorescence += itsEnergy;
576  else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
577  energyInAuger += itsEnergy;
578 
579  }
580  }
581 
582  }
583  }
584 
585 
586  /*
587  if(DeexcitationFlag() && Z > 5 && shellId>0) {
588 
589  const G4ProductionCutsTable* theCoupleTable=
590  G4ProductionCutsTable::GetProductionCutsTable();
591 
592  size_t index = couple->GetIndex();
593  G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
594  G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
595 
596  // Generation of fluorescence
597  // Data in EADL are available only for Z > 5
598  // Protection to avoid generating photons in the unphysical case of
599  // shell binding energy > photon energy
600  if (localEnergyDeposit > cutg || localEnergyDeposit > cute)
601  {
602  G4DynamicParticle* aPhoton;
603  deexcitationManager.SetCutForSecondaryPhotons(cutg);
604  deexcitationManager.SetCutForAugerElectrons(cute);
605 
606  photonVector = deexcitationManager.GenerateParticles(Z,shellId);
607  if(photonVector)
608  {
609  size_t nPhotons = photonVector->size();
610  for (size_t k=0; k<nPhotons; k++)
611  {
612  aPhoton = (*photonVector)[k];
613  if (aPhoton)
614  {
615  G4double itsEnergy = aPhoton->GetKineticEnergy();
616  G4bool keepIt = false;
617  if (itsEnergy <= localEnergyDeposit)
618  {
619  //check if good!
620  if(aPhoton->GetDefinition() == G4Gamma::Gamma()
621  && itsEnergy >= cutg)
622  {
623  keepIt = true;
624  energyInFluorescence += itsEnergy;
625  }
626  if (aPhoton->GetDefinition() == G4Electron::Electron() &&
627  itsEnergy >= cute)
628  {
629  energyInAuger += itsEnergy;
630  keepIt = true;
631  }
632  }
633  //good secondary, register it
634  if (keepIt)
635  {
636  localEnergyDeposit -= itsEnergy;
637  fvect->push_back(aPhoton);
638  }
639  else
640  {
641  delete aPhoton;
642  (*photonVector)[k] = 0;
643  }
644  }
645  }
646  delete photonVector;
647  }
648  }
649  }
650  */
651 
652 
653  //Always produce explicitely the electron
654  G4DynamicParticle* electron = 0;
655 
656  G4double xEl = sinThetaE * std::cos(phi+pi);
657  G4double yEl = sinThetaE * std::sin(phi+pi);
658  G4double zEl = cosThetaE;
659  G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
660  eDirection.rotateUz(photonDirection0);
661  electron = new G4DynamicParticle (G4Electron::Electron(),
662  eDirection,eKineticEnergy) ;
663  fvect->push_back(electron);
664 
665 
666  if (localEnergyDeposit < 0)
667  {
668  G4cout << "WARNING-"
669  << "G4PenelopeComptonModel::SampleSecondaries - Negative energy deposit"
670  << G4endl;
671  localEnergyDeposit=0.;
672  }
673  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
674 
675  G4double electronEnergy = 0.;
676  if (electron)
677  electronEnergy = eKineticEnergy;
678  if (verboseLevel > 1)
679  {
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;
688  if (energyInAuger)
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 <<
693  " keV" << G4endl;
694  G4cout << "-----------------------------------------------------------" << G4endl;
695  }
696  if (verboseLevel > 0)
697  {
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;
706  }
707 }
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
MeV
static const double MeV
Definition: G4SIunits.hh:211
G4VAtomDeexcitation::CheckDeexcitationActiveRegion
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:261
S
double S(double temp)
Definition: G4DNAElectronHoleRecombination.cc:76
G4AtomicTransitionManager::Shell
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
Definition: G4AtomicTransitionManager.cc:95
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
index
Int_t index
Definition: brachytherapy/macro.C:7
a1
static const G4double a1
Definition: G4BetheHeitlerModel.cc:72
G4Electron::Definition
static G4Electron * Definition()
Definition: G4Electron.cc:49
G4int
int G4int
Definition: G4Types.hh:78
G4PenelopeOscillatorManager::GetOscillatorTableCompton
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:364
eplot.material
string material
Definition: eplot.py:19
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:59
Z
Float_t Z
Definition: advanced/microbeam/plot.C:39
G4PenelopeComptonModel::fTransitionManager
const G4AtomicTransitionManager * fTransitionManager
Definition: G4PenelopeComptonModel.hh:132
G4bool
bool G4bool
Definition: G4Types.hh:79
twopi
static const double twopi
Definition: G4SIunits.hh:75
b2
static const G4double b2
Definition: G4BetheHeitlerModel.cc:76
nmax
const G4int nmax
Definition: G4UrQMD1_3Interface.hh:69
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
eV
static const double eV
Definition: G4SIunits.hh:212
G4PenelopeComptonModel::fParticleChange
G4ParticleChangeForGamma * fParticleChange
Definition: G4PenelopeComptonModel.hh:100
pi
static const double pi
Definition: G4SIunits.hh:74
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle *> *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:274
G4PenelopeOscillatorTable
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
Definition: G4PenelopeOscillatorManager.hh:63
b1
static const G4double b1
Definition: G4BetheHeitlerModel.cc:76
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4PenelopeComptonModel::oscManager
G4PenelopeOscillatorManager * oscManager
Definition: G4PenelopeComptonModel.hh:134
keV
static const double keV
Definition: G4SIunits.hh:213
G4double
double G4double
Definition: G4Types.hh:76
G4AtomicShell::BindingEnergy
G4double BindingEnergy() const
Definition: G4AtomicShell.hh:79
G4PenelopeOscillatorManager::GetTotalZ
G4double GetTotalZ(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:267
epsilon
double epsilon(double density, double temperature)
Definition: G4DNAElectronHoleRecombination.cc:95
G4PenelopeComptonModel::fAtomDeexcitation
G4VAtomDeexcitation * fAtomDeexcitation
Definition: G4PenelopeComptonModel.hh:131
a2
static const G4double a2
Definition: G4BetheHeitlerModel.cc:72
G4Gamma::Definition
static G4Gamma * Definition()
Definition: G4Gamma.cc:49
Here is the call graph for this function:

◆ SetParticle()

void G4PenelopeComptonModel::SetParticle ( const G4ParticleDefinition *  p)
private

Definition at line 948 of file G4PenelopeComptonModel.cc.

949 {
950  if(!fParticle) {
951  fParticle = p;
952  }
953 }
G4PenelopeComptonModel::fParticle
const G4ParticleDefinition * fParticle
Definition: G4PenelopeComptonModel.hh:105
Here is the caller graph for this function:

◆ SetVerbosityLevel()

void G4PenelopeComptonModel::SetVerbosityLevel ( G4int  lev)
inline

Definition at line 99 of file G4PenelopeComptonModel.hh.

99 {verboseLevel = lev;};

Member Data Documentation

◆ fAtomDeexcitation

G4VAtomDeexcitation* G4PenelopeComptonModel::fAtomDeexcitation
private

Definition at line 131 of file G4PenelopeComptonModel.hh.

◆ fIntrinsicHighEnergyLimit

G4double G4PenelopeComptonModel::fIntrinsicHighEnergyLimit
private

Definition at line 125 of file G4PenelopeComptonModel.hh.

◆ fIntrinsicLowEnergyLimit

G4double G4PenelopeComptonModel::fIntrinsicLowEnergyLimit
private

Definition at line 124 of file G4PenelopeComptonModel.hh.

◆ fParticle

const G4ParticleDefinition* G4PenelopeComptonModel::fParticle
protected

Definition at line 105 of file G4PenelopeComptonModel.hh.

◆ fParticleChange

G4ParticleChangeForGamma* G4PenelopeComptonModel::fParticleChange
protected

Definition at line 100 of file G4PenelopeComptonModel.hh.

◆ fTransitionManager

const G4AtomicTransitionManager* G4PenelopeComptonModel::fTransitionManager
private

Definition at line 132 of file G4PenelopeComptonModel.hh.

◆ isInitialised

G4bool G4PenelopeComptonModel::isInitialised
private

Definition at line 129 of file G4PenelopeComptonModel.hh.

◆ oscManager

G4PenelopeOscillatorManager* G4PenelopeComptonModel::oscManager
private

Definition at line 134 of file G4PenelopeComptonModel.hh.

◆ verboseLevel

G4int G4PenelopeComptonModel::verboseLevel
private

Definition at line 127 of file G4PenelopeComptonModel.hh.

The documentation for this class was generated from the following files: