Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Protected Attributes
G4LowEPComptonModel Class Reference

#include <G4LowEPComptonModel.hh>

Inheritance diagram for G4LowEPComptonModel:
G4VEmModel

Public Member Functions

 G4LowEPComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="LowEPComptonModel")
 
virtual ~G4LowEPComptonModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
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 InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
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 CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, 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 *)
 
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
 

Protected Attributes

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

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 *)
 

Detailed Description

Definition at line 78 of file G4LowEPComptonModel.hh.

Constructor & Destructor Documentation

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

Definition at line 87 of file G4LowEPComptonModel.cc.

References python.hepunit::eV, G4cout, G4endl, python.hepunit::GeV, and G4VEmModel::SetDeexcitationFlag().

89  :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
90  scatterFunctionData(0),crossSectionHandler(0),fAtomDeexcitation(0)
91 {
92  lowEnergyLimit = 250 * eV;
93  highEnergyLimit = 100 * GeV;
94 
95  verboseLevel=0 ;
96  // Verbosity scale:
97  // 0 = nothing
98  // 1 = warning for energy non-conservation
99  // 2 = details of energy budget
100  // 3 = calculation of cross sections, file openings, sampling of atoms
101  // 4 = entering in methods
102 
103  if( verboseLevel>0 ) {
104  G4cout << "Low energy photon Compton model is constructed " << G4endl
105  << "Energy range: "
106  << lowEnergyLimit / eV << " eV - "
107  << highEnergyLimit / GeV << " GeV"
108  << G4endl;
109  }
110 
111  //Mark this model as "applicable" for atomic deexcitation
112  SetDeexcitationFlag(true);
113 
114 }
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:65
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4LowEPComptonModel::fParticleChange
G4ParticleChangeForGamma * fParticleChange
Definition: G4LowEPComptonModel.hh:105
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:739
G4LowEPComptonModel::~G4LowEPComptonModel ( )
virtual

Definition at line 118 of file G4LowEPComptonModel.cc.

119 {
120  delete crossSectionHandler;
121  delete scatterFunctionData;
122 }

Member Function Documentation

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

Reimplemented from G4VEmModel.

Definition at line 179 of file G4LowEPComptonModel.cc.

References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.

184 {
185  if (verboseLevel > 3) {
186  G4cout << "Calling ComputeCrossSectionPerAtom() of G4LowEPComptonModel" << G4endl;
187  }
188  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) { return 0.0; }
189 
190  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
191  return cs;
192 }
G4int
int G4int
Definition: G4Types.hh:78
G4VCrossSectionHandler::FindValue
G4double FindValue(G4int Z, G4double e) const
Definition: G4VCrossSectionHandler.cc:379
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
void G4LowEPComptonModel::Initialise ( const G4ParticleDefinition *  particle,
const G4DataVector &  cuts 
)
virtual

Implements G4VEmModel.

Definition at line 126 of file G4LowEPComptonModel.cc.

References G4LossTableManager::AtomDeexcitation(), G4VCrossSectionHandler::Clear(), python.hepunit::eV, fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), python.hepunit::GeV, G4VEmModel::HighEnergyLimit(), G4VEmModel::InitialiseElementSelectors(), G4LossTableManager::Instance(), G4ShellData::LoadData(), G4VEMDataSet::LoadData(), G4VCrossSectionHandler::LoadData(), G4VEmModel::LowEnergyLimit(), and G4ShellData::SetOccupancyData().

