Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4LowEPPolarizedComptonModel Class Reference

#include <G4LowEPPolarizedComptonModel.hh>

Inheritance diagram for G4LowEPPolarizedComptonModel:
Inheritance graph
[legend]
Collaboration diagram for G4LowEPPolarizedComptonModel:
Collaboration graph
[legend]

Public Member Functions

 G4LowEPPolarizedComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="LowEPComptonModel")
 
virtual ~G4LowEPPolarizedComptonModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=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 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=nullptr)
 
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)
 

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 *)
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 88 of file G4LowEPPolarizedComptonModel.hh.

Constructor & Destructor Documentation

G4LowEPPolarizedComptonModel::G4LowEPPolarizedComptonModel ( const G4ParticleDefinition *  p = 0,
const G4String &  nam = "LowEPComptonModel" 
)

Definition at line 80 of file G4LowEPPolarizedComptonModel.cc.

82  : G4VEmModel(nam),isInitialised(false)
83 {
84  verboseLevel=1 ;
85  // Verbosity scale:
86  // 0 = nothing
87  // 1 = warning for energy non-conservation
88  // 2 = details of energy budget
89  // 3 = calculation of cross sections, file openings, sampling of atoms
90  // 4 = entering in methods
91 
92  if( verboseLevel>1 ) {
93  G4cout << "Low energy photon Compton model is constructed " << G4endl;
94  }
95 
96  //Mark this model as "applicable" for atomic deexcitation
97  SetDeexcitationFlag(true);
98 
99  fParticleChange = 0;
100  fAtomDeexcitation = 0;
101 }
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:780

Here is the call graph for this function:

G4LowEPPolarizedComptonModel::~G4LowEPPolarizedComptonModel ( )
virtual

Definition at line 105 of file G4LowEPPolarizedComptonModel.cc.

106 {
107  if(IsMaster()) {
108  delete shellData;
109  shellData = 0;
110  delete profileData;
111  profileData = 0;
112  }
113 }
G4VEmModel::IsMaster
G4bool IsMaster() const
Definition: G4VEmModel.hh:717

Here is the call graph for this function:

Member Function Documentation

G4double G4LowEPPolarizedComptonModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  ,
G4double  kinEnergy,
G4double  Z,
G4double  A = 0,
G4double  cut = 0,
G4double  emax = DBL_MAX 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 246 of file G4LowEPPolarizedComptonModel.cc.

250 {
251  if (verboseLevel > 3) {
252  G4cout << "G4LowEPPolarizedComptonModel::ComputeCrossSectionPerAtom()"
253  << G4endl;
254  }
255  G4double cs = 0.0;
256 
257  if (GammaEnergy < LowEnergyLimit()) { return 0.0; }
258 
259  G4int intZ = G4lrint(Z);
260  if(intZ < 1 || intZ > maxZ) { return cs; }
261 
262  G4LPhysicsFreeVector* pv = data[intZ];
263 
264  // if element was not initialised
265  // do initialisation safely for MT mode
266  if(!pv)
267  {
268  InitialiseForElement(0, intZ);
269  pv = data[intZ];
270  if(!pv) { return cs; }
271  }
272 
273  G4int n = pv->GetVectorLength() - 1;
274  G4double e1 = pv->Energy(0);
275  G4double e2 = pv->Energy(n);
276 
277  if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*e1)*pv->Value(e1); }
278  else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; }
279  else if(GammaEnergy > e2) { cs = pv->Value(e2)/GammaEnergy; }
280 
281  return cs;
282 }
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
G4LowEPPolarizedComptonModel::InitialiseForElement
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
Definition: G4LowEPPolarizedComptonModel.cc:709
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4int
int G4int
Definition: G4Types.hh:78
data
const XML_Char const XML_Char * data
Definition: expat.h:268
G4cout
G4GLOB_DLL std::ostream G4cout
G4PhysicsVector::Energy
G4double Energy(size_t index) const
G4PhysicsVector::Value
G4double Value(G4double theEnergy, size_t &lastidx) const
Definition: G4PhysicsVector.cc:498
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

void G4LowEPPolarizedComptonModel::Initialise ( const G4ParticleDefinition *  particle,
const G4DataVector &  cuts 
)
virtual

Implements G4VEmModel.

Definition at line 117 of file G4LowEPPolarizedComptonModel.cc.

