ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEPComptonModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LowEPComptonModel.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 // | G4LowEPComptonModel-- Geant4 Monash University |
29 // | low energy Compton scattering model. |
30 // | J. M. C. Brown, Monash University, Australia |
31 // | ## Unpolarised photons only ## |
32 // | |
33 // | |
34 // *********************************************************************
35 // | |
36 // | The following is a Geant4 class to simulate the process of |
37 // | bound electron Compton scattering. General code structure is |
38 // | based on G4LowEnergyCompton.cc and G4LivermoreComptonModel.cc. |
39 // | Algorithms for photon energy, and ejected Compton electron |
40 // | direction taken from: |
41 // | |
42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin, |
43 // | "A low energy bound atomic electron Compton scattering model |
44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014. |
45 // | |
46 // | The author acknowledges the work of the Geant4 collaboration |
47 // | in developing the following algorithms that have been employed |
48 // | or adapeted for the present software: |
49 // | |
50 // | # sampling of photon scattering angle, |
51 // | # target element selection in composite materials, |
52 // | # target shell selection in element, |
53 // | # and sampling of bound electron momentum from Compton profiles. |
54 // | |
55 // *********************************************************************
56 // | |
57 // | History: |
58 // | -------- |
59 // | |
60 // | Nov. 2011 JMCB - First version |
61 // | Feb. 2012 JMCB - Migration to Geant4 9.5 |
62 // | Sep. 2012 JMCB - Final fixes for Geant4 9.6 |
63 // | Feb. 2013 JMCB - Geant4 9.6 FPE fix for bug 1426 |
64 // | Jan. 2015 JMCB - Migration to MT (Based on Livermore |
65 // | implementation) |
66 // | Feb. 2016 JMCB - Geant4 10.2 FPE fix for bug 1676 |
67 // | |
68 // *********************************************************************
69 
70 #include <limits>
71 #include "G4LowEPComptonModel.hh"
72 #include "G4PhysicalConstants.hh"
73 #include "G4SystemOfUnits.hh"
74 #include "G4Electron.hh"
75 #include "G4ParticleChangeForGamma.hh"
76 #include "G4LossTableManager.hh"
77 #include "G4VAtomDeexcitation.hh"
78 #include "G4AtomicShell.hh"
79 #include "G4Gamma.hh"
80 #include "G4ShellData.hh"
81 #include "G4DopplerProfile.hh"
82 #include "G4Log.hh"
83 #include "G4Exp.hh"
84 
85 //****************************************************************************
86 
87 using namespace std;
88 
89 G4int G4LowEPComptonModel::maxZ = 99;
90 G4LPhysicsFreeVector* G4LowEPComptonModel::data[] = {0};
91 G4ShellData* G4LowEPComptonModel::shellData = 0;
92 G4DopplerProfile* G4LowEPComptonModel::profileData = 0;
93 
94 static const G4double ln10 = G4Log(10.);
95 
96 G4LowEPComptonModel::G4LowEPComptonModel(const G4ParticleDefinition*,
97  const G4String& nam)
98  : G4VEmModel(nam),isInitialised(false)
99 {
100  verboseLevel=1 ;
101  // Verbosity scale:
102  // 0 = nothing
103  // 1 = warning for energy non-conservation
104  // 2 = details of energy budget
105  // 3 = calculation of cross sections, file openings, sampling of atoms
106  // 4 = entering in methods
107 
108  if( verboseLevel>1 ) {
109  G4cout << "Low energy photon Compton model is constructed " << G4endl;
110  }
111 
112  //Mark this model as "applicable" for atomic deexcitation
113  SetDeexcitationFlag(true);
114 
115  fParticleChange = 0;
116  fAtomDeexcitation = 0;
117 }
118 
119 //****************************************************************************
120 
121 G4LowEPComptonModel::~G4LowEPComptonModel()
122 {
123  if(IsMaster()) {
124  delete shellData;
125  shellData = 0;
126  delete profileData;
127  profileData = 0;
128  }
129 }
130 
131 //****************************************************************************
132 
133 void G4LowEPComptonModel::Initialise(const G4ParticleDefinition* particle,
134  const G4DataVector& cuts)
135 {
136  if (verboseLevel > 1) {
137  G4cout << "Calling G4LowEPComptonModel::Initialise()" << G4endl;
138  }
139 
140  // Initialise element selector
141 
142  if(IsMaster()) {
143 
144  // Access to elements
145 
146  char* path = std::getenv("G4LEDATA");
147 
148  G4ProductionCutsTable* theCoupleTable =
149  G4ProductionCutsTable::GetProductionCutsTable();
150  G4int numOfCouples = theCoupleTable->GetTableSize();
151 
152  for(G4int i=0; i<numOfCouples; ++i) {
153  const G4Material* material =
154  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
155  const G4ElementVector* theElementVector = material->GetElementVector();
156  G4int nelm = material->GetNumberOfElements();
157 
158  for (G4int j=0; j<nelm; ++j) {
159  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
160  if(Z < 1) { Z = 1; }
161  else if(Z > maxZ){ Z = maxZ; }
162 
163  if( (!data[Z]) ) { ReadData(Z, path); }
164  }
165  }
166 
167  // For Doppler broadening
168  if(!shellData) {
169  shellData = new G4ShellData();
170  shellData->SetOccupancyData();
171  G4String file = "/doppler/shell-doppler";
172  shellData->LoadData(file);
173  }
174  if(!profileData) { profileData = new G4DopplerProfile(); }
175 
176  InitialiseElementSelectors(particle, cuts);
177  }
178 
179  if (verboseLevel > 2) {
180  G4cout << "Loaded cross section files" << G4endl;
181  }
182 
183  if( verboseLevel>1 ) {
184  G4cout << "G4LowEPComptonModel is initialized " << G4endl
185  << "Energy range: "
186  << LowEnergyLimit() / eV << " eV - "
187  << HighEnergyLimit() / GeV << " GeV"
188  << G4endl;
189  }
190 
191  if(isInitialised) { return; }
192 
193  fParticleChange = GetParticleChangeForGamma();
194  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
195  isInitialised = true;
196 }
197 
198 //****************************************************************************
199 
200 void G4LowEPComptonModel::InitialiseLocal(const G4ParticleDefinition*,
201  G4VEmModel* masterModel)
202 {
203  SetElementSelectors(masterModel->GetElementSelectors());
204 }
205 
206 //****************************************************************************
207 
208 void G4LowEPComptonModel::ReadData(size_t Z, const char* path)
209 {
210  if (verboseLevel > 1)
211  {
212  G4cout << "G4LowEPComptonModel::ReadData()"
213  << G4endl;
214  }
215  if(data[Z]) { return; }
216  const char* datadir = path;
217  if(!datadir)
218  {
219  datadir = std::getenv("G4LEDATA");
220  if(!datadir)
221  {
222  G4Exception("G4LowEPComptonModel::ReadData()",
223  "em0006",FatalException,
224  "Environment variable G4LEDATA not defined");
225  return;
226  }
227  }
228 
229  data[Z] = new G4LPhysicsFreeVector();
230 
231  // Activation of spline interpolation
232  data[Z]->SetSpline(false);
233 
234  std::ostringstream ost;
235  ost << datadir << "/livermore/comp/ce-cs-" << Z <<".dat";
236  std::ifstream fin(ost.str().c_str());
237 
238  if( !fin.is_open())
239  {
240  G4ExceptionDescription ed;
241  ed << "G4LowEPComptonModel data file <" << ost.str().c_str()
242  << "> is not opened!" << G4endl;
243  G4Exception("G4LowEPComptonModel::ReadData()",
244  "em0003",FatalException,
245  ed,"G4LEDATA version should be G4EMLOW6.34 or later");
246  return;
247  } else {
248  if(verboseLevel > 3) {
249  G4cout << "File " << ost.str()
250  << " is opened by G4LowEPComptonModel" << G4endl;
251  }
252  data[Z]->Retrieve(fin, true);
253  data[Z]->ScaleVector(MeV, MeV*barn);
254  }
255  fin.close();
256 }
257 
258 //****************************************************************************
259 
260 
261 G4double
262 G4LowEPComptonModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
263  G4double GammaEnergy,
264  G4double Z, G4double,
265  G4double, G4double)
266 {
267  if (verboseLevel > 3) {
268  G4cout << "G4LowEPComptonModel::ComputeCrossSectionPerAtom()"
269  << G4endl;
270  }
271  G4double cs = 0.0;
272 
273  if (GammaEnergy < LowEnergyLimit()) { return 0.0; }
274 
275  G4int intZ = G4lrint(Z);
276  if(intZ < 1 || intZ > maxZ) { return cs; }
277 
278  G4LPhysicsFreeVector* pv = data[intZ];
279 
280  // if element was not initialised
281  // do initialisation safely for MT mode
282  if(!pv)
283  {
284  InitialiseForElement(0, intZ);
285  pv = data[intZ];
286  if(!pv) { return cs; }
287  }
288 
289  G4int n = pv->GetVectorLength() - 1;
290  G4double e1 = pv->Energy(0);
291  G4double e2 = pv->Energy(n);
292 
293  if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*e1)*pv->Value(e1); }
294  else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; }
295  else if(GammaEnergy > e2) { cs = pv->Value(e2)/GammaEnergy; }
296 
297  return cs;
298 }
299 
300 //****************************************************************************
301 
302 void G4LowEPComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
303  const G4MaterialCutsCouple* couple,
304  const G4DynamicParticle* aDynamicGamma,
305  G4double, G4double)
306 {
307 
308  // The scattered gamma energy is sampled according to Klein - Nishina formula.
309  // then accepted or rejected depending on the Scattering Function multiplied
310  // by factor from Klein - Nishina formula.
311  // Expression of the angular distribution as Klein Nishina
312  // angular and energy distribution and Scattering fuctions is taken from
313  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
314  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
315  // data are interpolated while in the article they are fitted.
316  // Reference to the article is from J. Stepanek New Photon, Positron
317  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
318  // TeV (draft).
319  // The random number techniques of Butcher & Messel are used
320  // (Nucl Phys 20(1960),15).
321 
322 
323  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV;
324 
325  if (verboseLevel > 3) {
326  G4cout << "G4LowEPComptonModel::SampleSecondaries() E(MeV)= "
327  << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
328  << G4endl;
329  }
330  // do nothing below the threshold
331  // should never get here because the XS is zero below the limit
332  if (photonEnergy0 < LowEnergyLimit())
333  return ;
334 
335  G4double e0m = photonEnergy0 / electron_mass_c2 ;
336  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
337 
338  // Select randomly one element in the current material
339 
340  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
341  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
342  G4int Z = (G4int)elm->GetZ();
343 
344  G4double LowEPCepsilon0 = 1. / (1. + 2. * e0m);
345  G4double LowEPCepsilon0Sq = LowEPCepsilon0 * LowEPCepsilon0;
346  G4double alpha1 = -std::log(LowEPCepsilon0);
347  G4double alpha2 = 0.5 * (1. - LowEPCepsilon0Sq);
348 
349  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
350 
351  // Sample the energy of the scattered photon
352  G4double LowEPCepsilon;
353  G4double LowEPCepsilonSq;
354  G4double oneCosT;
355  G4double sinT2;
356  G4double gReject;
357 
358  if (verboseLevel > 3) {
359  G4cout << "Started loop to sample gamma energy" << G4endl;
360  }
361 
362  do
363  {
364  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
365  {
366  LowEPCepsilon = G4Exp(-alpha1 * G4UniformRand());
367  LowEPCepsilonSq = LowEPCepsilon * LowEPCepsilon;
368  }
369  else
370  {
371  LowEPCepsilonSq = LowEPCepsilon0Sq + (1. - LowEPCepsilon0Sq) * G4UniformRand();
372  LowEPCepsilon = std::sqrt(LowEPCepsilonSq);
373  }
374 
375  oneCosT = (1. - LowEPCepsilon) / ( LowEPCepsilon * e0m);
376  sinT2 = oneCosT * (2. - oneCosT);
377  G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
378  G4double scatteringFunction = ComputeScatteringFunction(x, Z);
379  gReject = (1. - LowEPCepsilon * sinT2 / (1. + LowEPCepsilonSq)) * scatteringFunction;
380 
381  } while(gReject < G4UniformRand()*Z);
382 
383  G4double cosTheta = 1. - oneCosT;
384  G4double sinTheta = std::sqrt(sinT2);
385  G4double phi = twopi * G4UniformRand();
386  G4double dirx = sinTheta * std::cos(phi);
387  G4double diry = sinTheta * std::sin(phi);
388  G4double dirz = cosTheta ;
389 
390 
391  // Scatter photon energy and Compton electron direction - Method based on:
392  // J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin'
393  // "A low energy bound atomic electron Compton scattering model for Geant4"
394  // NIMB, Vol. 338, 77-88, 2014.
395 
396  // Set constants and initialize scattering parameters
397 
398  const G4double vel_c = c_light / (m/s);
399  const G4double momentum_au_to_nat = halfpi* hbar_Planck / Bohr_radius / (kg*m/s);
400  const G4double e_mass_kg = electron_mass_c2 / c_squared / kg ;
401 
402  const G4int maxDopplerIterations = 1000;
403  G4double bindingE = 0.;
404  G4double pEIncident = photonEnergy0 ;
405  G4double pERecoil = -1.;
406  G4double eERecoil = -1.;
407  G4double e_alpha =0.;
408  G4double e_beta = 0.;
409 
410  G4double CE_emission_flag = 0.;
411  G4double ePAU = -1;
412  G4int shellIdx = 0;
413  G4double u_temp = 0;
414  G4double cosPhiE =0;
415  G4double sinThetaE =0;
416  G4double cosThetaE =0;
417  G4int iteration = 0;
418 
419  if (verboseLevel > 3) {
420  G4cout << "Started loop to sample photon energy and electron direction" << G4endl;
421  }
422 
423  do{
424 
425 
426  // ******************************************
427  // | Determine scatter photon energy |
428  // ******************************************
429 
430  do
431  {
432  iteration++;
433 
434 
435  // ********************************************
436  // | Sample bound electron information |
437  // ********************************************
438 
439  // Select shell based on shell occupancy
440 
441  shellIdx = shellData->SelectRandomShell(Z);
442  bindingE = shellData->BindingEnergy(Z,shellIdx)/MeV;
443 
444 
445  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
446  ePAU = profileData->RandomSelectMomentum(Z,shellIdx);
447 
448  // Convert to SI units
449  G4double ePSI = ePAU * momentum_au_to_nat;
450 
451  //Calculate bound electron velocity and normalise to natural units
452  u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)) )/vel_c;
453 
454  // Sample incident electron direction, amorphous material, to scattering photon scattering plane
455 
456  e_alpha = pi*G4UniformRand();
457  e_beta = twopi*G4UniformRand();
458 
459  // Total energy of system
460 
461  G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
462  G4double systemE = eEIncident + pEIncident;
463 
464 
465  G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
466  G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha));
467  G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
468  G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
469  G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident);
470  pERecoil = (numerator/denominator);
471  eERecoil = systemE - pERecoil;
472  CE_emission_flag = pEIncident - pERecoil;
473  } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
474 
475 // End of recalculation of photon energy with Doppler broadening
476 
477 
478 
479  // *******************************************************
480  // | Determine ejected Compton electron direction |
481  // *******************************************************
482 
483  // Calculate velocity of ejected Compton electron
484 
485  G4double a_temp = eERecoil / electron_mass_c2;
486  G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
487 
488  // Coefficients and terms from simulatenous equations
489 
490  G4double sinAlpha = std::sin(e_alpha);
491  G4double cosAlpha = std::cos(e_alpha);
492  G4double sinBeta = std::sin(e_beta);
493  G4double cosBeta = std::cos(e_beta);
494 
495  G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
496  G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
497 
498  G4double var_A = pERecoil*u_p_temp*sinTheta;
499  G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
500  G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta);
501 
502  G4double var_D1 = gamma*electron_mass_c2*pERecoil;
503  G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
504  G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil);
505  G4double var_D = var_D1*var_D2 + var_D3;
506 
507  G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha);
508  G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
509  G4double var_E = var_E1 - var_E2;
510 
511  G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
512  G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
513  G4double var_F = var_F1 - var_F2;
514 
515  G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
516 
517  // Two equations form a quadratic form of Wx^2 + Yx + Z = 0
518  // Coefficents and solution to quadratic
519 
520  G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
521  G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
522  G4double var_W = var_W1 + var_W2;
523 
524  G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C));
525 
526  G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
527  G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
528  G4double var_Z = var_Z1 + var_Z2;
529  G4double diff1 = var_Y*var_Y;
530  G4double diff2 = 4*var_W*var_Z;
531  G4double diff = diff1 - diff2;
532 
533 
534  // Check if diff is less than zero, if so ensure it is due to FPE
535 
536  //Determine number of digits (in decimal base) that G4double can accurately represent
537  G4double g4d_order = G4double(numeric_limits<G4double>::digits10);
538  G4double g4d_limit = std::pow(10.,-g4d_order);
539  //Confirm that diff less than zero is due FPE, i.e if abs of diff / diff1 and diff/ diff2 is less
540  //than 10^(-g4d_order), then set diff to zero
541 
542  if ((diff < 0.0) && (abs(diff / diff1) < g4d_limit) && (abs(diff / diff2) < g4d_limit) )
543  {
544  diff = 0.0;
545  }
546 
547 
548  // Plus and minus of quadratic
549  G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
550  G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
551 
552 
553  // Floating point precision protection
554  // Check if X_p and X_m are greater than or less than 1 or -1, if so clean up FPE
555  // Issue due to propagation of FPE and only impacts 8th sig fig onwards
556 
557  if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;}
558  if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;}
559 
560  // End of FP protection
561 
562  G4double ThetaE = 0.;
563 
564 
565  // Randomly sample one of the two possible solutions and determin theta angle of ejected Compton electron
566  G4double sol_select = G4UniformRand();
567 
568  if (sol_select < 0.5)
569  {
570  ThetaE = std::acos(X_p);
571  }
572  if (sol_select > 0.5)
573  {
574  ThetaE = std::acos(X_m);
575  }
576 
577 
578  cosThetaE = std::cos(ThetaE);
579  sinThetaE = std::sin(ThetaE);
580  G4double Theta = std::acos(cosTheta);
581 
582  //Calculate electron Phi
583  G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE));
584  G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta));
585  G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
586  // Trigs
587  cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
588 
589  // End of calculation of ejection Compton electron direction
590 
591  //Fix for floating point errors
592 
593  } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1));
594 
595  // Revert to original if maximum number of iterations threshold has been reached
596  if (iteration >= maxDopplerIterations)
597  {
598  pERecoil = photonEnergy0 ;
599  bindingE = 0.;
600  dirx=0.0;
601  diry=0.0;
602  dirz=1.0;
603  }
604 
605  // Set "scattered" photon direction and energy
606 
607  G4ThreeVector photonDirection1(dirx,diry,dirz);
608  photonDirection1.rotateUz(photonDirection0);
609  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
610 
611  if (pERecoil > 0.)
612  {
613  fParticleChange->SetProposedKineticEnergy(pERecoil) ;
614 
615  // Set ejected Compton electron direction and energy
616  G4double PhiE = std::acos(cosPhiE);
617  G4double eDirX = sinThetaE * std::cos(phi+PhiE);
618  G4double eDirY = sinThetaE * std::sin(phi+PhiE);
619  G4double eDirZ = cosThetaE;
620 
621  G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
622 
623  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
624  eDirection.rotateUz(photonDirection0);
625  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
626  eDirection,eKineticEnergy) ;
627  fvect->push_back(dp);
628 
629  }
630  else
631  {
632  fParticleChange->SetProposedKineticEnergy(0.);
633  fParticleChange->ProposeTrackStatus(fStopAndKill);
634  }
635 
636  // sample deexcitation
637  //
638 
639  if (verboseLevel > 3) {
640  G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl;
641  }
642 
643  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
644  G4int index = couple->GetIndex();
645  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
646  size_t nbefore = fvect->size();
647  G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
648  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
649  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
650  size_t nafter = fvect->size();
651  if(nafter > nbefore) {
652  for (size_t i=nbefore; i<nafter; ++i) {
653  //Check if there is enough residual energy
654  if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
655  {
656  //Ok, this is a valid secondary: keep it
657  bindingE -= ((*fvect)[i])->GetKineticEnergy();
658  }
659  else
660  {
661  //Invalid secondary: not enough energy to create it!
662  //Keep its energy in the local deposit
663  delete (*fvect)[i];
664  (*fvect)[i]=0;
665  }
666  }
667  }
668  }
669  }
670 
671  //This should never happen
672  if(bindingE < 0.0)
673  G4Exception("G4LowEPComptonModel::SampleSecondaries()",
674  "em2051",FatalException,"Negative local energy deposit");
675 
676  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
677 
678 }
679 
680 //****************************************************************************
681 
682 G4double
683 G4LowEPComptonModel::ComputeScatteringFunction(G4double x, G4int Z)
684 {
685  G4double value = Z;
686  if (x <= ScatFuncFitParam[Z][2]) {
687 
688  G4double lgq = G4Log(x)/ln10;
689 
690  if (lgq < ScatFuncFitParam[Z][1]) {
691  value = ScatFuncFitParam[Z][3] + lgq*ScatFuncFitParam[Z][4];
692  } else {
693  value = ScatFuncFitParam[Z][5] + lgq*ScatFuncFitParam[Z][6] +
694  lgq*lgq*ScatFuncFitParam[Z][7] + lgq*lgq*lgq*ScatFuncFitParam[Z][8];
695  }
696  value = G4Exp(value*ln10);
697  }
698  return value;
699 }
700 
701 
702 //****************************************************************************
703 
704 #include "G4AutoLock.hh"
705 namespace { G4Mutex LowEPComptonModelMutex = G4MUTEX_INITIALIZER; }
706 
707 void
708 G4LowEPComptonModel::InitialiseForElement(const G4ParticleDefinition*,
709  G4int Z)
710 {
711  G4AutoLock l(&LowEPComptonModelMutex);
712  if(!data[Z]) { ReadData(Z); }
713  l.unlock();
714 }
715 
716 //****************************************************************************
717 
718 //Fitting data to compute scattering function
719 
720 const G4double G4LowEPComptonModel::ScatFuncFitParam[101][9] = {
721 { 0, 0., 0., 0., 0., 0., 0., 0., 0.},
722 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -143.374, 50.787, -5.951, 0.2304418},
723 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -53.649, 13.892, -0.948, 0.006996759},
724 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -62.090, 21.462, -2.453, 0.093416},
725 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -127.906, 46.491, -5.614, 0.2262103},
726 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -131.153, 47.132, -5.619, 0.2233819},
727 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -128.143, 45.379, -5.325, 0.2083009},
728 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -131.048, 46.314, -5.421, 0.2114925},
729 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -128.225, 44.818, -5.183, 0.1997155},
730 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -122.112, 42.103, -4.796, 0.1819099},
731 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -110.143, 37.225, -4.143, 0.1532094},
732 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -41.137, 12.313, -1.152, 0.03384553},
733 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -53.549, 17.420, -1.840, 0.06431849},
734 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -66.243, 22.297, -2.460, 0.09045854},
735 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -78.271, 26.757, -3.008, 0.1128195},
736 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -85.069, 29.164, -3.291, 0.1239113},
737 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -93.640, 32.274, -3.665, 0.1388633},
738 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -98.534, 33.958, -3.857, 0.1461557},
739 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -100.077, 34.379, -3.891, 0.1468902},
740 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -53.819, 17.528, -1.851, 0.0648722},
741 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -52.221, 17.169, -1.832, 0.06502094},
742 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -55.365, 18.276, -1.961, 0.07002778},
743 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -57.412, 18.957, -2.036, 0.07278856},
744 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -58.478, 19.270, -2.065, 0.07362722},
745 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -62.192, 20.358, -2.167, 0.07666583},
746 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -58.007, 18.924, -2.003, 0.0704305},
747 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -61.176, 20.067, -2.141, 0.0760269},
748 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -61.909, 20.271, -2.159, 0.07653559},
749 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -62.402, 20.391, -2.167, 0.07664243},
750 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -67.305, 21.954, -2.331, 0.0823267},
751 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -62.064, 20.136, -2.122, 0.07437589},
752 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -61.068, 19.808, -2.086, 0.07307488},
753 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -63.126, 20.553, -2.175, 0.07660222},
754 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -65.674, 21.445, -2.278, 0.08054694},
755 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -69.457, 22.811, -2.442, 0.08709536},
756 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -72.283, 23.808, -2.558, 0.09156808},
757 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -74.696, 24.641, -2.653, 0.09516597},
758 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -46.235, 14.519, -1.458, 0.04837659},
759 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -41.784, 13.065, -1.300, 0.04267703},
760 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -44.609, 14.114, -1.429, 0.0479348},
761 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -47.142, 15.019, -1.536, 0.0521347},
762 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -53.114, 17.052, -1.766, 0.06079426},
763 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -54.590, 17.550, -1.822, 0.06290335},
764 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -51.272, 16.423, -1.694, 0.05806108},
765 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -58.314, 18.839, -1.969, 0.0684608},
766 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -59.674, 19.295, -2.020, 0.07037294},
767 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -72.228, 23.693, -2.532, 0.09017969},
768 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -60.890, 19.647, -2.053, 0.07138694},
769 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -56.152, 18.002, -1.863, 0.06410123},
770 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -56.208, 18.052, -1.872, 0.06456884},
771 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -58.478, 18.887, -1.973, 0.06860356},
772 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -60.708, 19.676, -2.066, 0.07225841},
773 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -63.341, 20.632, -2.180, 0.0767412},
774 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -66.339, 21.716, -2.310, 0.08191981},
775 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -67.649, 22.151, -2.357, 0.08357825},
776 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -45.302, 14.219, -1.423, 0.04712317},
777 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -39.825, 12.363, -1.214, 0.03931009},
778 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -38.952, 11.982, -1.160, 0.03681554},
779 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -41.959, 13.118, -1.302, 0.04271291},
780 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -40.528, 12.555, -1.230, 0.03971407},
781 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -39.986, 12.329, -1.200, 0.03843737},
782 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -42.756, 13.362, -1.327, 0.0436124},
783 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -42.682, 13.314, -1.319, 0.04322509},
784 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -42.399, 13.185, -1.301, 0.04243861},
785 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -43.226, 13.475, -1.335, 0.04377341},
786 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -43.232, 13.456, -1.330, 0.04347536},
787 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -42.047, 12.990, -1.270, 0.04095499},
788 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -42.405, 13.106, -1.283, 0.04146024},
789 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -41.974, 12.926, -1.259, 0.040435},
790 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -42.132, 12.967, -1.262, 0.04048908},
791 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -42.582, 13.122, -1.280, 0.04119599},
792 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -42.586, 13.115, -1.278, 0.04107587},
793 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -43.708, 13.509, -1.324, 0.04286491},
794 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -44.838, 13.902, -1.369, 0.04457132},
795 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -45.545, 14.137, -1.395, 0.04553459},
796 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -46.426, 14.431, -1.428, 0.04678218},
797 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -47.513, 14.806, -1.471, 0.04842566},
798 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -48.225, 15.042, -1.497, 0.04938364},
799 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -53.256, 16.739, -1.687, 0.05645173},
800 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -53.900, 16.946, -1.709, 0.05723134},
801 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -49.801, 15.536, -1.547, 0.05103522},
802 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -49.651, 15.512, -1.548, 0.05123203},
803 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -51.021, 16.018, -1.609, 0.05364831},
804 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -52.693, 16.612, -1.679, 0.05638698},
805 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -54.415, 17.238, -1.754, 0.05935566},
806 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -56.083, 17.834, -1.824, 0.06206068},
807 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -57.860, 18.463, -1.898, 0.0649633},
808 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -39.973, 12.164, -1.162, 0.0364598},
809 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -34.591, 10.338, -0.956, 0.0287409},
810 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -36.002, 10.867, -1.021, 0.03136835},
811 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -37.660, 11.475, -1.095, 0.03435334},
812 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -37.268, 11.301, -1.071, 0.0330539},
813 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -37.695, 11.438, -1.085, 0.03376669},
814 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -39.010, 11.927, -1.146, 0.03630848},
815 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -37.192, 11.229, -1.057, 0.0325621},
816 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -37.589, 11.363, -1.072, 0.03312393},
817 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -39.573, 12.095, -1.162, 0.03680527},
818 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -40.007, 12.242, -1.178, 0.03737377},
819 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -39.509, 12.041, -1.152, 0.03629023},
820 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -39.939, 12.183, -1.168, 0.03690464},
821 {100, 6.412, 1.00000E+10, -12.900, 1.993, -39.973, 12.180, -1.166, 0.036773}
822  };
823 
824 
825