128 {
129  if (verboseLevel > 2) {
130  G4cout << "Calling G4LowEPComptonModel::Initialise()" << G4endl;
131  }
132 
133  if (crossSectionHandler)
134  {
135  crossSectionHandler->Clear();
136  delete crossSectionHandler;
137  }
138  delete scatterFunctionData;
139 
140  // Reading of data files - all materials are read
141  crossSectionHandler = new G4CrossSectionHandler;
142  G4String crossSectionFile = "comp/ce-cs-";
143  crossSectionHandler->LoadData(crossSectionFile);
144 
145  G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
146  G4String scatterFile = "comp/ce-sf-";
147  scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
148  scatterFunctionData->LoadData(scatterFile);
149 
150  // For Doppler broadening
151  shellData.SetOccupancyData();
152  G4String file = "/doppler/shell-doppler";
153  shellData.LoadData(file);
154 
155  InitialiseElementSelectors(particle,cuts);
156 
157  if (verboseLevel > 2) {
158  G4cout << "Loaded cross section files for low energy photon Compton model" << G4endl;
159  }
160 
161  if(isInitialised) { return; }
162  isInitialised = true;
163 
164  fParticleChange = GetParticleChangeForGamma();
165 
166  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
167 
168  if( verboseLevel>0 ) {
169  G4cout << "Low energy photon Compton model is initialized " << G4endl
170  << "Energy range: "
171  << LowEnergyLimit() / eV << " eV - "
172  << HighEnergyLimit() / GeV << " GeV"
173  << G4endl;
174  }
175 }
G4ShellData::SetOccupancyData
void SetOccupancyData()
Definition: G4ShellData.hh:70
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:112
G4VEmModel::InitialiseElementSelectors
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:135
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
G4ShellData::LoadData
void LoadData(const G4String &fileName)
Definition: G4ShellData.cc:234
G4cout
G4GLOB_DLL std::ostream G4cout
G4VCrossSectionHandler::LoadData
void LoadData(const G4String &dataFile)
Definition: G4VCrossSectionHandler.cc:187
G4VEMDataSet::LoadData
virtual G4bool LoadData(const G4String &fileName)=0
G4endl
#define G4endl
Definition: G4ios.hh:61
G4LowEPComptonModel::fParticleChange
G4ParticleChangeForGamma * fParticleChange
Definition: G4LowEPComptonModel.hh:105
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.cc:1300
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:121
void G4LowEPComptonModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  fvect,
const G4MaterialCutsCouple *  couple,
const G4DynamicParticle *  aDynamicGamma,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 201 of file G4LowEPComptonModel.cc.

