ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KleinNishinaCompton.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4KleinNishinaCompton.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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4KleinNishinaCompton
33 //
34 // Author: Vladimir Ivanchenko on base of Michel Maire code
35 //
36 // Creation date: 15.03.2005
37 //
38 // Modifications:
39 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivantchenko)
40 // 27-03-06 Remove upper limit of cross section (V.Ivantchenko)
41 //
42 // Class Description:
43 //
44 // -------------------------------------------------------------------
45 //
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 #include "G4KleinNishinaCompton.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4Electron.hh"
53 #include "G4Gamma.hh"
54 #include "Randomize.hh"
55 #include "G4DataVector.hh"
56 #include "G4ParticleChangeForGamma.hh"
57 #include "G4Log.hh"
58 #include "G4Exp.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 
62 using namespace std;
63 
64 G4KleinNishinaCompton::G4KleinNishinaCompton(const G4ParticleDefinition*,
65  const G4String& nam)
66  : G4VEmModel(nam)
67 {
68  theGamma = G4Gamma::Gamma();
69  theElectron = G4Electron::Electron();
70  lowestSecondaryEnergy = 100.0*eV;
71  fParticleChange = nullptr;
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
76 G4KleinNishinaCompton::~G4KleinNishinaCompton()
77 {}
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
81 void G4KleinNishinaCompton::Initialise(const G4ParticleDefinition* p,
82  const G4DataVector& cuts)
83 {
84  if(IsMaster()) { InitialiseElementSelectors(p, cuts); }
85  if(nullptr == fParticleChange) {
86  fParticleChange = GetParticleChangeForGamma();
87  }
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
92 void G4KleinNishinaCompton::InitialiseLocal(const G4ParticleDefinition*,
93  G4VEmModel* masterModel)
94 {
95  SetElementSelectors(masterModel->GetElementSelectors());
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99 
100 G4double G4KleinNishinaCompton::ComputeCrossSectionPerAtom(
101  const G4ParticleDefinition*,
102  G4double GammaEnergy,
103  G4double Z, G4double,
104  G4double, G4double)
105 {
106  G4double xSection = 0.0 ;
107  if (GammaEnergy <= LowEnergyLimit()) { return xSection; }
108 
109  static const G4double a = 20.0 , b = 230.0 , c = 440.0;
110 
111  static const G4double
112  d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLHEP::barn,
113  d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLHEP::barn,
114  e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLHEP::barn,
115  e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLHEP::barn,
116  f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLHEP::barn,
117  f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLHEP::barn;
118 
119  G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
120  p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
121 
122  G4double T0 = 15.0*keV;
123  if (Z < 1.5) { T0 = 40.0*keV; }
124 
125  G4double X = max(GammaEnergy, T0) / electron_mass_c2;
126  xSection = p1Z*G4Log(1.+2.*X)/X
127  + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
128 
129  // modification for low energy. (special case for Hydrogen)
130  if (GammaEnergy < T0) {
131  static const G4double dT0 = keV;
132  X = (T0+dT0) / electron_mass_c2 ;
133  G4double sigma = p1Z*G4Log(1.+2*X)/X
134  + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
135  G4double c1 = -T0*(sigma-xSection)/(xSection*dT0);
136  G4double c2 = 0.150;
137  if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); }
138  G4double y = G4Log(GammaEnergy/T0);
139  xSection *= G4Exp(-y*(c1+c2*y));
140  }
141  // G4cout<<"e= "<< GammaEnergy<<" Z= "<<Z<<" cross= " << xSection << G4endl;
142  return xSection;
143 }
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 
147 void G4KleinNishinaCompton::SampleSecondaries(
148  std::vector<G4DynamicParticle*>* fvect,
149  const G4MaterialCutsCouple*,
150  const G4DynamicParticle* aDynamicGamma,
151  G4double,
152  G4double)
153 {
154  // The scattered gamma energy is sampled according to Klein - Nishina formula.
155  // The random number techniques of Butcher & Messel are used
156  // (Nuc Phys 20(1960),15).
157  // Note : Effects due to binding of atomic electrons are negliged.
158 
159  G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
160 
161  // do nothing below the threshold
162  if(gamEnergy0 <= LowEnergyLimit()) { return; }
163 
164  G4double E0_m = gamEnergy0 / electron_mass_c2 ;
165 
166  G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
167 
168  //
169  // sample the energy rate of the scattered gamma
170  //
171 
172  G4double epsilon, epsilonsq, onecost, sint2, greject ;
173 
174  G4double eps0 = 1./(1. + 2.*E0_m);
175  G4double epsilon0sq = eps0*eps0;
176  G4double alpha1 = - G4Log(eps0);
177  G4double alpha2 = alpha1 + 0.5*(1.- epsilon0sq);
178 
179  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
180  G4double rndm[3];
181 
182  static const G4int nlooplim = 1000;
183  G4int nloop = 0;
184  do {
185  ++nloop;
186  // false interaction if too many iterations
187  if(nloop > nlooplim) { return; }
188 
189  // 3 random numbers to sample scattering
190  rndmEngineMod->flatArray(3, rndm);
191 
192  if ( alpha1 > alpha2*rndm[0] ) {
193  epsilon = G4Exp(-alpha1*rndm[1]); // eps0**r
194  epsilonsq = epsilon*epsilon;
195 
196  } else {
197  epsilonsq = epsilon0sq + (1.- epsilon0sq)*rndm[1];
198  epsilon = sqrt(epsilonsq);
199  };
200 
201  onecost = (1.- epsilon)/(epsilon*E0_m);
202  sint2 = onecost*(2.-onecost);
203  greject = 1. - epsilon*sint2/(1.+ epsilonsq);
204 
205  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
206  } while (greject < rndm[2]);
207 
208  //
209  // scattered gamma angles. ( Z - axis along the parent gamma)
210  //
211 
212  if(sint2 < 0.0) { sint2 = 0.0; }
213  G4double cosTeta = 1. - onecost;
214  G4double sinTeta = sqrt (sint2);
215  G4double Phi = twopi * rndmEngineMod->flat();
216 
217  //
218  // update G4VParticleChange for the scattered gamma
219  //
220 
221  G4ThreeVector gamDirection1(sinTeta*cos(Phi), sinTeta*sin(Phi), cosTeta);
222  gamDirection1.rotateUz(gamDirection0);
223  G4double gamEnergy1 = epsilon*gamEnergy0;
224  G4double edep = 0.0;
225  if(gamEnergy1 > lowestSecondaryEnergy) {
226  fParticleChange->ProposeMomentumDirection(gamDirection1);
227  fParticleChange->SetProposedKineticEnergy(gamEnergy1);
228  } else {
229  fParticleChange->ProposeTrackStatus(fStopAndKill);
230  fParticleChange->SetProposedKineticEnergy(0.0);
231  edep = gamEnergy1;
232  }
233 
234  //
235  // kinematic of the scattered electron
236  //
237 
238  G4double eKinEnergy = gamEnergy0 - gamEnergy1;
239 
240  if(eKinEnergy > lowestSecondaryEnergy) {
241  G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
242  eDirection = eDirection.unit();
243 
244  // create G4DynamicParticle object for the electron.
245  G4DynamicParticle* dp = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
246  fvect->push_back(dp);
247  } else {
248  edep += eKinEnergy;
249  }
250  // energy balance
251  if(edep > 0.0) {
252  fParticleChange->ProposeLocalEnergyDeposit(edep);
253  }
254 }
255 
256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
257 
258