source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreComptonModel.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 12.6 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: G4LivermoreComptonModel.cc,v 1.7 2009/06/10 13:32:36 mantero Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29//
30// Author: Sebastien Inserti
31//         30 October 2008
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
42#include "G4LivermoreComptonModel.hh"
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45
46using namespace std;
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50G4LivermoreComptonModel::G4LivermoreComptonModel(const G4ParticleDefinition*,
51                                                 const G4String& nam)
52  :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
53   scatterFunctionData(0),crossSectionHandler(0)
54{
55  lowEnergyLimit = 250 * eV; 
56  highEnergyLimit = 100 * GeV;
57  //  SetLowEnergyLimit(lowEnergyLimit);
58  SetHighEnergyLimit(highEnergyLimit);
59
60  verboseLevel=0 ;
61  // Verbosity scale:
62  // 0 = nothing
63  // 1 = warning for energy non-conservation
64  // 2 = details of energy budget
65  // 3 = calculation of cross sections, file openings, sampling of atoms
66  // 4 = entering in methods
67
68  if(  verboseLevel>0 ) { 
69    G4cout << "Livermore Compton model is constructed " << G4endl
70           << "Energy range: "
71           << lowEnergyLimit / eV << " eV - "
72           << highEnergyLimit / GeV << " GeV"
73           << G4endl;
74  }
75}
76
77//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78
79G4LivermoreComptonModel::~G4LivermoreComptonModel()
80{ 
81  if (crossSectionHandler) delete crossSectionHandler;
82  if (scatterFunctionData) delete scatterFunctionData;
83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87void G4LivermoreComptonModel::Initialise(const G4ParticleDefinition* particle,
88                                         const G4DataVector& cuts)
89{
90  if (verboseLevel > 3)
91    G4cout << "Calling G4LivermoreComptonModel::Initialise()" << G4endl;
92
93  if (crossSectionHandler)
94  {
95    crossSectionHandler->Clear();
96    delete crossSectionHandler;
97  }
98 
99  // Reading of data files - all materials are read
100 
101  crossSectionHandler = new G4CrossSectionHandler;
102  crossSectionHandler->Clear();
103  G4String crossSectionFile = "comp/ce-cs-";
104  crossSectionHandler->LoadData(crossSectionFile);
105
106  G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
107  G4String scatterFile = "comp/ce-sf-";
108  scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
109  scatterFunctionData->LoadData(scatterFile);
110
111  // For Doppler broadening
112  shellData.SetOccupancyData();
113  G4String file = "/doppler/shell-doppler";
114  shellData.LoadData(file);
115
116  if (verboseLevel > 2) 
117    G4cout << "Loaded cross section files for Livermore Compton model" << G4endl;
118
119  InitialiseElementSelectors(particle,cuts);
120
121  if(  verboseLevel>0 ) { 
122    G4cout << "Livermore Compton model is initialized " << G4endl
123           << "Energy range: "
124           << LowEnergyLimit() / eV << " eV - "
125           << HighEnergyLimit() / GeV << " GeV"
126           << G4endl;
127  }
128  // 
129  if(isInitialised) return;
130  fParticleChange = GetParticleChangeForGamma();
131  isInitialised = true;
132}
133
134//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135
136G4double G4LivermoreComptonModel::ComputeCrossSectionPerAtom(
137                                       const G4ParticleDefinition*,
138                                             G4double GammaEnergy,
139                                             G4double Z, G4double,
140                                             G4double, G4double)
141{
142  if (verboseLevel > 3)
143    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModel" << G4endl;
144
145  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
146   
147  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy); 
148  return cs;
149}
150
151//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152
153void G4LivermoreComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
154                                                const G4MaterialCutsCouple* couple,
155                                                const G4DynamicParticle* aDynamicGamma,
156                                                G4double, G4double)
157{
158
159  // The scattered gamma energy is sampled according to Klein - Nishina formula.
160  // then accepted or rejected depending on the Scattering Function multiplied
161  // by factor from Klein - Nishina formula.
162  // Expression of the angular distribution as Klein Nishina
163  // angular and energy distribution and Scattering fuctions is taken from
164  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
165  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
166  // data are interpolated while in the article they are fitted.
167  // Reference to the article is from J. Stepanek New Photon, Positron
168  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
169  // TeV (draft).
170  // The random number techniques of Butcher & Messel are used
171  // (Nucl Phys 20(1960),15).
172
173  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
174
175  if (verboseLevel > 3) {
176    G4cout << "G4LivermoreComptonModel::SampleSecondaries() E(MeV)= " 
177           << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName() 
178           << G4endl;
179  }
180 
181  // low-energy gamma is absorpted by this process
182  if (photonEnergy0 <= lowEnergyLimit) 
183    {
184      fParticleChange->ProposeTrackStatus(fStopAndKill);
185      fParticleChange->SetProposedKineticEnergy(0.);
186      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
187      return ;
188    }
189
190  G4double e0m = photonEnergy0 / electron_mass_c2 ;
191  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
192
193  // Select randomly one element in the current material
194  //  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
195  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
196  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
197  G4int Z = (G4int)elm->GetZ();
198
199  G4double epsilon0 = 1. / (1. + 2. * e0m);
200  G4double epsilon0Sq = epsilon0 * epsilon0;
201  G4double alpha1 = -std::log(epsilon0);
202  G4double alpha2 = 0.5 * (1. - epsilon0Sq);
203
204  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
205
206  // Sample the energy of the scattered photon
207  G4double epsilon;
208  G4double epsilonSq;
209  G4double oneCosT;
210  G4double sinT2;
211  G4double gReject;
212 
213  do
214  {
215      if ( alpha1/(alpha1+alpha2) > G4UniformRand())
216      {
217        // std::pow(epsilon0,G4UniformRand())
218        epsilon = std::exp(-alpha1 * G4UniformRand()); 
219        epsilonSq = epsilon * epsilon;
220      }
221      else
222      {
223        epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
224        epsilon = std::sqrt(epsilonSq);
225      }
226
227      oneCosT = (1. - epsilon) / ( epsilon * e0m);
228      sinT2 = oneCosT * (2. - oneCosT);
229      G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
230      G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
231      gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
232
233  } while(gReject < G4UniformRand()*Z);
234
235  G4double cosTheta = 1. - oneCosT;
236  G4double sinTheta = std::sqrt (sinT2);
237  G4double phi = twopi * G4UniformRand() ;
238  G4double dirx = sinTheta * std::cos(phi);
239  G4double diry = sinTheta * std::sin(phi);
240  G4double dirz = cosTheta ;
241
242  // Doppler broadening -  Method based on:
243  // Y. Namito, S. Ban and H. Hirayama,
244  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon
245  // into the EGS4 Code", NIM A 349, pp. 489-494, 1994
246 
247  // Maximum number of sampling iterations
248  G4int maxDopplerIterations = 1000;
249  G4double bindingE = 0.;
250  G4double photonEoriginal = epsilon * photonEnergy0;
251  G4double photonE = -1.;
252  G4int iteration = 0;
253  G4double eMax = photonEnergy0;
254  do
255    {
256      iteration++;
257      // Select shell based on shell occupancy
258      G4int shell = shellData.SelectRandomShell(Z);
259      bindingE = shellData.BindingEnergy(Z,shell);
260     
261      eMax = photonEnergy0 - bindingE;
262     
263      // Randomly sample bound electron momentum
264      // (memento: the data set is in Atomic Units)
265      G4double pSample = profileData.RandomSelectMomentum(Z,shell);
266      // Rescale from atomic units
267      G4double pDoppler = pSample * fine_structure_const;
268      G4double pDoppler2 = pDoppler * pDoppler;
269      G4double var2 = 1. + oneCosT * e0m;
270      G4double var3 = var2*var2 - pDoppler2;
271      G4double var4 = var2 - pDoppler2 * cosTheta;
272      G4double var = var4*var4 - var3 + pDoppler2 * var3;
273      if (var > 0.)
274        {
275          G4double varSqrt = std::sqrt(var);       
276          G4double scale = photonEnergy0 / var3; 
277          // Random select either root
278          if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;               
279          else photonE = (var4 + varSqrt) * scale;
280        } 
281      else
282        {
283          photonE = -1.;
284        }
285   } while ( iteration <= maxDopplerIterations && 
286             (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
287 
288  // End of recalculation of photon energy with Doppler broadening
289  // Revert to original if maximum number of iterations threshold has been reached
290
291  if (iteration >= maxDopplerIterations)
292    {
293      photonE = photonEoriginal;
294      bindingE = 0.;
295    }
296
297  // Update G4VParticleChange for the scattered photon
298
299  G4ThreeVector photonDirection1(dirx,diry,dirz);
300  photonDirection1.rotateUz(photonDirection0);
301  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
302
303  G4double photonEnergy1 = photonE;
304
305  if (photonEnergy1 > 0.)
306    {
307      fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
308    }
309  else
310    {
311      photonEnergy1 = 0.;
312      fParticleChange->SetProposedKineticEnergy(0.) ;
313      fParticleChange->ProposeTrackStatus(fStopAndKill);   
314    }
315
316  // Kinematics of the scattered electron
317  G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
318
319  // protection against negative final energy: no e- is created
320  if(eKineticEnergy < 0.0) {
321    fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0 - photonEnergy1);
322    return;
323  }
324  G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
325
326  G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2; 
327  G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2;
328  G4double sinThetaE = -1.;
329  G4double cosThetaE = 0.;
330  if (electronP2 > 0.)
331    {
332      cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2);
333      sinThetaE = -1. * sqrt(1. - cosThetaE * cosThetaE); 
334    }
335 
336  G4double eDirX = sinThetaE * std::cos(phi);
337  G4double eDirY = sinThetaE * std::sin(phi);
338  G4double eDirZ = cosThetaE;
339
340  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
341  eDirection.rotateUz(photonDirection0);
342
343  // SI - The range test has been removed wrt original G4LowEnergyCompton class
344
345  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
346 
347  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
348                                                 eDirection,eKineticEnergy) ;
349  fvect->push_back(dp);
350}
351
Note: See TracBrowser for help on using the repository browser.