ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KleinNishinaModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4KleinNishinaModel.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: G4KleinNishinaModel
33 //
34 // Author: Vladimir Ivanchenko on base of G4KleinNishinaCompton
35 //
36 // Creation date: 13.06.2010
37 //
38 // Modifications:
39 //
40 // Class Description:
41 //
42 // -------------------------------------------------------------------
43 //
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 
47 #include "G4KleinNishinaModel.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4Electron.hh"
51 #include "G4Gamma.hh"
52 #include "Randomize.hh"
53 #include "G4RandomDirection.hh"
54 #include "G4DataVector.hh"
55 #include "G4ParticleChangeForGamma.hh"
56 #include "G4VAtomDeexcitation.hh"
57 #include "G4AtomicShells.hh"
58 #include "G4LossTableManager.hh"
59 #include "G4Log.hh"
60 #include "G4Exp.hh"
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 
64 using namespace std;
65 
66 G4KleinNishinaModel::G4KleinNishinaModel(const G4String& nam)
67  : G4VEmModel(nam),
68  lv1(0.,0.,0.,0.),
69  lv2(0.,0.,0.,0.),
70  bst(0.,0.,0.)
71 {
72  theGamma = G4Gamma::Gamma();
73  theElectron = G4Electron::Electron();
74  lowestSecondaryEnergy = 10*eV;
75  limitFactor = 4;
76  fProbabilities.resize(9,0.0);
77  SetDeexcitationFlag(true);
78  fParticleChange = nullptr;
79  fAtomDeexcitation = nullptr;
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
84 G4KleinNishinaModel::~G4KleinNishinaModel()
85 {}
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
89 void G4KleinNishinaModel::Initialise(const G4ParticleDefinition* p,
90  const G4DataVector& cuts)
91 {
92  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
93  if(IsMaster()) { InitialiseElementSelectors(p, cuts); }
94  if(nullptr == fParticleChange) {
95  fParticleChange = GetParticleChangeForGamma();
96  }
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
101 void G4KleinNishinaModel::InitialiseLocal(const G4ParticleDefinition*,
102  G4VEmModel* masterModel)
103 {
104  SetElementSelectors(masterModel->GetElementSelectors());
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 
109 G4double
110 G4KleinNishinaModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
111  G4double gammaEnergy,
112  G4double Z, G4double,
113  G4double, G4double)
114 {
115  G4double xSection = 0.0 ;
116  if (gammaEnergy <= LowEnergyLimit()) { return xSection; }
117 
118  static const G4double a = 20.0 , b = 230.0 , c = 440.0;
119 
120 static const G4double
121  d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLHEP::barn,
122  d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLHEP::barn,
123  e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLHEP::barn,
124  e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLHEP::barn,
125  f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLHEP::barn,
126  f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLHEP::barn;
127 
128  G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
129  p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
130 
131  G4double T0 = 15.0*keV;
132  if (Z < 1.5) { T0 = 40.0*keV; }
133 
134  G4double X = max(gammaEnergy, T0) / electron_mass_c2;
135  xSection = p1Z*G4Log(1.+2.*X)/X
136  + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
137 
138  // modification for low energy. (special case for Hydrogen)
139  static const G4double dT0 = keV;
140  if (gammaEnergy < T0) {
141  X = (T0+dT0) / electron_mass_c2 ;
142  G4double sigma = p1Z*G4Log(1.+2*X)/X
143  + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
144  G4double c1 = -T0*(sigma-xSection)/(xSection*dT0);
145  G4double c2 = 0.150;
146  if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); }
147  G4double y = G4Log(gammaEnergy/T0);
148  xSection *= G4Exp(-y*(c1+c2*y));
149  }
150 
151  if(xSection < 0.0) { xSection = 0.0; }
152  // G4cout << "e= " << GammaEnergy << " Z= " << Z
153  // << " cross= " << xSection << G4endl;
154  return xSection;
155 }
156 
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158 
159 void G4KleinNishinaModel::SampleSecondaries(
160  std::vector<G4DynamicParticle*>* fvect,
161  const G4MaterialCutsCouple* couple,
162  const G4DynamicParticle* aDynamicGamma,
163  G4double,
164  G4double)
165 {
166  // primary gamma
167  G4double energy = aDynamicGamma->GetKineticEnergy();
168 
169  // do nothing below the threshold
170  if(energy <= LowEnergyLimit()) { return; }
171 
172  G4ThreeVector direction = aDynamicGamma->GetMomentumDirection();
173 
174  // select atom
175  const G4Element* elm = SelectRandomAtom(couple, theGamma, energy);
176 
177  // select shell first
178  G4int nShells = elm->GetNbOfAtomicShells();
179  if(nShells > (G4int)fProbabilities.size()) { fProbabilities.resize(nShells); }
180  G4double totprob = 0.0;
181  G4int i;
182  for(i=0; i<nShells; ++i) {
183  //G4double bindingEnergy = elm->GetAtomicShell(i);
184  totprob += elm->GetNbOfShellElectrons(i);
185  //totprob += elm->GetNbOfShellElectrons(i)/(bindingEnergy*bindingEnergy);
186  fProbabilities[i] = totprob;
187  }
188 
189  // Loop on sampling
190  static const G4int nlooplim = 1000;
191  G4int nloop = 0;
192 
193  G4double bindingEnergy, ePotEnergy, eKinEnergy;
194  G4double gamEnergy0, gamEnergy1;
195 
196  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
197  G4double rndm[4];
198 
199  do {
200  ++nloop;
201 
202  // 4 random numbers to select e-
203  rndmEngineMod->flatArray(4, rndm);
204  G4double xprob = totprob*rndm[0];
205 
206  // select shell
207  for(i=0; i<nShells; ++i) { if(xprob <= fProbabilities[i]) { break; } }
208 
209  bindingEnergy = elm->GetAtomicShell(i);
210  lv1.set(0.0,0.0,energy,energy);
211  /*
212  G4cout << "nShells= " << nShells << " i= " << i
213  << " Egamma= " << energy << " Ebind= " << bindingEnergy
214  << G4endl;
215  */
216  // for rest frame of the electron
217  G4double x = -G4Log(rndm[1]);
218  eKinEnergy = bindingEnergy*x;
219  ePotEnergy = bindingEnergy*(1.0 + x);
220 
221  // for rest frame of the electron
222  G4double eTotMomentum = sqrt(eKinEnergy*(eKinEnergy + 2*electron_mass_c2));
223  G4double phi = rndm[2]*twopi;
224  G4double costet = 2*rndm[3] - 1;
225  G4double sintet = sqrt((1 - costet)*(1 + costet));
226  lv2.set(eTotMomentum*sintet*cos(phi),eTotMomentum*sintet*sin(phi),
227  eTotMomentum*costet,eKinEnergy + electron_mass_c2);
228  bst = lv2.boostVector();
229  lv1.boost(-bst);
230 
231  gamEnergy0 = lv1.e();
232 
233  // In the rest frame of the electron
234  // The scattered gamma energy is sampled according to Klein-Nishina formula
235  // The random number techniques of Butcher & Messel are used
236  // (Nuc Phys 20(1960),15).
237  G4double E0_m = gamEnergy0/electron_mass_c2;
238 
239  //G4cout << "Nloop= "<< nloop << " Ecm(keV)= " << gamEnergy0/keV << G4endl;
240  //
241  // sample the energy rate of the scattered gamma
242  //
243 
244  G4double epsilon, epsilonsq, onecost, sint2, greject ;
245 
246  G4double eps0 = 1./(1 + 2*E0_m);
247  G4double epsilon0sq = eps0*eps0;
248  G4double alpha1 = - G4Log(eps0);
249  G4double alpha2 = alpha1 + 0.5*(1 - epsilon0sq);
250 
251  do {
252  ++nloop;
253  // false interaction if too many iterations
254  if(nloop > nlooplim) { return; }
255 
256  // 3 random numbers to sample scattering
257  rndmEngineMod->flatArray(3, rndm);
258 
259  if ( alpha1 > alpha2*rndm[0] ) {
260  epsilon = G4Exp(-alpha1*rndm[1]); // epsilon0**r
261  epsilonsq = epsilon*epsilon;
262 
263  } else {
264  epsilonsq = epsilon0sq + (1.- epsilon0sq)*rndm[1];
265  epsilon = sqrt(epsilonsq);
266  }
267 
268  onecost = (1.- epsilon)/(epsilon*E0_m);
269  sint2 = onecost*(2.-onecost);
270  greject = 1. - epsilon*sint2/(1.+ epsilonsq);
271 
272  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
273  } while (greject < rndm[2]);
274  gamEnergy1 = epsilon*gamEnergy0;
275 
276  // before scattering total 4-momentum in e- system
277  lv2.set(0.0,0.0,0.0,electron_mass_c2);
278  lv2 += lv1;
279 
280  //
281  // scattered gamma angles. ( Z - axis along the parent gamma)
282  //
283  if(sint2 < 0.0) { sint2 = 0.0; }
284  costet = 1. - onecost;
285  sintet = sqrt(sint2);
286  phi = twopi * rndmEngineMod->flat();
287 
288  // e- recoil
289  //
290  // in rest frame of the electron
291  G4ThreeVector gamDir = lv1.vect().unit();
292  G4ThreeVector v = G4ThreeVector(sintet*cos(phi),sintet*sin(phi),costet);
293  v.rotateUz(gamDir);
294  lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(),gamEnergy1*v.z(),gamEnergy1);
295  lv2 -= lv1;
296  //G4cout<<"Egam(keV)= " << lv1.e()/keV
297  // <<" Ee(keV)= " << (lv2.e()-electron_mass_c2)/keV << G4endl;
298  lv2.boost(bst);
299  eKinEnergy = lv2.e() - electron_mass_c2 - ePotEnergy;
300  //G4cout << "Nloop= " << nloop << " eKinEnergy= " << eKinEnergy << G4endl;
301 
302  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
303  } while ( eKinEnergy < 0.0 );
304 
305  //
306  // update G4VParticleChange for the scattered gamma
307  //
308 
309  lv1.boost(bst);
310  gamEnergy1 = lv1.e();
311  if(gamEnergy1 > lowestSecondaryEnergy) {
312  G4ThreeVector gamDirection1 = lv1.vect().unit();
313  gamDirection1.rotateUz(direction);
314  fParticleChange->ProposeMomentumDirection(gamDirection1);
315  } else {
316  fParticleChange->ProposeTrackStatus(fStopAndKill);
317  gamEnergy1 = 0.0;
318  }
319  fParticleChange->SetProposedKineticEnergy(gamEnergy1);
320 
321  //
322  // kinematic of the scattered electron
323  //
324 
325  if(eKinEnergy > lowestSecondaryEnergy) {
326  G4ThreeVector eDirection = lv2.vect().unit();
327  eDirection.rotateUz(direction);
328  G4DynamicParticle* dp =
329  new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
330  fvect->push_back(dp);
331  } else { eKinEnergy = 0.0; }
332 
333  G4double edep = energy - gamEnergy1 - eKinEnergy;
334  G4double esec = 0.0;
335 
336  // sample deexcitation
337  //
338  if(fAtomDeexcitation) {
339  G4int index = couple->GetIndex();
340  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
341  G4int Z = elm->GetZasInt();
342  G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i);
343  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
344  G4int nbefore = fvect->size();
345  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
346  G4int nafter = fvect->size();
347  //G4cout << "N1= " << nbefore << " N2= " << nafter << G4endl;
348  for (G4int j=nbefore; j<nafter; ++j) {
349  G4double e = ((*fvect)[j])->GetKineticEnergy();
350  if(esec + e > edep) {
351  // correct energy in order to have energy balance
352  e = edep - esec;
353  ((*fvect)[j])->SetKineticEnergy(e);
354  esec += e;
355  /*
356  G4cout << "### G4KleinNishinaModel Edep(eV)= " << edep/eV
357  << " Esec(eV)= " << esec/eV
358  << " E["<< j << "](eV)= " << e/eV
359  << " N= " << nafter
360  << " Z= " << Z << " shell= " << i
361  << " Ebind(keV)= " << bindingEnergy/keV
362  << " Eshell(keV)= " << shell->BindingEnergy()/keV
363  << G4endl;
364  */
365  // delete the rest of secondaries (should not happens)
366  for (G4int jj=nafter-1; jj>j; --jj) {
367  delete (*fvect)[jj];
368  fvect->pop_back();
369  }
370  break;
371  }
372  esec += e;
373  }
374  edep -= esec;
375  }
376  }
377  if(std::abs(energy - gamEnergy1 - eKinEnergy - esec - edep) > eV) {
378  G4cout << "### G4KleinNishinaModel dE(eV)= "
379  << (energy - gamEnergy1 - eKinEnergy - esec - edep)/eV
380  << " shell= " << i
381  << " E(keV)= " << energy/keV
382  << " Ebind(keV)= " << bindingEnergy/keV
383  << " Eg(keV)= " << gamEnergy1/keV
384  << " Ee(keV)= " << eKinEnergy/keV
385  << " Esec(keV)= " << esec/keV
386  << " Edep(keV)= " << edep/keV
387  << G4endl;
388  }
389  // energy balance
390  if(edep > 0.0) {
391  fParticleChange->ProposeLocalEnergyDeposit(edep);
392  }
393 }
394 
395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
396