Geant4  10.03.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: G4LivermoreComptonModifiedModel.cc 95950 2016-03-03 10:42:48Z gcosmo $
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 
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4Electron.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 #include "G4Exp.hh"
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59 
60 using namespace std;
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 
65  const G4String& nam)
66  :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
67  scatterFunctionData(0),
68  crossSectionHandler(0),fAtomDeexcitation(0)
69 {
70  verboseLevel=0 ;
71  // Verbosity scale:
72  // 0 = nothing
73  // 1 = warning for energy non-conservation
74  // 2 = details of energy budget
75  // 3 = calculation of cross sections, file openings, sampling of atoms
76  // 4 = entering in methods
77 
78  if( verboseLevel>0 )
79  G4cout << "Livermore Modified Compton model is constructed " << G4endl;
80 
81  //Mark this model as "applicable" for atomic deexcitation
82  SetDeexcitationFlag(true);
83 
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 
89 {
90  delete crossSectionHandler;
91  delete scatterFunctionData;
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
95 
97  const G4DataVector& cuts)
98 {
99  if (verboseLevel > 2) {
100  G4cout << "Calling G4LivermoreComptonModifiedModel::Initialise()" << G4endl;
101  }
102 
103  if (crossSectionHandler)
104  {
105  crossSectionHandler->Clear();
106  delete crossSectionHandler;
107  }
108  delete scatterFunctionData;
109 
110  // Reading of data files - all materials are read
111  crossSectionHandler = new G4CrossSectionHandler;
112  G4String crossSectionFile = "comp/ce-cs-";
113  crossSectionHandler->LoadData(crossSectionFile);
114 
115  G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
116  G4String scatterFile = "comp/ce-sf-";
117  scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
118  scatterFunctionData->LoadData(scatterFile);
119 
120  // For Doppler broadening
121  shellData.SetOccupancyData();
122  G4String file = "/doppler/shell-doppler";
123  shellData.LoadData(file);
124 
125  InitialiseElementSelectors(particle,cuts);
126 
127  if (verboseLevel > 2) {
128  G4cout << "Loaded cross section files for Livermore Modified Compton model" << G4endl;
129  }
130 
131  if(isInitialised) { return; }
132  isInitialised = true;
133 
135 
136  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
137 
138  if( verboseLevel>0 ) {
139  G4cout << "Livermore modified Compton model is initialized " << G4endl
140  << "Energy range: "
141  << LowEnergyLimit() / eV << " eV - "
142  << HighEnergyLimit() / GeV << " GeV"
143  << G4endl;
144  }
145 }
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148 
150  const G4ParticleDefinition*,
151  G4double GammaEnergy,
154 {
155  if (verboseLevel > 3) {
156  G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModifiedModel" << G4endl;
157  }
158  if (GammaEnergy < LowEnergyLimit())
159  { return 0.0; }
160 
161  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
162  return cs;
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 
167 void G4LivermoreComptonModifiedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
168  const G4MaterialCutsCouple* couple,
169  const G4DynamicParticle* aDynamicGamma,
171 {
172 
173  // The scattered gamma energy is sampled according to Klein - Nishina formula.
174  // then accepted or rejected depending on the Scattering Function multiplied
175  // by factor from Klein - Nishina formula.
176  // Expression of the angular distribution as Klein Nishina
177  // angular and energy distribution and Scattering fuctions is taken from
178  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
179  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
180  // data are interpolated while in the article they are fitted.
181  // Reference to the article is from J. Stepanek New Photon, Positron
182  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
183  // TeV (draft).
184  // The random number techniques of Butcher & Messel are used
185  // (Nucl Phys 20(1960),15).
186 
187  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
188 
189  if (verboseLevel > 3) {
190  G4cout << "G4LivermoreComptonModifiedModel::SampleSecondaries() E(MeV)= "
191  << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
192  << G4endl;
193  }
194 
195  // do nothing below the threshold
196  // should never get here because the XS is zero below the limit
197  if (photonEnergy0 < LowEnergyLimit())
198  return ;
199 
200  G4double e0m = photonEnergy0 / electron_mass_c2 ;
201  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
202 
203  // Select randomly one element in the current material
204  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
205  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
206  G4int Z = (G4int)elm->GetZ();
207 
208  G4double epsilon0Local = 1. / (1. + 2. * e0m);
209  G4double epsilon0Sq = epsilon0Local * epsilon0Local;
210  G4double alpha1 = -std::log(epsilon0Local);
211  G4double alpha2 = 0.5 * (1. - epsilon0Sq);
212 
213  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
214 
215  // Sample the energy of the scattered photon
217  G4double epsilonSq;
218  G4double oneCosT;
219  G4double sinT2;
220  G4double gReject;
221 
222  do
223  {
224  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
225  {
226  // std::pow(epsilon0Local,G4UniformRand())
227  epsilon = G4Exp(-alpha1 * G4UniformRand());
228  epsilonSq = epsilon * epsilon;
229  }
230  else
231  {
232  epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
233  epsilon = std::sqrt(epsilonSq);
234  }
235 
236  oneCosT = (1. - epsilon) / ( epsilon * e0m);
237  sinT2 = oneCosT * (2. - oneCosT);
238  G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
239  G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
240  gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
241 
242  } while(gReject < G4UniformRand()*Z);
243 
244  G4double cosTheta = 1. - oneCosT;
245  G4double sinTheta = std::sqrt (sinT2);
246  G4double phi = twopi * G4UniformRand() ;
247  G4double dirx = sinTheta * std::cos(phi);
248  G4double diry = sinTheta * std::sin(phi);
249  G4double dirz = cosTheta ;
250 
251  // Doppler broadening - Method based on:
252  // Y. Namito, S. Ban and H. Hirayama,
253  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon
254  // into the EGS4 Code", NIM A 349, pp. 489-494, 1994
255 
256  // Maximum number of sampling iterations
257  G4int maxDopplerIterations = 1000;
258  G4double bindingE = 0.;
259  G4double photonEoriginal = epsilon * photonEnergy0;
260  G4double photonE = -1.;
261  G4int iteration = 0;
262  G4double systemE = 0;
263  G4double ePAU = -1;
264  G4int shellIdx = 0;
265  G4double vel_c = 299792458;
266  G4double momentum_au_to_nat = 1.992851740*std::pow(10.,-24.);
267  G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
268  G4double eMax = -1;
269  G4double Alpha=0;
270  do
271  {
272  ++iteration;
273  // Select shell based on shell occupancy
274  shellIdx = shellData.SelectRandomShell(Z);
275  bindingE = shellData.BindingEnergy(Z,shellIdx);
276 
277 
278 
279  // Randomly sample bound electron momentum
280  // (memento: the data set is in Atomic Units)
281  G4double pSample = profileData.RandomSelectMomentum(Z,shellIdx);
282  // Rescale from atomic units
283 
284 
285  //Kinetic energy of target electron
286 
287 
288  // Reverse vector projection onto scattering vector
289 
290  do {
291  Alpha = G4UniformRand()*pi/2.0;
292  } while(Alpha >= (pi/2.0));
293 
294  ePAU = pSample / std::cos(Alpha);
295 
296  // Convert to SI and the calculate electron energy in natural units
297 
298  G4double ePSI = ePAU * momentum_au_to_nat;
299  G4double u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
300  G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
301 
302  //Total energy of the system
303  systemE = eEIncident+photonEnergy0;
304 
305  eMax = systemE - bindingE - electron_mass_c2;
306  G4double pDoppler = pSample * fine_structure_const;
307  G4double pDoppler2 = pDoppler * pDoppler;
308  G4double var2 = 1. + oneCosT * e0m;
309  G4double var3 = var2*var2 - pDoppler2;
310  G4double var4 = var2 - pDoppler2 * cosTheta;
311  G4double var = var4*var4 - var3 + pDoppler2 * var3;
312  if (var > 0.)
313  {
314  G4double varSqrt = std::sqrt(var);
315  G4double scale = photonEnergy0 / var3;
316  // Random select either root
317  if (G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
318  else { photonE = (var4 + varSqrt) * scale; }
319  }
320  else
321  {
322  photonE = -1.;
323  }
324  } while ( iteration <= maxDopplerIterations &&
325  (photonE < 0. || photonE > eMax ) );
326 
327  // End of recalculation of photon energy with Doppler broadening
328  // Kinematics of the scattered electron
329  G4double eKineticEnergy = systemE - photonE - bindingE - electron_mass_c2;
330 
331  // protection against negative final energy: no e- is created
332  G4double eDirX = 0.0;
333  G4double eDirY = 0.0;
334  G4double eDirZ = 1.0;
335 
336  if(eKineticEnergy < 0.0) {
337  G4cout << "Error, kinetic energy of electron less than zero" << G4endl;
338  }
339 
340  else{
341  // Estimation of Compton electron polar angle taken from:
342  // The EGSnrc Code System: Monte Carlo Simulation of Electron and Photon Transport
343  // Eqn 2.2.25 Pg 42, NRCC Report PIRS-701
344  G4double E_num = photonEnergy0 - photonE*cosTheta;
345  G4double E_dom = sqrt(photonEnergy0*photonEnergy0 + photonE*photonE -2*photonEnergy0*photonE*cosTheta);
346  G4double cosThetaE = E_num / E_dom;
347  G4double sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE));
348 
349  eDirX = sinThetaE * std::cos(phi);
350  eDirY = sinThetaE * std::sin(phi);
351  eDirZ = cosThetaE;
352 
353  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
354  eDirection.rotateUz(photonDirection0);
356  eDirection,eKineticEnergy) ;
357  fvect->push_back(dp);
358  }
359 
360 
361  // Revert to original if maximum number of iterations threshold has been reached
362 
363  if (iteration >= maxDopplerIterations)
364  {
365  photonE = photonEoriginal;
366  bindingE = 0.;
367  }
368 
369  // Update G4VParticleChange for the scattered photon
370 
371  G4ThreeVector photonDirection1(dirx,diry,dirz);
372  photonDirection1.rotateUz(photonDirection0);
373  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
374 
375  G4double photonEnergy1 = photonE;
376 
377  if (photonEnergy1 > 0.)
378  {
379  fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
380 
381  if (iteration < maxDopplerIterations)
382  {
383  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
384  eDirection.rotateUz(photonDirection0);
386  eDirection,eKineticEnergy) ;
387  fvect->push_back(dp);
388  }
389  }
390  else
391  {
392  photonEnergy1 = 0.;
395  }
396 
397  // sample deexcitation
398  //
399  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
400  G4int index = couple->GetIndex();
401  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
402  size_t nbefore = fvect->size();
404  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
405  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
406  size_t nafter = fvect->size();
407  if(nafter > nbefore) {
408  for (size_t i=nbefore; i<nafter; ++i) {
409  bindingE -= ((*fvect)[i])->GetKineticEnergy();
410  }
411  }
412  }
413  }
414  if(bindingE < 0.0) { bindingE = 0.0; }
416 }
void SetOccupancyData()
Definition: G4ShellData.hh:70
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
G4double GetKineticEnergy() const
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
float h_Planck
Definition: hepunit.py:263
const G4String & GetName() const
Definition: G4Material.hh:178
G4LivermoreComptonModifiedModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermoreModifiedCompton")
G4double GetZ() const
Definition: G4Element.hh:131
G4ParticleDefinition * GetDefinition() const
tuple x
Definition: test.py:50
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int SelectRandomShell(G4int Z) const
Definition: G4ShellData.cc:363
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
void LoadData(const G4String &fileName)
Definition: G4ShellData.cc:234
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4double FindValue(G4int Z, G4double e) const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
G4double BindingEnergy(G4int Z, G4int shellIndex) const
Definition: G4ShellData.cc:166
static constexpr double cm
Definition: G4SIunits.hh:119
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
static constexpr double eV
Definition: G4SIunits.hh:215
float electron_mass_c2
Definition: hepunit.py:274
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void LoadData(const G4String &dataFile)
static constexpr double GeV
Definition: G4SIunits.hh:217
virtual G4bool LoadData(const G4String &fileName)=0
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double pi
Definition: G4SIunits.hh:75
G4VAtomDeexcitation * AtomDeexcitation()
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:788
double epsilon(double density, double temperature)
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
float c_light
Definition: hepunit.py:257
G4AtomicShellEnumerator
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
const G4Material * GetMaterial() const