Geant4  10.02.p01
G4LowEnergyCompton.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id$
27 // GEANT4 tag $Name: $
28 //
29 // Author: A. Forti
30 // Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
31 //
32 // History:
33 // --------
34 // Added Livermore data table construction methods A. Forti
35 // Modified BuildMeanFreePath to read new data tables A. Forti
36 // Modified PostStepDoIt to insert sampling with EPDL97 data A. Forti
37 // Added SelectRandomAtom A. Forti
38 // Added map of the elements A. Forti
39 // 24.04.2001 V.Ivanchenko - Remove RogueWave
40 // 06.08.2001 MGP - Revised according to a design iteration
41 // 22.01.2003 V.Ivanchenko - Cut per region
42 // 10.03.2003 V.Ivanchenko - Remove CutPerMaterial warning
43 // 24.04.2003 V.Ivanchenko - Cut per region mfpt
44 //
45 // -------------------------------------------------------------------
46 
47 #include "G4LowEnergyCompton.hh"
48 #include "Randomize.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4ParticleDefinition.hh"
52 #include "G4Track.hh"
53 #include "G4Step.hh"
54 #include "G4ForceCondition.hh"
55 #include "G4Gamma.hh"
56 #include "G4Electron.hh"
57 #include "G4DynamicParticle.hh"
58 #include "G4VParticleChange.hh"
59 #include "G4ThreeVector.hh"
60 #include "G4EnergyLossTables.hh"
63 #include "G4RDVEMDataSet.hh"
65 #include "G4RDVDataSetAlgorithm.hh"
67 #include "G4RDVRangeTest.hh"
68 #include "G4RDRangeTest.hh"
69 #include "G4RDRangeNoTest.hh"
70 #include "G4MaterialCutsCouple.hh"
71 
72 
74  : G4VDiscreteProcess(processName),
75  lowEnergyLimit(250*eV),
76  highEnergyLimit(100*GeV),
77  intrinsicLowEnergyLimit(10*eV),
78  intrinsicHighEnergyLimit(100*GeV)
79 {
82  {
83  G4Exception("G4LowEnergyCompton::G4LowEnergyCompton()",
84  "OutOfRange", FatalException,
85  "Energy outside intrinsic process validity range!");
86  }
87 
89 
90  G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation;
91  G4String scatterFile = "comp/ce-sf-";
92  scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation, 1., 1.);
93  scatterFunctionData->LoadData(scatterFile);
94 
96 
98 
99  // For Doppler broadening
101 
102  if (verboseLevel > 0)
103  {
104  G4cout << GetProcessName() << " is created " << G4endl
105  << "Energy range: "
106  << lowEnergyLimit / keV << " keV - "
107  << highEnergyLimit / GeV << " GeV"
108  << G4endl;
109  }
110 }
111 
113 {
114  delete meanFreePathTable;
115  delete crossSectionHandler;
116  delete scatterFunctionData;
117  delete rangeTest;
118 }
119 
121 {
122 
124  G4String crossSectionFile = "comp/ce-cs-";
125  crossSectionHandler->LoadData(crossSectionFile);
126 
127  delete meanFreePathTable;
129 
130  // For Doppler broadening
131  G4String file = "/doppler/shell-doppler";
132  shellData.LoadData(file);
133 }
134 
136  const G4Step& aStep)
137 {
138  // The scattered gamma energy is sampled according to Klein - Nishina formula.
139  // then accepted or rejected depending on the Scattering Function multiplied
140  // by factor from Klein - Nishina formula.
141  // Expression of the angular distribution as Klein Nishina
142  // angular and energy distribution and Scattering fuctions is taken from
143  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
144  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
145  // data are interpolated while in the article they are fitted.
146  // Reference to the article is from J. Stepanek New Photon, Positron
147  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
148  // TeV (draft).
149  // The random number techniques of Butcher & Messel are used
150  // (Nucl Phys 20(1960),15).
151 
152  aParticleChange.Initialize(aTrack);
153 
154  // Dynamic particle quantities
155  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
156  G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
157 
158  if (photonEnergy0 <= lowEnergyLimit)
159  {
163  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
164  }
165 
166  G4double e0m = photonEnergy0 / electron_mass_c2 ;
167  G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
168  G4double epsilon0 = 1. / (1. + 2. * e0m);
169  G4double epsilon0Sq = epsilon0 * epsilon0;
170  G4double alpha1 = -std::log(epsilon0);
171  G4double alpha2 = 0.5 * (1. - epsilon0Sq);
172  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
173 
174  // Select randomly one element in the current material
175  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
176  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
177 
178  // Sample the energy of the scattered photon
180  G4double epsilonSq;
181  G4double oneCosT;
182  G4double sinT2;
183  G4double gReject;
184  do
185  {
186  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
187  {
188  epsilon = std::exp(-alpha1 * G4UniformRand()); // std::pow(epsilon0,G4UniformRand())
189  epsilonSq = epsilon * epsilon;
190  }
191  else
192  {
193  epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
194  epsilon = std::sqrt(epsilonSq);
195  }
196 
197  oneCosT = (1. - epsilon) / ( epsilon * e0m);
198  sinT2 = oneCosT * (2. - oneCosT);
199  G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
200  G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
201  gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
202 
203  } while(gReject < G4UniformRand()*Z);
204 
205  G4double cosTheta = 1. - oneCosT;
206  G4double sinTheta = std::sqrt (sinT2);
207  G4double phi = twopi * G4UniformRand() ;
208  G4double dirX = sinTheta * std::cos(phi);
209  G4double dirY = sinTheta * std::sin(phi);
210  G4double dirZ = cosTheta ;
211 
212  // Doppler broadening - Method based on:
213  // Y. Namito, S. Ban and H. Hirayama,
214  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
215  // NIM A 349, pp. 489-494, 1994
216 
217  // Maximum number of sampling iterations
218  G4int maxDopplerIterations = 1000;
219  G4double bindingE = 0.;
220  G4double photonEoriginal = epsilon * photonEnergy0;
221  G4double photonE = -1.;
222  G4int iteration = 0;
223  G4double eMax = photonEnergy0;
224  do
225  {
226  iteration++;
227  // Select shell based on shell occupancy
228  G4int shell = shellData.SelectRandomShell(Z);
229  bindingE = shellData.BindingEnergy(Z,shell);
230 
231  eMax = photonEnergy0 - bindingE;
232 
233  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
234  G4double pSample = profileData.RandomSelectMomentum(Z,shell);
235  // Rescale from atomic units
236  G4double pDoppler = pSample * fine_structure_const;
237  G4double pDoppler2 = pDoppler * pDoppler;
238  G4double var2 = 1. + oneCosT * e0m;
239  G4double var3 = var2*var2 - pDoppler2;
240  G4double var4 = var2 - pDoppler2 * cosTheta;
241  G4double var = var4*var4 - var3 + pDoppler2 * var3;
242  if (var > 0.)
243  {
244  G4double varSqrt = std::sqrt(var);
245  G4double scale = photonEnergy0 / var3;
246  // Random select either root
247  if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
248  else photonE = (var4 + varSqrt) * scale;
249  }
250  else
251  {
252  photonE = -1.;
253  }
254  } while ( iteration <= maxDopplerIterations &&
255  (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
256 
257  // End of recalculation of photon energy with Doppler broadening
258  // Revert to original if maximum number of iterations threshold has been reached
259  if (iteration >= maxDopplerIterations)
260  {
261  photonE = photonEoriginal;
262  bindingE = 0.;
263  }
264 
265  // Update G4VParticleChange for the scattered photon
266 
267  G4ThreeVector photonDirection1(dirX,dirY,dirZ);
268  photonDirection1.rotateUz(photonDirection0);
269  aParticleChange.ProposeMomentumDirection(photonDirection1);
270 
271  G4double photonEnergy1 = photonE;
272  //G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
273 
274  if (photonEnergy1 > 0.)
275  {
276  aParticleChange.ProposeEnergy(photonEnergy1) ;
277  }
278  else
279  {
282  }
283 
284  // Kinematics of the scattered electron
285  G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
286  G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
287 
288  G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2;
289  G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2;
290  G4double sinThetaE = -1.;
291  G4double cosThetaE = 0.;
292  if (electronP2 > 0.)
293  {
294  cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2);
295  sinThetaE = -1. * std::sqrt(1. - cosThetaE * cosThetaE);
296  }
297 
298  G4double eDirX = sinThetaE * std::cos(phi);
299  G4double eDirY = sinThetaE * std::sin(phi);
300  G4double eDirZ = cosThetaE;
301 
302  // Generate the electron only if with large enough range w.r.t. cuts and safety
303 
304  G4double safety = aStep.GetPostStepPoint()->GetSafety();
305 
306  if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
307  {
308  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
309  eDirection.rotateUz(photonDirection0);
310 
311  G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),eDirection,eKineticEnergy) ;
313  aParticleChange.AddSecondary(electron);
314  // Binding energy deposited locally
316  }
317  else
318  {
320  aParticleChange.ProposeLocalEnergyDeposit(eKineticEnergy + bindingE);
321  }
322 
323  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
324 }
325 
327 {
328  return ( &particle == G4Gamma::Gamma() );
329 }
330 
332  G4double, // previousStepSize
334 {
336  G4double energy = photon->GetKineticEnergy();
337  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
338  size_t materialIndex = couple->GetIndex();
339 
340  G4double meanFreePath;
341  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
342  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
343  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
344  return meanFreePath;
345 }
static const double cm
Definition: G4SIunits.hh:118
void SetOccupancyData()
G4int verboseLevel
Definition: G4VProcess.hh:368
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
G4LowEnergyCompton(const G4String &processName="LowEnCompton")
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4RDVEMDataSet * meanFreePathTable
G4double BindingEnergy(G4int Z, G4int shellIndex) const
G4RDVEMDataSet * BuildMeanFreePathForMaterials(const G4DataVector *energyCuts=0)
int G4int
Definition: G4Types.hh:78
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void LoadData(const G4String &fileName)
void BuildPhysicsTable(const G4ParticleDefinition &definition)
G4RDVRangeTest * rangeTest
G4RDVCrossSectionHandler * crossSectionHandler
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
static const double twopi
Definition: G4SIunits.hh:75
G4bool IsApplicable(const G4ParticleDefinition &definition)
static const double GeV
Definition: G4SIunits.hh:214
Definition: G4Step.hh:76
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4RDVEMDataSet * scatterFunctionData
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
virtual G4bool Escape(const G4ParticleDefinition *particle, const G4MaterialCutsCouple *couple, G4double energy, G4double safety) const =0
virtual void Initialize(const G4Track &)
static const double eV
Definition: G4SIunits.hh:212
G4double energy(const ThreeVector &p, const G4double m)
void SetNumberOfSecondaries(G4int totSecondaries)
G4StepPoint * GetPostStepPoint() const
const G4double intrinsicLowEnergyLimit
const G4double x[NPOINTSGL]
G4int SelectRandomAtom(const G4MaterialCutsCouple *couple, G4double e) const
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4double GetSafety() const
void AddSecondary(G4Track *aSecondary)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &step)
void LoadData(const G4String &dataFile)
static const double keV
Definition: G4SIunits.hh:213
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ForceCondition
G4ThreeVector G4ParticleMomentum
#define DBL_MAX
Definition: templates.hh:83
G4RDDopplerProfile profileData
G4int SelectRandomShell(G4int Z) const
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
double epsilon(double density, double temperature)
virtual G4bool LoadData(const G4String &fileName)=0
const G4double intrinsicHighEnergyLimit