119 {
120  if (verboseLevel > 1) {
121  G4cout << "Calling G4LowEPPolarizedComptonModel::Initialise()" << G4endl;
122  }
123 
124  // Initialise element selector
125 
126  if(IsMaster()) {
127 
128  // Access to elements
129 
130  char* path = getenv("G4LEDATA");
131 
132  G4ProductionCutsTable* theCoupleTable =
133  G4ProductionCutsTable::GetProductionCutsTable();
134  G4int numOfCouples = theCoupleTable->GetTableSize();
135 
136  for(G4int i=0; i<numOfCouples; ++i) {
137  const G4Material* material =
138  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
139  const G4ElementVector* theElementVector = material->GetElementVector();
140  G4int nelm = material->GetNumberOfElements();
141 
142  for (G4int j=0; j<nelm; ++j) {
143  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
144  if(Z < 1) { Z = 1; }
145  else if(Z > maxZ){ Z = maxZ; }
146 
147  if( (!data[Z]) ) { ReadData(Z, path); }
148  }
149  }
150 
151  // For Doppler broadening
152  if(!shellData) {
153  shellData = new G4ShellData();
154  shellData->SetOccupancyData();
155  G4String file = "/doppler/shell-doppler";
156  shellData->LoadData(file);
157  }
158  if(!profileData) { profileData = new G4DopplerProfile(); }
159 
160  InitialiseElementSelectors(particle, cuts);
161  }
162 
163  if (verboseLevel > 2) {
164  G4cout << "Loaded cross section files" << G4endl;
165  }
166 
167  if( verboseLevel>1 ) {
168  G4cout << "G4LowEPPolarizedComptonModel is initialized " << G4endl
169  << "Energy range: "
170  << LowEnergyLimit() / eV << " eV - "
171  << HighEnergyLimit() / GeV << " GeV"
172  << G4endl;
173  }
174 
175  if(isInitialised) { return; }
176 
177  fParticleChange = GetParticleChangeForGamma();
178  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
179  isInitialised = true;
180 }
G4ShellData::SetOccupancyData
void SetOccupancyData()
Definition: G4ShellData.hh:70
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:112
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
G4VEmModel::InitialiseElementSelectors
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
G4VEmModel::IsMaster
G4bool IsMaster() const
Definition: G4VEmModel.hh:717
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
G4int
int G4int
Definition: G4Types.hh:78
G4ShellData::LoadData
void LoadData(const G4String &fileName)
Definition: G4ShellData.cc:234
data
const XML_Char const XML_Char * data
Definition: expat.h:268
G4cout
G4GLOB_DLL std::ostream G4cout
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:444
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150

Here is the call graph for this function:

