ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEnergyCompton.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LowEnergyCompton.cc
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 //
27 // Author: A. Forti
28 // Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 //
30 // History:
31 // --------
32 // Added Livermore data table construction methods A. Forti
33 // Modified BuildMeanFreePath to read new data tables A. Forti
34 // Modified PostStepDoIt to insert sampling with EPDL97 data A. Forti
35 // Added SelectRandomAtom A. Forti
36 // Added map of the elements A. Forti
37 // 24.04.2001 V.Ivanchenko - Remove RogueWave
38 // 06.08.2001 MGP - Revised according to a design iteration
39 // 22.01.2003 V.Ivanchenko - Cut per region
40 // 10.03.2003 V.Ivanchenko - Remove CutPerMaterial warning
41 // 24.04.2003 V.Ivanchenko - Cut per region mfpt
42 //
43 // -------------------------------------------------------------------
44 
45 #include "G4LowEnergyCompton.hh"
46 #include "Randomize.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4ParticleDefinition.hh"
50 #include "G4Track.hh"
51 #include "G4Step.hh"
52 #include "G4ForceCondition.hh"
53 #include "G4Gamma.hh"
54 #include "G4Electron.hh"
55 #include "G4DynamicParticle.hh"
56 #include "G4VParticleChange.hh"
57 #include "G4ThreeVector.hh"
58 #include "G4EnergyLossTables.hh"
61 #include "G4RDVEMDataSet.hh"
63 #include "G4RDVDataSetAlgorithm.hh"
65 #include "G4RDVRangeTest.hh"
66 #include "G4RDRangeTest.hh"
67 #include "G4RDRangeNoTest.hh"
68 #include "G4MaterialCutsCouple.hh"
69 
70 
72  : G4VDiscreteProcess(processName),
73  lowEnergyLimit(250*eV),
74  highEnergyLimit(100*GeV),
75  intrinsicLowEnergyLimit(10*eV),
76  intrinsicHighEnergyLimit(100*GeV)
77 {
80  {
81  G4Exception("G4LowEnergyCompton::G4LowEnergyCompton()",
82  "OutOfRange", FatalException,
83  "Energy outside intrinsic process validity range!");
84  }
85 
87 
88  G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation;
89  G4String scatterFile = "comp/ce-sf-";
90  scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation, 1., 1.);
91  scatterFunctionData->LoadData(scatterFile);
92 
94 
96 
97  // For Doppler broadening
99 
100  if (verboseLevel > 0)
101  {
102  G4cout << GetProcessName() << " is created " << G4endl
103  << "Energy range: "
104  << lowEnergyLimit / keV << " keV - "
105  << highEnergyLimit / GeV << " GeV"
106  << G4endl;
107  }
108 }
109 
111 {
112  delete meanFreePathTable;
113  delete crossSectionHandler;
114  delete scatterFunctionData;
115  delete rangeTest;
116 }
117 
119 {
120 
122  G4String crossSectionFile = "comp/ce-cs-";
123  crossSectionHandler->LoadData(crossSectionFile);
124 
125  delete meanFreePathTable;
127 
128  // For Doppler broadening
129  G4String file = "/doppler/shell-doppler";
130  shellData.LoadData(file);
131 }
132 
134  const G4Step& aStep)
135 {
136  // The scattered gamma energy is sampled according to Klein - Nishina formula.
137  // then accepted or rejected depending on the Scattering Function multiplied
138  // by factor from Klein - Nishina formula.
139  // Expression of the angular distribution as Klein Nishina
140  // angular and energy distribution and Scattering fuctions is taken from
141  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
142  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
143  // data are interpolated while in the article they are fitted.
144  // Reference to the article is from J. Stepanek New Photon, Positron
145  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
146  // TeV (draft).
147  // The random number techniques of Butcher & Messel are used
148  // (Nucl Phys 20(1960),15).
149 
150  aParticleChange.Initialize(aTrack);
151 
152  // Dynamic particle quantities
153  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
154  G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
155 
156  if (photonEnergy0 <= lowEnergyLimit)
157  {
161  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
162  }
163 
164  G4double e0m = photonEnergy0 / electron_mass_c2 ;
165  G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
166  G4double epsilon0 = 1. / (1. + 2. * e0m);
167  G4double epsilon0Sq = epsilon0 * epsilon0;
168  G4double alpha1 = -std::log(epsilon0);
169  G4double alpha2 = 0.5 * (1. - epsilon0Sq);
170  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
171 
172  // Select randomly one element in the current material
173  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
174  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
175 
176  // Sample the energy of the scattered photon
178  G4double epsilonSq;
179  G4double oneCosT;
180  G4double sinT2;
181  G4double gReject;
182  do
183  {
184  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
185  {
186  epsilon = std::exp(-alpha1 * G4UniformRand()); // std::pow(epsilon0,G4UniformRand())
187  epsilonSq = epsilon * epsilon;
188  }
189  else
190  {
191  epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
192  epsilon = std::sqrt(epsilonSq);
193  }
194 
195  oneCosT = (1. - epsilon) / ( epsilon * e0m);
196  sinT2 = oneCosT * (2. - oneCosT);
197  G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
198  G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
199  gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
200 
201  } while(gReject < G4UniformRand()*Z);
202 
203  G4double cosTheta = 1. - oneCosT;
204  G4double sinTheta = std::sqrt (sinT2);
206  G4double dirX = sinTheta * std::cos(phi);
207  G4double dirY = sinTheta * std::sin(phi);
208  G4double dirZ = cosTheta ;
209 
210  // Doppler broadening - Method based on:
211  // Y. Namito, S. Ban and H. Hirayama,
212  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
213  // NIM A 349, pp. 489-494, 1994
214 
215  // Maximum number of sampling iterations
216  G4int maxDopplerIterations = 1000;
217  G4double bindingE = 0.;
218  G4double photonEoriginal = epsilon * photonEnergy0;
219  G4double photonE = -1.;
220  G4int iteration = 0;
221  G4double eMax = photonEnergy0;
222  do
223  {
224  iteration++;
225  // Select shell based on shell occupancy
226  G4int shell = shellData.SelectRandomShell(Z);
227  bindingE = shellData.BindingEnergy(Z,shell);
228 
229  eMax = photonEnergy0 - bindingE;
230 
231  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
232  G4double pSample = profileData.RandomSelectMomentum(Z,shell);
233  // Rescale from atomic units
234  G4double pDoppler = pSample * fine_structure_const;
235  G4double pDoppler2 = pDoppler * pDoppler;
236  G4double var2 = 1. + oneCosT * e0m;
237  G4double var3 = var2*var2 - pDoppler2;
238  G4double var4 = var2 - pDoppler2 * cosTheta;
239  G4double var = var4*var4 - var3 + pDoppler2 * var3;
240  if (var > 0.)
241  {
242  G4double varSqrt = std::sqrt(var);
243  G4double scale = photonEnergy0 / var3;
244  // Random select either root
245  if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
246  else photonE = (var4 + varSqrt) * scale;
247  }
248  else
249  {
250  photonE = -1.;
251  }
252  } while ( iteration <= maxDopplerIterations &&
253  (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
254 
255  // End of recalculation of photon energy with Doppler broadening
256  // Revert to original if maximum number of iterations threshold has been reached
257  if (iteration >= maxDopplerIterations)
258  {
259  photonE = photonEoriginal;
260  bindingE = 0.;
261  }
262 
263  // Update G4VParticleChange for the scattered photon
264 
265  G4ThreeVector photonDirection1(dirX,dirY,dirZ);
266  photonDirection1.rotateUz(photonDirection0);
267  aParticleChange.ProposeMomentumDirection(photonDirection1);
268 
269  G4double photonEnergy1 = photonE;
270  //G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
271 
272  if (photonEnergy1 > 0.)
273  {
274  aParticleChange.ProposeEnergy(photonEnergy1) ;
275  }
276  else
277  {
280  }
281 
282  // Kinematics of the scattered electron
283  G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
284  G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
285 
286  G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2;
287  G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2;
288  G4double sinThetaE = -1.;
289  G4double cosThetaE = 0.;
290  if (electronP2 > 0.)
291  {
292  cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2);
293  sinThetaE = -1. * std::sqrt(1. - cosThetaE * cosThetaE);
294  }
295 
296  G4double eDirX = sinThetaE * std::cos(phi);
297  G4double eDirY = sinThetaE * std::sin(phi);
298  G4double eDirZ = cosThetaE;
299 
300  // Generate the electron only if with large enough range w.r.t. cuts and safety
301 
302  G4double safety = aStep.GetPostStepPoint()->GetSafety();
303 
304  if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
305  {
306  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
307  eDirection.rotateUz(photonDirection0);
308 
309  G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),eDirection,eKineticEnergy) ;
311  aParticleChange.AddSecondary(electron);
312  // Binding energy deposited locally
314  }
315  else
316  {
318  aParticleChange.ProposeLocalEnergyDeposit(eKineticEnergy + bindingE);
319  }
320 
321  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
322 }
323 
325 {
326  return ( &particle == G4Gamma::Gamma() );
327 }
328 
330  G4double, // previousStepSize
332 {
334  G4double energy = photon->GetKineticEnergy();
335  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
336  size_t materialIndex = couple->GetIndex();
337 
338  G4double meanFreePath;
339  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
340  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
341  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
342  return meanFreePath;
343 }