Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermoreComptonModifiedModel.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 //
28 //
29 // Author: Sebastien Incerti
30 // 30 October 2008
31 // on base of G4LowEnergyCompton developed by A.Forti and M.G.Pia
32 //
33 // History:
34 // --------
35 // 18 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
36 // - apply internal high-energy limit only in constructor
37 // - do not apply low-energy limit (default is 0)
38 // - remove GetMeanFreePath method and table
39 // - added protection against numerical problem in energy sampling
40 // - use G4ElementSelector
41 // 26 Dec 2010 V Ivanchenko Load data tables only once to avoid memory leak
42 // 30 May 2011 V Ivanchenko Migration to model design for deexcitation
43 
44 #include "G4LivermoreComptonModifiedModel.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4Electron.hh"
48 #include "G4ParticleChangeForGamma.hh"
49 #include "G4LossTableManager.hh"
50 #include "G4VAtomDeexcitation.hh"
51 #include "G4AtomicShell.hh"
52 #include "G4CrossSectionHandler.hh"
53 #include "G4CompositeEMDataSet.hh"
54 #include "G4LogLogInterpolation.hh"
55 #include "G4Gamma.hh"
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58 
59 using namespace std;
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 
63 G4LivermoreComptonModifiedModel::G4LivermoreComptonModifiedModel(const G4ParticleDefinition*,
64  const G4String& nam)
65  :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
66  scatterFunctionData(0),
67  crossSectionHandler(0),fAtomDeexcitation(0)
68 {
69  lowEnergyLimit = 250 * eV;
70  highEnergyLimit = 100 * GeV;
71 
72  verboseLevel=0 ;
73  // Verbosity scale:
74  // 0 = nothing
75  // 1 = warning for energy non-conservation
76  // 2 = details of energy budget
77  // 3 = calculation of cross sections, file openings, sampling of atoms
78  // 4 = entering in methods
79 
80  if( verboseLevel>0 ) {
81  G4cout << "Livermore Modified Compton model is constructed " << G4endl
82  << "Energy range: "
83  << lowEnergyLimit / eV << " eV - "
84  << highEnergyLimit / GeV << " GeV"
85  << G4endl;
86  }
87 
88  //Mark this model as "applicable" for atomic deexcitation
89  SetDeexcitationFlag(true);
90 
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 
95 G4LivermoreComptonModifiedModel::~G4LivermoreComptonModifiedModel()
96 {
97  delete crossSectionHandler;
98  delete scatterFunctionData;
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
103 void G4LivermoreComptonModifiedModel::Initialise(const G4ParticleDefinition* particle,
104  const G4DataVector& cuts)
105 {
106  if (verboseLevel > 2) {
107  G4cout << "Calling G4LivermoreComptonModifiedModel::Initialise()" << G4endl;
108  }
109 
110  if (crossSectionHandler)
111  {
112  crossSectionHandler->Clear();
113  delete crossSectionHandler;
114  }
115  delete scatterFunctionData;
116 
117  // Reading of data files - all materials are read
118  crossSectionHandler = new G4CrossSectionHandler;
119  G4String crossSectionFile = "comp/ce-cs-";
120  crossSectionHandler->LoadData(crossSectionFile);
121 
122  G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
123  G4String scatterFile = "comp/ce-sf-";
124  scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
125  scatterFunctionData->LoadData(scatterFile);
126 
127  // For Doppler broadening
128  shellData.SetOccupancyData();
129  G4String file = "/doppler/shell-doppler";
130  shellData.LoadData(file);
131 
132  InitialiseElementSelectors(particle,cuts);
133 
134  if (verboseLevel > 2) {
135  G4cout << "Loaded cross section files for Livermore Modified Compton model" << G4endl;
136  }
137 
138  if(isInitialised) { return; }
139  isInitialised = true;
140 
141  fParticleChange = GetParticleChangeForGamma();
142 
143  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
144 
145  if( verboseLevel>0 ) {
146  G4cout << "Livermore Compton model is initialized " << G4endl
147  << "Energy range: "
148  << LowEnergyLimit() / eV << " eV - "
149  << HighEnergyLimit() / GeV << " GeV"
150  << G4endl;
151  }
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 
156 G4double G4LivermoreComptonModifiedModel::ComputeCrossSectionPerAtom(
157  const G4ParticleDefinition*,
158  G4double GammaEnergy,
159  G4double Z, G4double,
160  G4double, G4double)
161 {
162  if (verboseLevel > 3) {
163  G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModifiedModel" << G4endl;
164  }
165  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) { return 0.0; }
166 
167  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
168  return cs;
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172 
173 void G4LivermoreComptonModifiedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
174  const G4MaterialCutsCouple* couple,
175  const G4DynamicParticle* aDynamicGamma,
176  G4double, G4double)
177 {
178 
179  // The scattered gamma energy is sampled according to Klein - Nishina formula.
180  // then accepted or rejected depending on the Scattering Function multiplied
181  // by factor from Klein - Nishina formula.
182  // Expression of the angular distribution as Klein Nishina
183  // angular and energy distribution and Scattering fuctions is taken from
184  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
185  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
186  // data are interpolated while in the article they are fitted.
187  // Reference to the article is from J. Stepanek New Photon, Positron
188  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
189  // TeV (draft).
190  // The random number techniques of Butcher & Messel are used
191  // (Nucl Phys 20(1960),15).
192 
193  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
194 
195  if (verboseLevel > 3) {
196  G4cout << "G4LivermoreComptonModifiedModel::SampleSecondaries() E(MeV)= "
197  << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
198  << G4endl;
199  }
200 
201  // low-energy gamma is absorpted by this process
202  if (photonEnergy0 <= lowEnergyLimit)
203  {
204  fParticleChange->ProposeTrackStatus(fStopAndKill);
205  fParticleChange->SetProposedKineticEnergy(0.);
206  fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
207  return ;
208  }
209 
210  G4double e0m = photonEnergy0 / electron_mass_c2 ;
211  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
212 
213  // Select randomly one element in the current material
214  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
215  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
216  G4int Z = (G4int)elm->GetZ();
217 
218  G4double epsilon0Local = 1. / (1. + 2. * e0m);
219  G4double epsilon0Sq = epsilon0Local * epsilon0Local;
220  G4double alpha1 = -std::log(epsilon0Local);
221  G4double alpha2 = 0.5 * (1. - epsilon0Sq);
222 
223  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
224 
225  // Sample the energy of the scattered photon
226  G4double epsilon;
227  G4double epsilonSq;
228  G4double oneCosT;
229  G4double sinT2;
230  G4double gReject;
231 
232  do
233  {
234  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
235  {
236  // std::pow(epsilon0Local,G4UniformRand())
237  epsilon = std::exp(-alpha1 * G4UniformRand());
238  epsilonSq = epsilon * epsilon;
239  }
240  else
241  {
242  epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
243  epsilon = std::sqrt(epsilonSq);
244  }
245 
246  oneCosT = (1. - epsilon) / ( epsilon * e0m);
247  sinT2 = oneCosT * (2. - oneCosT);
248  G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
249  G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
250  gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
251 
252  } while(gReject < G4UniformRand()*Z);
253 
254  G4double cosTheta = 1. - oneCosT;
255  G4double sinTheta = std::sqrt (sinT2);
256  G4double phi = twopi * G4UniformRand() ;
257  G4double dirx = sinTheta * std::cos(phi);
258  G4double diry = sinTheta * std::sin(phi);
259  G4double dirz = cosTheta ;
260 
261  // Doppler broadening - Method based on:
262  // Y. Namito, S. Ban and H. Hirayama,
263  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon
264  // into the EGS4 Code", NIM A 349, pp. 489-494, 1994
265 
266  // Maximum number of sampling iterations
267  G4int maxDopplerIterations = 1000;
268  G4double bindingE = 0.;
269  G4double photonEoriginal = epsilon * photonEnergy0;
270  G4double photonE = -1.;
271  G4int iteration = 0;
272  G4double systemE = 0;
273  G4double ePAU = -1;
274  G4int shellIdx = 0;
275  G4double vel_c = 299792458;
276  G4double momentum_au_to_nat = 1.992851740*std::pow(10.,-24.);
277  G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
278  G4double eMax = -1;
279  G4double Alpha=0;
280  do
281  {
282  ++iteration;
283  // Select shell based on shell occupancy
284  shellIdx = shellData.SelectRandomShell(Z);
285  bindingE = shellData.BindingEnergy(Z,shellIdx);
286 
287 
288 
289  // Randomly sample bound electron momentum
290  // (memento: the data set is in Atomic Units)
291  G4double pSample = profileData.RandomSelectMomentum(Z,shellIdx);
292  // Rescale from atomic units
293 
294 
295  //Kinetic energy of target electron
296 
297 
298  // Reverse vector projection onto scattering vector
299 
300  do {
301  Alpha = G4UniformRand()*pi/2.0;
302  } while(Alpha >= (pi/2.0));
303 
304  ePAU = pSample / std::cos(Alpha);
305 
306  // Convert to SI and the calculate electron energy in natural units
307 
308  G4double ePSI = ePAU * momentum_au_to_nat;
309  G4double u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
310  G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
311 
312  //Total energy of the system
313  systemE = eEIncident+photonEnergy0;
314 
315  eMax = systemE - bindingE - electron_mass_c2;
316  G4double pDoppler = pSample * fine_structure_const;
317  G4double pDoppler2 = pDoppler * pDoppler;
318  G4double var2 = 1. + oneCosT * e0m;
319  G4double var3 = var2*var2 - pDoppler2;
320  G4double var4 = var2 - pDoppler2 * cosTheta;
321  G4double var = var4*var4 - var3 + pDoppler2 * var3;
322  if (var > 0.)
323  {
324  G4double varSqrt = std::sqrt(var);
325  G4double scale = photonEnergy0 / var3;
326  // Random select either root
327  if (G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
328  else { photonE = (var4 + varSqrt) * scale; }
329  }
330  else
331  {
332  photonE = -1.;
333  }
334  } while ( iteration <= maxDopplerIterations &&
335  (photonE < 0. || photonE > eMax ) );
336 
337  // End of recalculation of photon energy with Doppler broadening
338  // Kinematics of the scattered electron
339  G4double eKineticEnergy = systemE - photonE - bindingE - electron_mass_c2;
340 
341  // protection against negative final energy: no e- is created
342  G4double eDirX = 0.0;
343  G4double eDirY = 0.0;
344  G4double eDirZ = 1.0;
345 
346  if(eKineticEnergy < 0.0) {
347  G4cout << "Error, kinetic energy of electron less than zero" << G4endl;
348  }
349 
350  else{
351  // Estimation of Compton electron polar angle taken from:
352  // The EGSnrc Code System: Monte Carlo Simulation of Electron and Photon Transport
353  // Eqn 2.2.25 Pg 42, NRCC Report PIRS-701
354  G4double E_num = photonEnergy0 - photonE*cosTheta;
355  G4double E_dom = sqrt(photonEnergy0*photonEnergy0 + photonE*photonE -2*photonEnergy0*photonE*cosTheta);
356  G4double cosThetaE = E_num / E_dom;
357  G4double sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE));
358 
359  eDirX = sinThetaE * std::cos(phi);
360  eDirY = sinThetaE * std::sin(phi);
361  eDirZ = cosThetaE;
362 
363  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
364  eDirection.rotateUz(photonDirection0);
365  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
366  eDirection,eKineticEnergy) ;
367  fvect->push_back(dp);
368  }
369 
370 
371  // Revert to original if maximum number of iterations threshold has been reached
372 
373  if (iteration >= maxDopplerIterations)
374  {
375  photonE = photonEoriginal;
376  bindingE = 0.;
377  }
378 
379  // Update G4VParticleChange for the scattered photon
380 
381  G4ThreeVector photonDirection1(dirx,diry,dirz);
382  photonDirection1.rotateUz(photonDirection0);
383  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
384 
385  G4double photonEnergy1 = photonE;
386 
387  if (photonEnergy1 > 0.)
388  {
389  fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
390 
391  if (iteration < maxDopplerIterations)
392  {
393  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
394  eDirection.rotateUz(photonDirection0);
395  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
396  eDirection,eKineticEnergy) ;
397  fvect->push_back(dp);
398  }
399  }
400  else
401  {
402  photonEnergy1 = 0.;
403  fParticleChange->SetProposedKineticEnergy(0.) ;
404  fParticleChange->ProposeTrackStatus(fStopAndKill);
405  }
406 
407  // sample deexcitation
408  //
409  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
410  G4int index = couple->GetIndex();
411  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
412  size_t nbefore = fvect->size();
413  G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
414  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
415  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
416  size_t nafter = fvect->size();
417  if(nafter > nbefore) {
418  for (size_t i=nbefore; i<nafter; ++i) {
419  bindingE -= ((*fvect)[i])->GetKineticEnergy();
420  }
421  }
422  }
423  }
424  if(bindingE < 0.0) { bindingE = 0.0; }
425  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
426 }
427 
Generated on Sat May 25 2013 14:33:34 for Geant4 by   doxygen 1.8.4