Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LowEPComptonModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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