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