void G4LowEPPolarizedComptonModel::InitialiseForElement ( const G4ParticleDefinition *  ,
G4int  Z 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 709 of file G4LowEPPolarizedComptonModel.cc.

711 {
712  G4AutoLock l(&LowEPPolarizedComptonModelMutex);
713  if(!data[Z]) { ReadData(Z); }
714  l.unlock();
715 }
data
const XML_Char const XML_Char * data
Definition: expat.h:268
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

Here is the caller graph for this function:

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

Reimplemented from G4VEmModel.

Definition at line 184 of file G4LowEPPolarizedComptonModel.cc.

186 {
187  SetElementSelectors(masterModel->GetElementSelectors());
188 }
G4VEmModel::GetElementSelectors
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:801
G4VEmModel::SetElementSelectors
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:809

Here is the call graph for this function:

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

Implements G4VEmModel.

Definition at line 286 of file G4LowEPPolarizedComptonModel.cc.

290 {
291 
292  //Determine number of digits (in decimal base) that G4double can accurately represent
293  G4double g4d_order = G4double(numeric_limits<G4double>::digits10);
294  G4double g4d_limit = std::pow(10.,-g4d_order);
295 
296  // The scattered gamma energy is sampled according to Klein - Nishina formula.
297  // then accepted or rejected depending on the Scattering Function multiplied
298  // by factor from Klein - Nishina formula.
299  // Expression of the angular distribution as Klein Nishina
300  // angular and energy distribution and Scattering fuctions is taken from
301  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
302  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
303  // data are interpolated while in the article they are fitted.
304  // Reference to the article is from J. Stepanek New Photon, Positron
305  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
306  // TeV (draft).
307  // The random number techniques of Butcher & Messel are used
308  // (Nucl Phys 20(1960),15).
309 
310 
311  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV;
312 
313  if (verboseLevel > 3) {
314  G4cout << "G4LowEPPolarizedComptonModel::SampleSecondaries() E(MeV)= "
315  << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
316  << G4endl;
317  }
318  // do nothing below the threshold
319  // should never get here because the XS is zero below the limit
320  if (photonEnergy0 < LowEnergyLimit())
321  return ;
322 
323  G4double e0m = photonEnergy0 / electron_mass_c2 ;
324  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
325 
326 
327  // Polarisation: check orientation of photon propagation direction and polarisation
328  // Fix if needed
329 
330  G4ThreeVector photonPolarization0 = aDynamicGamma->GetPolarization();
331 
332  // Check if polarisation vector is perpendicular and fix if not
333 
334  if (!(photonPolarization0.isOrthogonal(photonDirection0, 1e-6))||(photonPolarization0.mag()==0))
335  {
336  photonPolarization0 = GetRandomPolarization(photonDirection0);
337  }
338 
339  else
340  {
341  if ((photonPolarization0.howOrthogonal(photonDirection0) !=0) && (photonPolarization0.howOrthogonal(photonDirection0) > g4d_limit))
342  {
343  photonPolarization0 = GetPerpendicularPolarization(photonDirection0,photonPolarization0);
344  }
345  }
346 
347  // Select randomly one element in the current material
348 
349  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
350  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
351  G4int Z = (G4int)elm->GetZ();
352 
353  G4double LowEPPCepsilon0 = 1. / (1. + 2. * e0m);
354  G4double LowEPPCepsilon0Sq = LowEPPCepsilon0 * LowEPPCepsilon0;
355  G4double alpha1 = -std::log(LowEPPCepsilon0);
356  G4double alpha2 = 0.5 * (1. - LowEPPCepsilon0Sq);
357 
358  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
359 
360  // Sample the energy of the scattered photon
361  G4double LowEPPCepsilon;
362  G4double LowEPPCepsilonSq;
363  G4double oneCosT;
364  G4double sinT2;
365  G4double gReject;
366 
367  if (verboseLevel > 3) {
368  G4cout << "Started loop to sample gamma energy" << G4endl;
369  }
370 
371  do
372  {
373  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
374  {
375  LowEPPCepsilon = G4Exp(-alpha1 * G4UniformRand());
376  LowEPPCepsilonSq = LowEPPCepsilon * LowEPPCepsilon;
377  }
378  else
379  {
380  LowEPPCepsilonSq = LowEPPCepsilon0Sq + (1. - LowEPPCepsilon0Sq) * G4UniformRand();
381  LowEPPCepsilon = std::sqrt(LowEPPCepsilonSq);
382  }
383 
384  oneCosT = (1. - LowEPPCepsilon) / ( LowEPPCepsilon * e0m);
385  sinT2 = oneCosT * (2. - oneCosT);
386  G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
387  G4double scatteringFunction = ComputeScatteringFunction(x, Z);
388  gReject = (1. - LowEPPCepsilon * sinT2 / (1. + LowEPPCepsilonSq)) * scatteringFunction;
389 
390  } while(gReject < G4UniformRand()*Z);
391 
392  G4double cosTheta = 1. - oneCosT;
393  G4double sinTheta = std::sqrt(sinT2);
394  G4double phi = SetPhi(LowEPPCepsilon,sinT2);
395  G4double dirx = sinTheta * std::cos(phi);
396  G4double diry = sinTheta * std::sin(phi);
397  G4double dirz = cosTheta ;
398 
399  // Set outgoing photon polarization
400 
401  G4ThreeVector photonPolarization1 = SetNewPolarization(LowEPPCepsilon,
402  sinT2,
403  phi,
404  cosTheta);
405 
406  // Scatter photon energy and Compton electron direction - Method based on:
407  // J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin'
408  // "A low energy bound atomic electron Compton scattering model for Geant4"
409  // NIMB, Vol. 338, 77-88, 2014.
410 
411  // Set constants and initialize scattering parameters
412 
413  const G4double vel_c = c_light / (m/s);
414  const G4double momentum_au_to_nat = halfpi* hbar_Planck / Bohr_radius / (kg*m/s);
415  const G4double e_mass_kg = electron_mass_c2 / c_squared / kg ;
416 
417  const G4int maxDopplerIterations = 1000;
418  G4double bindingE = 0.;
419  G4double pEIncident = photonEnergy0 ;
420  G4double pERecoil = -1.;
421  G4double eERecoil = -1.;
422  G4double e_alpha =0.;
423  G4double e_beta = 0.;
424 
425  G4double CE_emission_flag = 0.;
426  G4double ePAU = -1;
427  G4int shellIdx = 0;
428  G4double u_temp = 0;
429  G4double cosPhiE =0;
430  G4double sinThetaE =0;
431  G4double cosThetaE =0;
432  G4int iteration = 0;
433 
434  if (verboseLevel > 3) {
435  G4cout << "Started loop to sample photon energy and electron direction" << G4endl;
436  }
437 
438  do{
439 
440 
441  // ******************************************
442  // | Determine scatter photon energy |
443  // ******************************************
444 
445  do
446  {
447  iteration++;
448 
449 
450  // ********************************************
451  // | Sample bound electron information |
452  // ********************************************
453 
454  // Select shell based on shell occupancy
455 
456  shellIdx = shellData->SelectRandomShell(Z);
457  bindingE = shellData->BindingEnergy(Z,shellIdx)/MeV;
458 
459 
460  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
461  ePAU = profileData->RandomSelectMomentum(Z,shellIdx);
462 
463  // Convert to SI units
464  G4double ePSI = ePAU * momentum_au_to_nat;
465 
466  //Calculate bound electron velocity and normalise to natural units
467  u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)) )/vel_c;
468 
469  // Sample incident electron direction, amorphous material, to scattering photon scattering plane
470 
471  e_alpha = pi*G4UniformRand();
472  e_beta = twopi*G4UniformRand();
473 
474  // Total energy of system
475 
476  G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
477  G4double systemE = eEIncident + pEIncident;
478 
479 
480  G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
481  G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha));
482  G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
483  G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
484  G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident);
485  pERecoil = (numerator/denominator);
486  eERecoil = systemE - pERecoil;
487  CE_emission_flag = pEIncident - pERecoil;
488  } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
489 
490 // End of recalculation of photon energy with Doppler broadening
491 
492 
493 
494  // *******************************************************
495  // | Determine ejected Compton electron direction |
496  // *******************************************************
497 
498  // Calculate velocity of ejected Compton electron
499 
500  G4double a_temp = eERecoil / electron_mass_c2;
501  G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
502 
503  // Coefficients and terms from simulatenous equations
504 
505  G4double sinAlpha = std::sin(e_alpha);
506  G4double cosAlpha = std::cos(e_alpha);
507  G4double sinBeta = std::sin(e_beta);
508  G4double cosBeta = std::cos(e_beta);
509 
510  G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
511  G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
512 
513  G4double var_A = pERecoil*u_p_temp*sinTheta;
514  G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
515  G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta);
516 
517  G4double var_D1 = gamma*electron_mass_c2*pERecoil;
518  G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
519  G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil);
520  G4double var_D = var_D1*var_D2 + var_D3;
521 
522  G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha);
523  G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
524  G4double var_E = var_E1 - var_E2;
525 
526  G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
527  G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
528  G4double var_F = var_F1 - var_F2;
529 
530  G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
531 
532  // Two equations form a quadratic form of Wx^2 + Yx + Z = 0
533  // Coefficents and solution to quadratic
534 
535  G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
536  G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
537  G4double var_W = var_W1 + var_W2;
538 
539  G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C));
540 
541  G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
542  G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
543  G4double var_Z = var_Z1 + var_Z2;
544  G4double diff1 = var_Y*var_Y;
545  G4double diff2 = 4*var_W*var_Z;
546  G4double diff = diff1 - diff2;
547 
548 
549  // Check if diff is less than zero, if so ensure it is due to FPE
550 
551  //Confirm that diff less than zero is due FPE, i.e if abs of diff / diff1 and diff/ diff2 is less
552  //than 10^(-g4d_order), then set diff to zero
553 
554  if ((diff < 0.0) && (abs(diff / diff1) < g4d_limit) && (abs(diff / diff2) < g4d_limit) )
555  {
556  diff = 0.0;
557  }
558 
559 
560  // Plus and minus of quadratic
561  G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
562  G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
563 
564 
565  // Randomly sample one of the two possible solutions and determin theta angle of ejected Compton electron
566  G4double ThetaE = 0.;
567  G4double sol_select = G4UniformRand();
568 
569  if (sol_select < 0.5)
570  {
571  ThetaE = std::acos(X_p);
572  }
573  if (sol_select > 0.5)
574  {
575  ThetaE = std::acos(X_m);
576  }
577 
578  cosThetaE = std::cos(ThetaE);
579  sinThetaE = std::sin(ThetaE);
580  G4double Theta = std::acos(cosTheta);
581 
582  //Calculate electron Phi
583  G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE));
584  G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta));
585  G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
586  // Trigs
587  cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
588 
589  // End of calculation of ejection Compton electron direction
590 
591  //Fix for floating point errors
592 
593  } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1));
594 
595  // Revert to original if maximum number of iterations threshold has been reached
596  if (iteration >= maxDopplerIterations)
597  {
598  pERecoil = photonEnergy0 ;
599  bindingE = 0.;
600  dirx=0.0;
601  diry=0.0;
602  dirz=1.0;
603  }
604 
605  // Set "scattered" photon direction and energy
606 
607  G4ThreeVector photonDirection1(dirx,diry,dirz);
608  photonDirection1.rotateUz(photonDirection0);
609  photonPolarization1.rotateUz(photonDirection0);
610  if (pERecoil > 0.)
611  {
612  fParticleChange->SetProposedKineticEnergy(pERecoil) ;
613  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
614  fParticleChange->ProposePolarization(photonPolarization1);
615 
616  // Set ejected Compton electron direction and energy
617  G4double PhiE = std::acos(cosPhiE);
618  G4double eDirX = sinThetaE * std::cos(phi+PhiE);
619  G4double eDirY = sinThetaE * std::sin(phi+PhiE);
620  G4double eDirZ = cosThetaE;
621 
622  G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
623 
624  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
625  eDirection.rotateUz(photonDirection0);
626  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
627  eDirection,eKineticEnergy) ;
628  fvect->push_back(dp);
629 
630  }
631  else
632  {
633  fParticleChange->SetProposedKineticEnergy(0.);
634  fParticleChange->ProposeTrackStatus(fStopAndKill);
635  }
636 
637  // sample deexcitation
638  //
639 
640  if (verboseLevel > 3) {
641  G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl;
642  }
643 
644  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
645  G4int index = couple->GetIndex();
646  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
647  size_t nbefore = fvect->size();
648  G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
649  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
650  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
651  size_t nafter = fvect->size();
652  if(nafter > nbefore) {
653  for (size_t i=nbefore; i<nafter; ++i) {
654  //Check if there is enough residual energy
655  if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
656  {
657  //Ok, this is a valid secondary: keep it
658  bindingE -= ((*fvect)[i])->GetKineticEnergy();
659  }
660  else
661  {
662  //Invalid secondary: not enough energy to create it!
663  //Keep its energy in the local deposit
664  delete (*fvect)[i];
665  (*fvect)[i]=0;
666  }
667  }
668  }
669  }
670  }
671 
672  //This should never happen
673  if(bindingE < 0.0)
674  G4Exception("G4LowEPPolarizedComptonModel::SampleSecondaries()",
675  "em2051",FatalException,"Negative local energy deposit");
676 
677  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
678 
679 }
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
CLHEP::h_Planck
static constexpr double h_Planck
Definition: PhysicalConstants.h:63
G4VAtomDeexcitation::CheckDeexcitationActiveRegion
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:262
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
CLHEP::Bohr_radius
static constexpr double Bohr_radius
Definition: PhysicalConstants.h:99
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:178
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
CLHEP::Hep3Vector::isOrthogonal
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:237
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
Definition: G4ParticleChangeForGamma.hh:147
G4ShellData::SelectRandomShell
G4int SelectRandomShell(G4int Z) const
Definition: G4ShellData.cc:363
CLHEP::Hep3Vector::howOrthogonal
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:219
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
s
const XML_Char * s
Definition: expat.h:262
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
m
static constexpr double m
Definition: G4SIunits.hh:129
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4ShellData::BindingEnergy
G4double BindingEnergy(G4int Z, G4int shellIndex) const
Definition: G4ShellData.cc:166
cm
static constexpr double cm
Definition: G4SIunits.hh:119
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4ParticleChangeForGamma::ProposePolarization
void ProposePolarization(const G4ThreeVector &dir)
Definition: G4ParticleChangeForGamma.hh:166
kg
static constexpr double kg
Definition: G4SIunits.hh:182
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
CLHEP::c_squared
static constexpr double c_squared
Definition: PhysicalConstants.h:56
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
CLHEP::c_light
static constexpr double c_light
Definition: PhysicalConstants.h:55
G4DynamicParticle::GetPolarization
const G4ThreeVector & GetPolarization() const
FatalException
Definition: G4ExceptionSeverity.hh:60
fStopAndKill
Definition: G4TrackStatus.hh:56
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:129
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
pi
static constexpr double pi
Definition: G4SIunits.hh:75
halfpi
static constexpr double halfpi
Definition: G4SIunits.hh:77
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:274
CLHEP::hbar_Planck
static constexpr double hbar_Planck
Definition: PhysicalConstants.h:64
CLHEP::Hep3Vector::mag
double mag() const
G4DopplerProfile::RandomSelectMomentum
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
Definition: G4DopplerProfile.cc:227
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
G4VEmModel::SelectRandomAtom
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:541
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150

Here is the call graph for this function:

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