References G4ShellData::BindingEnergy(), python.hepunit::Bohr_radius, python.hepunit::c_light, python.hepunit::c_squared, G4VAtomDeexcitation::CheckDeexcitationActiveRegion(), python.hepunit::cm, G4Electron::Electron(), python.hepunit::electron_mass_c2, G4VEMDataSet::FindValue(), fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4VAtomDeexcitation::GenerateParticles(), G4VAtomDeexcitation::GetAtomicShell(), G4DynamicParticle::GetDefinition(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4Material::GetName(), G4Element::GetZ(), python.hepunit::h_Planck, python.hepunit::halfpi, python.hepunit::hbar_Planck, python.hepunit::kg, python.hepunit::m, python.hepunit::MeV, python.hepunit::pi, G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4DopplerProfile::RandomSelectMomentum(), CLHEP::Hep3Vector::rotateUz(), G4VEmModel::SelectRandomAtom(), G4ShellData::SelectRandomShell(), G4ParticleChangeForGamma::SetProposedKineticEnergy(), python.hepunit::twopi, and test::x.

205 {
206 
207  // The scattered gamma energy is sampled according to Klein - Nishina formula.
208  // then accepted or rejected depending on the Scattering Function multiplied
209  // by factor from Klein - Nishina formula.
210  // Expression of the angular distribution as Klein Nishina
211  // angular and energy distribution and Scattering fuctions is taken from
212  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
213  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
214  // data are interpolated while in the article they are fitted.
215  // Reference to the article is from J. Stepanek New Photon, Positron
216  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
217  // TeV (draft).
218  // The random number techniques of Butcher & Messel are used
219  // (Nucl Phys 20(1960),15).
220 
221  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV;
222 
223  if (verboseLevel > 3) {
224  G4cout << "G4LowEPComptonModel::SampleSecondaries() E(MeV)= "
225  << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
226  << G4endl;
227  }
228 
229  // low-energy gamma is absorpted by this process
230  if (photonEnergy0 <= lowEnergyLimit)
231  {
232  fParticleChange->ProposeTrackStatus(fStopAndKill);
233  fParticleChange->SetProposedKineticEnergy(0.);
234  fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
235  return ;
236  }
237 
238  G4double e0m = photonEnergy0 / electron_mass_c2 ;
239  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
240 
241  // Select randomly one element in the current material
242  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
243  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
244  G4int Z = (G4int)elm->GetZ();
245 
246  G4double LowEPCepsilon0 = 1. / (1. + 2. * e0m);
247  G4double LowEPCepsilon0Sq = LowEPCepsilon0 * LowEPCepsilon0;
248  G4double alpha1 = -std::log(LowEPCepsilon0);
249  G4double alpha2 = 0.5 * (1. - LowEPCepsilon0Sq);
250 
251  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
252 
253  // Sample the energy of the scattered photon
254  G4double LowEPCepsilon;
255  G4double LowEPCepsilonSq;
256  G4double oneCosT;
257  G4double sinT2;
258  G4double gReject;
259 
260  do
261  {
262  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
263  {
264  LowEPCepsilon = std::exp(-alpha1 * G4UniformRand());
265  LowEPCepsilonSq = LowEPCepsilon * LowEPCepsilon;
266  }
267  else
268  {
269  LowEPCepsilonSq = LowEPCepsilon0Sq + (1. - LowEPCepsilon0Sq) * G4UniformRand();
270  LowEPCepsilon = std::sqrt(LowEPCepsilonSq);
271  }
272 
273  oneCosT = (1. - LowEPCepsilon) / ( LowEPCepsilon * e0m);
274  sinT2 = oneCosT * (2. - oneCosT);
275  G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
276  G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
277  gReject = (1. - LowEPCepsilon * sinT2 / (1. + LowEPCepsilonSq)) * scatteringFunction;
278 
279  } while(gReject < G4UniformRand()*Z);
280 
281  G4double cosTheta = 1. - oneCosT;
282  G4double sinTheta = std::sqrt(sinT2);
283  G4double phi = twopi * G4UniformRand();
284  G4double dirx = sinTheta * std::cos(phi);
285  G4double diry = sinTheta * std::sin(phi);
286  G4double dirz = cosTheta ;
287 
288 
289  // Scatter photon energy and Compton electron direction - Method based on:
290  // J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin'
291  // "A low energy bound atomic electron Compton scattering model for Geant4"
292  // NIMA ISSUE, PG, 2013
293 
294  // Set constants and initialize scattering parameters
295 
296  const G4double vel_c = c_light / (m/s);
297  const G4double momentum_au_to_nat = halfpi* hbar_Planck / Bohr_radius / (kg*m/s);
298  const G4double e_mass_kg = electron_mass_c2 / c_squared / kg ;
299 
300  const G4int maxDopplerIterations = 1000;
301  G4double bindingE = 0.;
302  G4double pEIncident = photonEnergy0 ;
303  G4double pERecoil = -1.;
304  G4double eERecoil = -1.;
305  G4double e_alpha =0.;
306  G4double e_beta = 0.;
307 
308  G4double CE_emission_flag = 0.;
309  G4double ePAU = -1;
310  G4int shellIdx = 0;
311  G4double u_temp = 0;
312  G4double cosPhiE =0;
313  G4double sinThetaE =0;
314  G4double cosThetaE =0;
315  G4int iteration = 0;
316  do{
317 
318 
319  // ******************************************
320  // | Determine scatter photon energy |
321  // ******************************************
322 
323  do
324  {
325  iteration++;
326 
327 
328  // ********************************************
329  // | Sample bound electron information |
330  // ********************************************
331 
332  // Select shell based on shell occupancy
333 
334  shellIdx = shellData.SelectRandomShell(Z);
335  bindingE = shellData.BindingEnergy(Z,shellIdx)/MeV;
336 
337 
338  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
339  ePAU = profileData.RandomSelectMomentum(Z,shellIdx);
340 
341  // Convert to SI units
342  G4double ePSI = ePAU * momentum_au_to_nat;
343 
344  //Calculate bound electron velocity and normalise to natural units
345  u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)) )/vel_c;
346 
347  // Sample incident electron direction, amorphous material, to scattering photon scattering plane
348 
349  e_alpha = pi*G4UniformRand();
350  e_beta = twopi*G4UniformRand();
351 
352  // Total energy of system
353 
354  G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
355  G4double systemE = eEIncident + pEIncident;
356 
357 
358  G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
359  G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha));
360  G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
361  G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
362  G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident);
363  pERecoil = (numerator/denominator);
364  eERecoil = systemE - pERecoil;
365  CE_emission_flag = pEIncident - pERecoil;
366  } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
367 
368 
369 
370  // End of recalculation of photon energy with Doppler broadening
371 
372 
373 
374  // *******************************************************
375  // | Determine ejected Compton electron direction |
376  // *******************************************************
377 
378  // Calculate velocity of ejected Compton electron
379 
380  G4double a_temp = eERecoil / electron_mass_c2;
381  G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
382 
383  // Coefficients and terms from simulatenous equations
384 
385  G4double sinAlpha = std::sin(e_alpha);
386  G4double cosAlpha = std::cos(e_alpha);
387  G4double sinBeta = std::sin(e_beta);
388  G4double cosBeta = std::cos(e_beta);
389 
390  G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
391  G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
392 
393  G4double var_A = pERecoil*u_p_temp*sinTheta;
394  G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
395  G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta);
396 
397  G4double var_D1 = gamma*electron_mass_c2*pERecoil;
398  G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
399  G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil);
400  G4double var_D = var_D1*var_D2 + var_D3;
401 
402  G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha);
403  G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
404  G4double var_E = var_E1 - var_E2;
405 
406  G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
407  G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
408  G4double var_F = var_F1 - var_F2;
409 
410  G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
411 
412  // Two equations form a quadratic form of Wx^2 + Yx + Z = 0
413  // Coefficents and solution to quadratic
414 
415  G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
416  G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
417  G4double var_W = var_W1 + var_W2;
418 
419  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));
420 
421  G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
422  G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
423  G4double var_Z = var_Z1 + var_Z2;
424  G4double diff1 = var_Y*var_Y;
425  G4double diff2 = 4*var_W*var_Z;
426  G4double diff = diff1 - diff2;
427 
428 
429  // Check if diff is less than zero, if so ensure it is due to FPE
430 
431  //Determine number of digits (in decimal base) that G4double can accurately represent
432  G4double g4d_order = G4double(numeric_limits<G4double>::digits10);
433  G4double g4d_limit = std::pow(10.,-g4d_order);
434  //Confirm that diff less than zero is due FPE, i.e if abs of diff / diff1 and diff/ diff2 is less
435  //than 10^(-g4d_order), then set diff to zero
436 
437  if ((diff < 0.0) && (abs(diff / diff1) < g4d_limit) && (abs(diff / diff2) < g4d_limit) )
438  {
439  diff = 0.0;
440  }
441 
442 
443  // Plus and minus of quadratic
444  G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
445  G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
446 
447 
448  // Randomly sample one of the two possible solutions and determin theta angle of ejected Compton electron
449  G4double ThetaE = 0.;
450  G4double sol_select = G4UniformRand();
451 
452  if (sol_select < 0.5)
453  {
454  ThetaE = std::acos(X_p);
455  }
456  if (sol_select > 0.5)
457  {
458  ThetaE = std::acos(X_m);
459  }
460 
461  cosThetaE = std::cos(ThetaE);
462  sinThetaE = std::sin(ThetaE);
463  G4double Theta = std::acos(cosTheta);
464 
465  //Calculate electron Phi
466  G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE));
467  G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta));
468  G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
469  // Trigs
470  cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
471 
472  // End of calculation of ejection Compton electron direction
473 
474  //Fix for floating point errors
475 
476  } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1));
477 
478  // Revert to original if maximum number of iterations threshold has been reached
479 
480 
481  if (iteration >= maxDopplerIterations)
482  {
483  pERecoil = photonEnergy0 ;
484  bindingE = 0.;
485  dirx=0.0;
486  diry=0.0;
487  dirz=1.0;
488  }
489 
490  // Set "scattered" photon direction and energy
491 
492  G4ThreeVector photonDirection1(dirx,diry,dirz);
493  photonDirection1.rotateUz(photonDirection0);
494  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
495 
496  if (pERecoil > 0.)
497  {
498  fParticleChange->SetProposedKineticEnergy(pERecoil) ;
499 
500  // Set ejected Compton electron direction and energy
501  G4double PhiE = std::acos(cosPhiE);
502  G4double eDirX = sinThetaE * std::cos(phi+PhiE);
503  G4double eDirY = sinThetaE * std::sin(phi+PhiE);
504  G4double eDirZ = cosThetaE;
505 
506  G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
507 
508  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
509  eDirection.rotateUz(photonDirection0);
510  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
511  eDirection,eKineticEnergy) ;
512  fvect->push_back(dp);
513 
514  }
515  else
516  {
517  fParticleChange->SetProposedKineticEnergy(0.);
518  fParticleChange->ProposeTrackStatus(fStopAndKill);
519  }
520 
521 
522 
523 
524  // sample deexcitation
525 
526  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
527  G4int index = couple->GetIndex();
528  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
529  size_t nbefore = fvect->size();
530  G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
531  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
532  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
533  size_t nafter = fvect->size();
534  if(nafter > nbefore) {
535  for (size_t i=nbefore; i<nafter; ++i) {
536  bindingE -= ((*fvect)[i])->GetKineticEnergy();
537  }
538  }
539  }
540  }
541  if(bindingE < 0.0) { bindingE = 0.0; }
542  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
543 
544 }
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4VAtomDeexcitation::CheckDeexcitationActiveRegion
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:276
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
python.hepunit.h_Planck
float h_Planck
Definition: hepunit.py:263
s
const XML_Char * s
Definition: include/expat.h:262
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:176
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
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
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4ShellData::BindingEnergy
G4double BindingEnergy(G4int Z, G4int shellIndex) const
Definition: G4ShellData.cc:166
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
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
G4LowEPComptonModel::fParticleChange
G4ParticleChangeForGamma * fParticleChange
Definition: G4LowEPComptonModel.hh:105
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:289
G4DopplerProfile::RandomSelectMomentum
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
Definition: G4DopplerProfile.cc:227
python.hepunit.c_light
float c_light
Definition: hepunit.py:257
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:510
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150

Field Documentation

G4ParticleChangeForGamma* G4LowEPComptonModel::fParticleChange
protected

Definition at line 105 of file G4LowEPComptonModel.hh.

Referenced by Initialise(), and SampleSecondaries().

The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:40 for Geant4.10 by   doxygen 1.8.7