source: trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyCompton.cc@ 1285

Last change on this file since 1285 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 13.0 KB
Line 
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: G4LowEnergyCompton.cc,v 1.50 2009/06/11 15:47:08 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03 $
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 "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"
59#include "G4VCrossSectionHandler.hh"
60#include "G4CrossSectionHandler.hh"
61#include "G4VEMDataSet.hh"
62#include "G4CompositeEMDataSet.hh"
63#include "G4VDataSetAlgorithm.hh"
64#include "G4LogLogInterpolation.hh"
65#include "G4VRangeTest.hh"
66#include "G4RangeTest.hh"
67#include "G4RangeNoTest.hh"
68#include "G4MaterialCutsCouple.hh"
69
70
71G4LowEnergyCompton::G4LowEnergyCompton(const G4String& processName)
72 : G4VDiscreteProcess(processName),
73 lowEnergyLimit(250*eV),
74 highEnergyLimit(100*GeV),
75 intrinsicLowEnergyLimit(10*eV),
76 intrinsicHighEnergyLimit(100*GeV)
77{
78 if (lowEnergyLimit < intrinsicLowEnergyLimit ||
79 highEnergyLimit > intrinsicHighEnergyLimit)
80 {
81 G4Exception("G4LowEnergyCompton::G4LowEnergyCompton - energy outside intrinsic process validity range");
82 }
83
84 crossSectionHandler = new G4CrossSectionHandler;
85
86 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
87 G4String scatterFile = "comp/ce-sf-";
88 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
89 scatterFunctionData->LoadData(scatterFile);
90
91 meanFreePathTable = 0;
92
93 rangeTest = new G4RangeNoTest;
94
95 // For Doppler broadening
96 shellData.SetOccupancyData();
97
98 if (verboseLevel > 0)
99 {
100 G4cout << GetProcessName() << " is created " << G4endl
101 << "Energy range: "
102 << lowEnergyLimit / keV << " keV - "
103 << highEnergyLimit / GeV << " GeV"
104 << G4endl;
105 }
106
107 G4cout << G4endl;
108 G4cout << "*******************************************************************************" << G4endl;
109 G4cout << "*******************************************************************************" << G4endl;
110 G4cout << " The class G4LowEnergyCompton is NOT SUPPORTED ANYMORE. " << G4endl;
111 G4cout << " It will be REMOVED with the next major release of Geant4. " << G4endl;
112 G4cout << " Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
113 G4cout << "*******************************************************************************" << G4endl;
114 G4cout << "*******************************************************************************" << G4endl;
115 G4cout << G4endl;
116}
117
118G4LowEnergyCompton::~G4LowEnergyCompton()
119{
120 delete meanFreePathTable;
121 delete crossSectionHandler;
122 delete scatterFunctionData;
123 delete rangeTest;
124}
125
126void G4LowEnergyCompton::BuildPhysicsTable(const G4ParticleDefinition& )
127{
128
129 crossSectionHandler->Clear();
130 G4String crossSectionFile = "comp/ce-cs-";
131 crossSectionHandler->LoadData(crossSectionFile);
132
133 delete meanFreePathTable;
134 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
135
136 // For Doppler broadening
137 G4String file = "/doppler/shell-doppler";
138 shellData.LoadData(file);
139}
140
141G4VParticleChange* G4LowEnergyCompton::PostStepDoIt(const G4Track& aTrack,
142 const G4Step& aStep)
143{
144 // The scattered gamma energy is sampled according to Klein - Nishina formula.
145 // then accepted or rejected depending on the Scattering Function multiplied
146 // by factor from Klein - Nishina formula.
147 // Expression of the angular distribution as Klein Nishina
148 // angular and energy distribution and Scattering fuctions is taken from
149 // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
150 // Phys. Res. B 101 (1995). Method of sampling with form factors is different
151 // data are interpolated while in the article they are fitted.
152 // Reference to the article is from J. Stepanek New Photon, Positron
153 // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
154 // TeV (draft).
155 // The random number techniques of Butcher & Messel are used
156 // (Nucl Phys 20(1960),15).
157
158 aParticleChange.Initialize(aTrack);
159
160 // Dynamic particle quantities
161 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
162 G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
163
164 if (photonEnergy0 <= lowEnergyLimit)
165 {
166 aParticleChange.ProposeTrackStatus(fStopAndKill);
167 aParticleChange.ProposeEnergy(0.);
168 aParticleChange.ProposeLocalEnergyDeposit(photonEnergy0);
169 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
170 }
171
172 G4double e0m = photonEnergy0 / electron_mass_c2 ;
173 G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
174 G4double epsilon0 = 1. / (1. + 2. * e0m);
175 G4double epsilon0Sq = epsilon0 * epsilon0;
176 G4double alpha1 = -std::log(epsilon0);
177 G4double alpha2 = 0.5 * (1. - epsilon0Sq);
178 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
179
180 // Select randomly one element in the current material
181 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
182 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
183
184 // Sample the energy of the scattered photon
185 G4double epsilon;
186 G4double epsilonSq;
187 G4double oneCosT;
188 G4double sinT2;
189 G4double gReject;
190 do
191 {
192 if ( alpha1/(alpha1+alpha2) > G4UniformRand())
193 {
194 epsilon = std::exp(-alpha1 * G4UniformRand()); // std::pow(epsilon0,G4UniformRand())
195 epsilonSq = epsilon * epsilon;
196 }
197 else
198 {
199 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
200 epsilon = std::sqrt(epsilonSq);
201 }
202
203 oneCosT = (1. - epsilon) / ( epsilon * e0m);
204 sinT2 = oneCosT * (2. - oneCosT);
205 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
206 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
207 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
208
209 } while(gReject < G4UniformRand()*Z);
210
211 G4double cosTheta = 1. - oneCosT;
212 G4double sinTheta = std::sqrt (sinT2);
213 G4double phi = twopi * G4UniformRand() ;
214 G4double dirX = sinTheta * std::cos(phi);
215 G4double dirY = sinTheta * std::sin(phi);
216 G4double dirZ = cosTheta ;
217
218 // Doppler broadening - Method based on:
219 // Y. Namito, S. Ban and H. Hirayama,
220 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
221 // NIM A 349, pp. 489-494, 1994
222
223 // Maximum number of sampling iterations
224 G4int maxDopplerIterations = 1000;
225 G4double bindingE = 0.;
226 G4double photonEoriginal = epsilon * photonEnergy0;
227 G4double photonE = -1.;
228 G4int iteration = 0;
229 G4double eMax = photonEnergy0;
230 do
231 {
232 iteration++;
233 // Select shell based on shell occupancy
234 G4int shell = shellData.SelectRandomShell(Z);
235 bindingE = shellData.BindingEnergy(Z,shell);
236
237 eMax = photonEnergy0 - bindingE;
238
239 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
240 G4double pSample = profileData.RandomSelectMomentum(Z,shell);
241 // Rescale from atomic units
242 G4double pDoppler = pSample * fine_structure_const;
243 G4double pDoppler2 = pDoppler * pDoppler;
244 G4double var2 = 1. + oneCosT * e0m;
245 G4double var3 = var2*var2 - pDoppler2;
246 G4double var4 = var2 - pDoppler2 * cosTheta;
247 G4double var = var4*var4 - var3 + pDoppler2 * var3;
248 if (var > 0.)
249 {
250 G4double varSqrt = std::sqrt(var);
251 G4double scale = photonEnergy0 / var3;
252 // Random select either root
253 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
254 else photonE = (var4 + varSqrt) * scale;
255 }
256 else
257 {
258 photonE = -1.;
259 }
260 } while ( iteration <= maxDopplerIterations &&
261 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
262
263 // End of recalculation of photon energy with Doppler broadening
264 // Revert to original if maximum number of iterations threshold has been reached
265 if (iteration >= maxDopplerIterations)
266 {
267 photonE = photonEoriginal;
268 bindingE = 0.;
269 }
270
271 // Update G4VParticleChange for the scattered photon
272
273 G4ThreeVector photonDirection1(dirX,dirY,dirZ);
274 photonDirection1.rotateUz(photonDirection0);
275 aParticleChange.ProposeMomentumDirection(photonDirection1);
276
277 G4double photonEnergy1 = photonE;
278 //G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
279
280 if (photonEnergy1 > 0.)
281 {
282 aParticleChange.ProposeEnergy(photonEnergy1) ;
283 }
284 else
285 {
286 aParticleChange.ProposeEnergy(0.) ;
287 aParticleChange.ProposeTrackStatus(fStopAndKill);
288 }
289
290 // Kinematics of the scattered electron
291 G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
292 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
293
294 G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2;
295 G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2;
296 G4double sinThetaE = -1.;
297 G4double cosThetaE = 0.;
298 if (electronP2 > 0.)
299 {
300 cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2);
301 sinThetaE = -1. * std::sqrt(1. - cosThetaE * cosThetaE);
302 }
303
304 G4double eDirX = sinThetaE * std::cos(phi);
305 G4double eDirY = sinThetaE * std::sin(phi);
306 G4double eDirZ = cosThetaE;
307
308 // Generate the electron only if with large enough range w.r.t. cuts and safety
309
310 G4double safety = aStep.GetPostStepPoint()->GetSafety();
311
312 if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
313 {
314 G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
315 eDirection.rotateUz(photonDirection0);
316
317 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),eDirection,eKineticEnergy) ;
318 aParticleChange.SetNumberOfSecondaries(1);
319 aParticleChange.AddSecondary(electron);
320 // Binding energy deposited locally
321 aParticleChange.ProposeLocalEnergyDeposit(bindingE);
322 }
323 else
324 {
325 aParticleChange.SetNumberOfSecondaries(0);
326 aParticleChange.ProposeLocalEnergyDeposit(eKineticEnergy + bindingE);
327 }
328
329 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
330}
331
332G4bool G4LowEnergyCompton::IsApplicable(const G4ParticleDefinition& particle)
333{
334 return ( &particle == G4Gamma::Gamma() );
335}
336
337G4double G4LowEnergyCompton::GetMeanFreePath(const G4Track& track,
338 G4double, // previousStepSize
339 G4ForceCondition*)
340{
341 const G4DynamicParticle* photon = track.GetDynamicParticle();
342 G4double energy = photon->GetKineticEnergy();
343 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
344 size_t materialIndex = couple->GetIndex();
345
346 G4double meanFreePath;
347 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
348 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
349 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
350 return meanFreePath;
351}
Note: See TracBrowser for help on using the repository browser.