Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 level, const G4ParticleDefinition *, G4double kineticEnergy)
 
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 * > *)
 
virtual 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 G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
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 SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Attributes

G4ParticleChangeForGammafParticleChange
 
const G4ParticleDefinitionfParticle
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () 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::G4PenelopeComptonModel ( const G4ParticleDefinition p = 0,
const G4String processName = "PenCompton" 
)

Definition at line 64 of file G4PenelopeComptonModel.cc.

67  isInitialised(false),fAtomDeexcitation(0),
68  oscManager(0)
69 {
70  fIntrinsicLowEnergyLimit = 100.0*eV;
71  fIntrinsicHighEnergyLimit = 100.0*GeV;
72  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
73  //
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(const G4String &nam)
Definition: G4VEmModel.cc:68
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
static constexpr double eV
Definition: G4SIunits.hh:215
static G4PenelopeOscillatorManager * GetOscillatorManager()
G4ParticleChangeForGamma * fParticleChange
const G4ParticleDefinition * fParticle
static constexpr double GeV
Definition: G4SIunits.hh:217
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:788
static G4AtomicTransitionManager * Instance()

Here is the call graph for this function:

G4PenelopeComptonModel::~G4PenelopeComptonModel ( )
virtual

Definition at line 95 of file G4PenelopeComptonModel.cc.

96 {;}

Member Function Documentation

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 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
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
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 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
static constexpr double mm
Definition: G4SIunits.hh:115
const G4String & GetName() const
Definition: G4Material.hh:178
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:411
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
G4GLOB_DLL std::ostream G4cout
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
G4double energy(const ThreeVector &p, const G4double m)
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4int G4PenelopeComptonModel::GetVerbosityLevel ( )
inline

Definition at line 100 of file G4PenelopeComptonModel.hh.

100 {return verboseLevel;};
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  {
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;
147  isInitialised = true;
148 
149 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
static G4LossTableManager * Instance()
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
G4GLOB_DLL std::ostream G4cout
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ParticleChangeForGamma * fParticleChange
const G4ParticleDefinition * fParticle
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
G4VAtomDeexcitation * AtomDeexcitation()
static constexpr double keV
Definition: G4SIunits.hh:216
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132

Here is the call graph for this function:

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 }
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * fParticle
#define G4endl
Definition: G4ios.hh:61
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  {
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-
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 nullptr (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
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 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
double S(double temp)
G4double GetKineticEnergy() const
static G4Electron * Definition()
Definition: G4Electron.cc:49
G4double BindingEnergy() const
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
string material
Definition: eplot.py:19
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4double ek
static constexpr double eV
Definition: G4SIunits.hh:215
const G4int nmax
float electron_mass_c2
Definition: hepunit.py:274
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4ParticleChangeForGamma * fParticleChange
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
G4double bindingEnergy(G4int A, G4int Z)
static constexpr double keV
Definition: G4SIunits.hh:216
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
G4double GetTotalZ(const G4Material *)
double epsilon(double density, double temperature)
static G4Gamma * Definition()
Definition: G4Gamma.cc:49
const G4Material * GetMaterial() const

Here is the call graph for this function:

void G4PenelopeComptonModel::SetVerbosityLevel ( G4int  lev)
inline

Definition at line 99 of file G4PenelopeComptonModel.hh.

99 {verboseLevel = lev;};

Member Data Documentation

const G4ParticleDefinition* G4PenelopeComptonModel::fParticle
protected

Definition at line 105 of file G4PenelopeComptonModel.hh.

G4ParticleChangeForGamma* G4PenelopeComptonModel::fParticleChange
protected

Definition at line 100 of file G4PenelopeComptonModel.hh.


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