Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopeComptonModel.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 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4PenelopeComptonModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4PenelopeComptonModel.cc (Version 9.3.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4PenelopeComptonModel.cc,v 1.8 2009/10/23 09:29:24 pandola Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-03-patch-01 $
 26 //                                                 28 //
 27 // Author: Luciano Pandola                         29 // Author: Luciano Pandola
 28 //                                                 30 //
 29 // History:                                        31 // History:
 30 // --------                                        32 // --------
 31 // 15 Feb 2010   L Pandola  Implementation     <<  33 // 02 Oct 2008   L Pandola    Migration from process to model 
 32 // 18 Mar 2010   L Pandola  Removed GetAtomsPe <<  34 // 28 Oct 2008   L Pandola    Treat the database data from Penelope according to the 
 33 //                            to G4PenelopeOsc <<  35 //                            original model, namely merging levels below 15 eV in 
 34 // 01 Feb 2011   L Pandola  Suppress fake ener <<  36 //                            a single one. Still, it is not fully compliant with the 
 35 //                            active.          <<  37 //                            original Penelope model, because plasma excitation is not 
 36 //                          Make sure that flu <<  38 //                            considered.
 37 //                            if above thresho <<  39 // 22 Nov 2008   L Pandola    Make unit of measurements explicit for binding energies 
 38 // 24 May 2011   L Pandola  Renamed (make v200 <<  40 //            that are read from the external files.
 39 // 10 Jun 2011   L Pandola  Migrate atomic dee <<  41 // 24 Nov 2008   L Pandola    Find a cleaner way to delete vectors.
 40 // 09 Oct 2013   L Pandola  Migration to MT    <<  42 // 16 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
 41 // 25 Jul 2023   D Iuso     Fix for possible i <<  43 //                  - apply internal high-energy limit only in constructor 
                                                   >>  44 //                  - do not apply low-energy limit (default is 0)
                                                   >>  45 //                  - do not apply production threshold on level of the model
                                                   >>  46 // 21 Oct 2009   L Pandola    Remove un-necessary fUseAtomicDeexcitation flag - now managed by
                                                   >>  47 //                            G4VEmModel::DeexcitationFlag()
                                                   >>  48 //                            Add ActivateAuger() method
 42 //                                                 49 //
                                                   >>  50 
 43 #include "G4PenelopeComptonModel.hh"               51 #include "G4PenelopeComptonModel.hh"
 44 #include "G4PhysicalConstants.hh"              << 
 45 #include "G4SystemOfUnits.hh"                  << 
 46 #include "G4ParticleDefinition.hh"                 52 #include "G4ParticleDefinition.hh"
 47 #include "G4MaterialCutsCouple.hh"                 53 #include "G4MaterialCutsCouple.hh"
                                                   >>  54 #include "G4ProductionCutsTable.hh"
 48 #include "G4DynamicParticle.hh"                    55 #include "G4DynamicParticle.hh"
 49 #include "G4VEMDataSet.hh"                         56 #include "G4VEMDataSet.hh"
 50 #include "G4PhysicsTable.hh"                       57 #include "G4PhysicsTable.hh"
                                                   >>  58 #include "G4ElementTable.hh"
                                                   >>  59 #include "G4Element.hh"
 51 #include "G4PhysicsLogVector.hh"                   60 #include "G4PhysicsLogVector.hh"
                                                   >>  61 #include "G4PenelopeIntegrator.hh"
 52 #include "G4AtomicTransitionManager.hh"            62 #include "G4AtomicTransitionManager.hh"
 53 #include "G4AtomicShell.hh"                        63 #include "G4AtomicShell.hh"
 54 #include "G4Gamma.hh"                              64 #include "G4Gamma.hh"
 55 #include "G4Electron.hh"                           65 #include "G4Electron.hh"
 56 #include "G4PenelopeOscillatorManager.hh"      << 
 57 #include "G4PenelopeOscillator.hh"             << 
 58 #include "G4LossTableManager.hh"               << 
 59 #include "G4Exp.hh"                            << 
 60                                                    66 
 61 //....oooOO0OOooo........oooOO0OOooo........oo     67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 62                                                    68 
 63 G4PenelopeComptonModel::G4PenelopeComptonModel <<  69 
 64                  const G4String& nam)          <<  70 G4PenelopeComptonModel::G4PenelopeComptonModel(const G4ParticleDefinition*,
 65   :G4VEmModel(nam),fParticleChange(nullptr),fP <<  71                                              const G4String& nam)
 66    fAtomDeexcitation(nullptr),                 <<  72   :G4VEmModel(nam),ionizationEnergy(0),hartreeFunction(0),
 67    fOscManager(nullptr),fIsInitialised(false)  <<  73    occupationNumber(0),isInitialised(false)
 68 {                                                  74 {
 69   fIntrinsicLowEnergyLimit = 100.0*eV;             75   fIntrinsicLowEnergyLimit = 100.0*eV;
 70   fIntrinsicHighEnergyLimit = 100.0*GeV;           76   fIntrinsicHighEnergyLimit = 100.0*GeV;
                                                   >>  77   //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
 71   SetHighEnergyLimit(fIntrinsicHighEnergyLimit     78   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
 72   //                                               79   //
 73   fOscManager = G4PenelopeOscillatorManager::G <<  80   energyForIntegration = 0.0;
                                                   >>  81   ZForIntegration = 1;
 74                                                    82 
 75   if (part)                                    <<  83   //by default, the model will use atomic deexcitation
 76     SetParticle(part);                         <<  84   SetDeexcitationFlag(true);
                                                   >>  85   ActivateAuger(false);
 77                                                    86 
 78   fVerboseLevel= 0;                            <<  87   verboseLevel= 0;
 79   // Verbosity scale:                              88   // Verbosity scale:
 80   // 0 = nothing                               <<  89   // 0 = nothing 
 81   // 1 = warning for energy non-conservation   <<  90   // 1 = warning for energy non-conservation 
 82   // 2 = details of energy budget                  91   // 2 = details of energy budget
 83   // 3 = calculation of cross sections, file o     92   // 3 = calculation of cross sections, file openings, sampling of atoms
 84   // 4 = entering in methods                       93   // 4 = entering in methods
 85                                                    94 
 86   //Mark this model as "applicable" for atomic <<  95   //These vectors do not change when materials or cut change.
 87   SetDeexcitationFlag(true);                   <<  96   //Therefore I can read it at the constructor
                                                   >>  97   ionizationEnergy = new std::map<G4int,G4DataVector*>;
                                                   >>  98   hartreeFunction  = new std::map<G4int,G4DataVector*>;
                                                   >>  99   occupationNumber = new std::map<G4int,G4DataVector*>;
                                                   >> 100 
                                                   >> 101   ReadData(); //Read data from file
 88                                                   102 
 89   fTransitionManager = G4AtomicTransitionManag << 
 90 }                                                 103 }
 91                                                   104 
 92 //....oooOO0OOooo........oooOO0OOooo........oo    105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 93                                                   106 
 94 G4PenelopeComptonModel::~G4PenelopeComptonMode    107 G4PenelopeComptonModel::~G4PenelopeComptonModel()
 95 {;}                                            << 108 {  
                                                   >> 109   std::map <G4int,G4DataVector*>::iterator i;
                                                   >> 110   for (i=ionizationEnergy->begin();i != ionizationEnergy->end();i++)
                                                   >> 111     if (i->second) delete i->second;
                                                   >> 112   for (i=hartreeFunction->begin();i != hartreeFunction->end();i++)
                                                   >> 113     if (i->second) delete i->second;
                                                   >> 114   for (i=occupationNumber->begin();i != occupationNumber->end();i++)
                                                   >> 115     if (i->second) delete i->second;
                                                   >> 116 
                                                   >> 117 
                                                   >> 118   if (ionizationEnergy)
                                                   >> 119     delete ionizationEnergy;
                                                   >> 120   if (hartreeFunction)
                                                   >> 121     delete hartreeFunction;
                                                   >> 122   if (occupationNumber)
                                                   >> 123     delete occupationNumber;
                                                   >> 124 }
 96                                                   125 
 97 //....oooOO0OOooo........oooOO0OOooo........oo    126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 98                                                   127 
 99 void G4PenelopeComptonModel::Initialise(const  << 128 void G4PenelopeComptonModel::Initialise(const G4ParticleDefinition* particle,
100             const G4DataVector&)               << 129           const G4DataVector& cuts)
101 {                                                 130 {
102   if (fVerboseLevel > 3)                       << 131   if (verboseLevel > 3)
103     G4cout << "Calling G4PenelopeComptonModel:    132     G4cout << "Calling G4PenelopeComptonModel::Initialise()" << G4endl;
104                                                   133 
105   fAtomDeexcitation = G4LossTableManager::Inst << 134   InitialiseElementSelectors(particle,cuts);
106   //Issue warning if the AtomicDeexcitation ha << 
107   if (!fAtomDeexcitation)                      << 
108     {                                          << 
109       G4cout << G4endl;                        << 
110       G4cout << "WARNING from G4PenelopeCompto << 
111       G4cout << "Atomic de-excitation module i << 
112       G4cout << "any fluorescence/Auger emissi << 
113       G4cout << "Please make sure this is inte << 
114     }                                          << 
115                                                << 
116   SetParticle(part);                           << 
117                                                << 
118   if (IsMaster() && part == fParticle)         << 
119     {                                          << 
120                                                   135 
121       if (fVerboseLevel > 0)                   << 136   if (verboseLevel > 0) {
122   {                                            << 137     G4cout << "Penelope Compton model is initialized " << G4endl
123     G4cout << "Penelope Compton model v2008 is << 138      << "Energy range: "
124      << "Energy range: "                       << 139      << LowEnergyLimit() / keV << " keV - "
125      << LowEnergyLimit() / keV << " keV - "    << 140      << HighEnergyLimit() / GeV << " GeV"
126      << HighEnergyLimit() / GeV << " GeV";     << 141      << G4endl;
127   }                                            << 142   }
128       //Issue a warning, if the model is going << 
129       //energy which is outside the validity o << 
130       if (LowEnergyLimit() < fIntrinsicLowEner << 
131   {                                            << 
132     G4ExceptionDescription ed;                 << 
133     ed << "Using the Penelope Compton model ou << 
134        << G4endl;                              << 
135     ed << "-> LowEnergyLimit() in process = "  << 
136     ed << "-> Instrinsic low-energy limit = "  << 
137        << G4endl;                              << 
138     ed << "Result of the simulation have to be << 
139     G4Exception("G4PenelopeComptonModel::Initi << 
140           "em2100",JustWarning,ed);            << 
141   }                                            << 
142     }                                          << 
143                                                   143 
144   if(fIsInitialised) return;                   << 144   if(isInitialised) return;
145   fParticleChange = GetParticleChangeForGamma(    145   fParticleChange = GetParticleChangeForGamma();
146   fIsInitialised = true;                       << 146   isInitialised = true; 
147                                                << 
148 }                                              << 
149                                                << 
150 //....oooOO0OOooo........oooOO0OOooo........oo << 
151                                                << 
152 void G4PenelopeComptonModel::InitialiseLocal(c << 
153                                                << 
154 {                                              << 
155   if (fVerboseLevel > 3)                       << 
156     G4cout << "Calling  G4PenelopeComptonModel << 
157   //                                           << 
158   //Check that particle matches: one might hav << 
159   //for e+ and e-).                            << 
160   //                                           << 
161   if (part == fParticle)                       << 
162     {                                          << 
163       //Get the const table pointers from the  << 
164       const G4PenelopeComptonModel* theModel = << 
165         static_cast<G4PenelopeComptonModel*> ( << 
166                                                << 
167       //Same verbosity for all workers, as the << 
168       fVerboseLevel = theModel->fVerboseLevel; << 
169     }                                          << 
170   return;                                      << 
171 }                                                 147 }
172                                                   148 
173                                                << 
174 //....oooOO0OOooo........oooOO0OOooo........oo    149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
175                                                   150 
176 G4double G4PenelopeComptonModel::CrossSectionP << 151 G4double G4PenelopeComptonModel::ComputeCrossSectionPerAtom(
177                                                << 152                                        const G4ParticleDefinition*,
178                                                << 153                                              G4double energy,
179                                                << 154                                              G4double Z, G4double,
180                                                << 155                                              G4double, G4double)
181 {                                                 156 {
182   // Penelope model v2008 to calculate the Com << 157   // Penelope model to calculate the Compton scattering cross section:
183   // D. Brusa et al., Nucl. Instrum. Meth. A 3    158   // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
184   //                                           << 159   // 
185   // The cross section for Compton scattering  << 160   // The cross section for Compton scattering is calculated according to the Klein-Nishina 
186   // formula for energy > 5 MeV.                  161   // formula for energy > 5 MeV.
187   // For E < 5 MeV it is used a parametrizatio    162   // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega,
188   // which is integrated numerically in cos(th    163   // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection().
189   // The parametrization includes the J(p)     << 164   // The parametrization includes the J(p) 
190   // distribution profiles for the atomic shel << 165   // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations 
191   // from F. Biggs et al., At. Data Nucl. Data    166   // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
192   //                                              167   //
193   if (fVerboseLevel > 3)                       << 168   if (verboseLevel > 3)
194     G4cout << "Calling CrossSectionPerVolume() << 169     G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeComptonModel" << G4endl;
195   SetupForMaterial(p, material, energy);       << 
196                                                << 
197   G4double cs = 0;                             << 
198   //Force null cross-section if below the low- << 
199   if (energy < LowEnergyLimit())               << 
200     return cs;                                 << 
201                                                   170 
202   //Retrieve the oscillator table for this mat << 171   G4int iZ = (G4int) Z;
203   G4PenelopeOscillatorTable* theTable = fOscMa << 172   G4double cs=0.0;
204                                                << 173   energyForIntegration=energy; 
205   if (energy < 5*MeV) //explicit calculation f << 174   ZForIntegration = iZ;
206     {                                          << 175   if (energy< 5*MeV)
207       size_t numberOfOscillators = theTable->s << 176     {
208       for (size_t i=0;i<numberOfOscillators;i+ << 177       // numerically integrate differential cross section dSigma/dOmega
209   {                                            << 178       G4PenelopeIntegrator<G4PenelopeComptonModel,G4double (G4PenelopeComptonModel::*)(G4double)> 
210     G4PenelopeOscillator* theOsc = (*theTable) << 179   theIntegrator;
211     //sum contributions coming from each oscil << 180       cs = theIntegrator.Calculate(this,&G4PenelopeComptonModel::DifferentialCrossSection,-1.0,1.0,1e-05);
212     cs += OscillatorTotalCrossSection(energy,t << 181     }
                                                   >> 182   else 
                                                   >> 183     {
                                                   >> 184       // use Klein-Nishina formula
                                                   >> 185       G4double ki=energy/electron_mass_c2;
                                                   >> 186       G4double ki3=ki*ki;
                                                   >> 187       G4double ki2=1.0+2*ki;
                                                   >> 188       G4double ki1=ki3-ki2-1.0;
                                                   >> 189       G4double t0=1.0/(ki2);
                                                   >> 190       G4double csl = 0.5*ki3*t0*t0+ki2*t0+ki1*std::log(t0)-(1.0/t0);
                                                   >> 191       G4int nosc = occupationNumber->find(iZ)->second->size();
                                                   >> 192       for (G4int i=0;i<nosc;i++)
                                                   >> 193   {
                                                   >> 194     G4double ionEnergy = (*(ionizationEnergy->find(iZ)->second))[i];
                                                   >> 195     G4double tau=(energy-ionEnergy)/energy;
                                                   >> 196     if (tau > t0)
                                                   >> 197       {
                                                   >> 198         G4double csu = 0.5*ki3*tau*tau+ki2*tau+ki1*std::log(tau)-(1.0/tau);
                                                   >> 199         G4int f = (G4int) (*(occupationNumber->find(iZ)->second))[i];
                                                   >> 200         cs = cs + f*(csu-csl);
                                                   >> 201       }
213   }                                               202   }
                                                   >> 203       cs=pi*classic_electr_radius*classic_electr_radius*cs/(ki*ki3);
214     }                                             204     }
215   else //use Klein-Nishina for E>5 MeV         << 205   
216     cs = KleinNishinaCrossSection(energy,mater << 206   if (verboseLevel > 2)
217                                                << 207     G4cout << "Compton cross Section at " << energy/keV << " keV for Z=" << Z << 
218   //cross sections are in units of pi*classic_ << 208       " = " << cs/barn << " barn" << G4endl;
219   cs *= pi*classic_electr_radius*classic_elect << 209   return cs;
220                                                << 
221   //Now, cs is the cross section *per molecule << 
222   //cross section per volume                   << 
223   G4double atomDensity = material->GetTotNbOfA << 
224   G4double atPerMol =  fOscManager->GetAtomsPe << 
225                                                << 
226   if (fVerboseLevel > 3)                       << 
227     G4cout << "Material " << material->GetName << 
228       "atoms per molecule" << G4endl;          << 
229                                                << 
230   G4double moleculeDensity = 0.;               << 
231                                                << 
232   if (atPerMol)                                << 
233     moleculeDensity = atomDensity/atPerMol;    << 
234                                                << 
235   G4double csvolume = cs*moleculeDensity;      << 
236                                                << 
237   if (fVerboseLevel > 2)                       << 
238     G4cout << "Compton mean free path at " <<  << 
239             material->GetName() << " = " << (1 << 
240   return csvolume;                             << 
241 }                                              << 
242                                                << 
243 //....oooOO0OOooo........oooOO0OOooo........oo << 
244                                                << 
245 //This is a dummy method. Never inkoved by the << 
246 //a warning if one tries to get Cross Sections << 
247 //G4EmCalculator.                              << 
248 G4double G4PenelopeComptonModel::ComputeCrossS << 
249                                                << 
250                                                << 
251                                                << 
252                                                << 
253                                                << 
254 {                                              << 
255   G4cout << "*** G4PenelopeComptonModel -- WAR << 
256   G4cout << "Penelope Compton model v2008 does << 
257   G4cout << "so the result is always zero. For << 
258   G4cout << "GetCrossSectionPerVolume() or Get << 
259   return 0;                                    << 
260 }                                                 210 }
261                                                   211 
262 //....oooOO0OOooo........oooOO0OOooo........oo    212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263                                                   213 
264 void G4PenelopeComptonModel::SampleSecondaries    214 void G4PenelopeComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
265                  const G4MaterialCutsCouple* c << 215                 const G4MaterialCutsCouple* couple,
266                 const G4DynamicParticle* aDyna    216                 const G4DynamicParticle* aDynamicGamma,
267                 G4double,                         217                 G4double,
268                 G4double)                         218                 G4double)
269 {                                                 219 {
270   // Penelope model v2008 to sample the Compto << 220   
                                                   >> 221   // Penelope model to sample the Compton scattering final state.
271   // D. Brusa et al., Nucl. Instrum. Meth. A 3    222   // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
272   // The model determines also the original sh << 223   // The model determines also the original shell from which the electron is expelled, 
273   // in order to produce fluorescence de-excit    224   // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
274   //                                           << 225   // 
275   // The final state for Compton scattering is << 226   // The final state for Compton scattering is calculated according to the Klein-Nishina 
276   // formula for energy > 5 MeV. In this case, << 227   // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and 
277   // one can assume that the target electron i    228   // one can assume that the target electron is at rest.
278   // For E < 5 MeV it is used the parametrizat    229   // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega,
279   // to sample the scattering angle and the en << 230   // to sample the scattering angle and the energy of the emerging electron, which is  
280   // G4PenelopeComptonModel::DifferentialCross << 231   // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is 
281   // used to sample cos(theta). The efficiency << 232   // used to sample cos(theta). The efficiency increases monotonically with photon energy and is 
282   // nearly independent on the Z; typical valu << 233   // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV, 
283   // respectively.                                234   // respectively.
284   // The parametrization includes the J(p) dis << 235   // The parametrization includes the J(p) distribution profiles for the atomic shells, that are 
285   // tabulated                                 << 236   // tabulated 
286   // from Hartree-Fock calculations from F. Bi << 237   // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201. 
287   // Doppler broadening is included.              238   // Doppler broadening is included.
288   //                                              239   //
289                                                   240 
290   if (fVerboseLevel > 3)                       << 241   if (verboseLevel > 3)
291     G4cout << "Calling SampleSecondaries() of     242     G4cout << "Calling SampleSecondaries() of G4PenelopeComptonModel" << G4endl;
292                                                   243 
293   G4double photonEnergy0 = aDynamicGamma->GetK    244   G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
294                                                   245 
295   // do nothing below the threshold            << 246   if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
296   // should never get here because the XS is z << 247     {
297   if(photonEnergy0 < LowEnergyLimit())         << 248       fParticleChange->ProposeTrackStatus(fStopAndKill);
298     return;                                    << 249       fParticleChange->SetProposedKineticEnergy(0.);
                                                   >> 250       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
                                                   >> 251       return ;
                                                   >> 252     }
299                                                   253 
300   G4ParticleMomentum photonDirection0 = aDynam    254   G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
301   const G4Material* material = couple->GetMate << 
302                                                   255 
303   G4PenelopeOscillatorTable* theTable = fOscMa << 256   // Select randomly one element in the current material
                                                   >> 257   if (verboseLevel > 2)
                                                   >> 258     G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
                                                   >> 259   // atom can be selected effitiantly if element selectors are initialised 
                                                   >> 260   const G4Element* anElement = 
                                                   >> 261     SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy0);
                                                   >> 262   G4int Z = (G4int) anElement->GetZ();
                                                   >> 263   if (verboseLevel > 2)
                                                   >> 264     G4cout << "Selected " << anElement->GetName() << G4endl;
304                                                   265 
305   const G4int nmax = 64;                          266   const G4int nmax = 64;
306   G4double rn[nmax]={0.0};                     << 267   G4double rn[nmax],pac[nmax];
307   G4double pac[nmax]={0.0};                    << 268   
308                                                << 269   G4double ki,ki1,ki2,ki3,taumin,a1,a2;
                                                   >> 270   G4double tau,TST;
309   G4double S=0.0;                                 271   G4double S=0.0;
310   G4double epsilon = 0.0;                      << 272   G4double epsilon,cosTheta;
311   G4double cosTheta = 1.0;                     << 273   G4double harFunc = 0.0;
312   G4double hartreeFunc = 0.0;                  << 274   G4int occupNb= 0;
313   G4double oscStren = 0.0;                     << 275   G4double ionEnergy=0.0;
314   size_t numberOfOscillators = theTable->size( << 276   G4int nosc = occupationNumber->find(Z)->second->size();
315   size_t targetOscillator = 0;                 << 277   G4int iosc = nosc;
316   G4double ionEnergy = 0.0*eV;                 << 278   ki = photonEnergy0/electron_mass_c2;
317                                                << 279   ki2 = 2*ki+1.0;
318   G4double ek = photonEnergy0/electron_mass_c2 << 280   ki3 = ki*ki;
319   G4double ek2 = 2.*ek+1.0;                    << 281   ki1 = ki3-ki2-1.0;
320   G4double eks = ek*ek;                        << 282   taumin = 1.0/ki2;
321   G4double ek1 = eks-ek2-1.0;                  << 283   a1 = std::log(ki2);
322                                                << 284   a2 = a1+2.0*ki*(1.0+ki)/(ki2*ki2);
323   G4double taumin = 1.0/ek2;                   << 285   //If the incoming photon is above 5 MeV, the quicker approach based on the 
324   //This is meant to fix a possible (rare) flo << 
325   //causing an infinite loop. The maximum of t << 
326   //be represented (i.e. ~ 1. - 1e-16). Fix by << 
327   static G4double taumax = std::nexttoward(1.0 << 
328   if (fVerboseLevel > 3)                       << 
329     G4cout << "G4PenelopeComptonModel: maximum << 
330   //To here.                                   << 
331   G4double a1 = G4Log(ek2);                    << 
332   G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2);  << 
333                                                << 
334   G4double TST = 0;                            << 
335   G4double tau = 0.;                           << 
336                                                << 
337   //If the incoming photon is above 5 MeV, the << 
338   //pure Klein-Nishina formula is used            286   //pure Klein-Nishina formula is used
339   if (photonEnergy0 > 5*MeV)                      287   if (photonEnergy0 > 5*MeV)
340     {                                             288     {
341       do{                                         289       do{
342   do{                                             290   do{
343     if ((a2*G4UniformRand()) < a1)                291     if ((a2*G4UniformRand()) < a1)
344       tau = std::pow(taumin,G4UniformRand());  << 292       {
                                                   >> 293         tau = std::pow(taumin,G4UniformRand());
                                                   >> 294       }
345     else                                          295     else
346       tau = std::sqrt(1.0+G4UniformRand()*(tau << 296       {
                                                   >> 297         tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
                                                   >> 298       }
347     //rejection function                          299     //rejection function
348     TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(e << 300     TST = (1+tau*(ki1+tau*(ki2+tau*ki3)))/(ki3*tau*(1.0+tau*tau));
349   }while (G4UniformRand()> TST);                  301   }while (G4UniformRand()> TST);
350   if (tau > taumax) tau = taumax; //prevent FP << 
351   epsilon=tau;                                    302   epsilon=tau;
352   cosTheta = 1.0 - (1.0-tau)/(ek*tau);         << 303   cosTheta = 1.0 - (1.0-tau)/(ki*tau);
353                                                << 
354   //Target shell electrons                        304   //Target shell electrons
355   TST = fOscManager->GetTotalZ(material)*G4Uni << 305   TST = Z*G4UniformRand();
356   targetOscillator = numberOfOscillators-1; // << 306   iosc = nosc;
357   S=0.0;                                          307   S=0.0;
358   G4bool levelFound = false;                   << 308   for (G4int j=0;j<nosc;j++)
359   for (size_t j=0;j<numberOfOscillators && !le << 
360     {                                             309     {
361       S += (*theTable)[j]->GetOscillatorStreng << 310       occupNb = (G4int) (*(occupationNumber->find(Z)->second))[j];
362       if (S > TST)                             << 311       S = S + occupNb;
363         {                                      << 312       if (S > TST) iosc = j;
364     targetOscillator = j;                      << 313       if (S > TST) break; 
365     levelFound = true;                         << 
366         }                                      << 
367     }                                             314     }
368   //check whether the level is valid           << 315   ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
369   ionEnergy = (*theTable)[targetOscillator]->G << 
370       }while((epsilon*photonEnergy0-photonEner    316       }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
371     }                                             317     }
372   else //photonEnergy0 < 5 MeV                    318   else //photonEnergy0 < 5 MeV
373     {                                             319     {
374       //Incoherent scattering function for the    320       //Incoherent scattering function for theta=PI
375       G4double s0=0.0;                            321       G4double s0=0.0;
376       G4double pzomc=0.0;                      << 322       G4double pzomc=0.0,rni=0.0;
377       G4double rni=0.0;                        << 
378       G4double aux=0.0;                           323       G4double aux=0.0;
379       for (size_t i=0;i<numberOfOscillators;i+ << 324       for (G4int i=0;i<nosc;i++){
380   {                                            << 325   ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
381     ionEnergy = (*theTable)[i]->GetIonisationE << 326   if (photonEnergy0 > ionEnergy)
382     if (photonEnergy0 > ionEnergy)             << 327     {
383       {                                        << 328       G4double aux = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
384         G4double aux2 = photonEnergy0*(photonE << 329       harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
385         hartreeFunc = (*theTable)[i]->GetHartr << 330       occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
386         oscStren = (*theTable)[i]->GetOscillat << 331       pzomc = harFunc*(aux-electron_mass_c2*ionEnergy)/
387         pzomc = hartreeFunc*(aux2-electron_mas << 332          (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
388     (electron_mass_c2*std::sqrt(2.0*aux2+ionEn << 333       if (pzomc > 0) 
389         if (pzomc > 0)                         << 334         {
390     rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+st << 335     rni = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
391                (std::sqrt(0.5)+std::sqrt(2.0)*    336                (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
392         else                                   << 337         }
393     rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::s << 338       else
                                                   >> 339         {
                                                   >> 340     rni = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
394            (std::sqrt(0.5)-std::sqrt(2.0)*pzom    341            (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
395         s0 += oscStren*rni;                    << 342         }
396       }                                        << 343       s0 = s0 + occupNb*rni;
397   }                                            << 344     }
                                                   >> 345       }
                                                   >> 346       
398       //Sampling tau                              347       //Sampling tau
399       G4double cdt1 = 0.;                      << 348       G4double cdt1;
400       do                                          349       do
401   {                                               350   {
402     if ((G4UniformRand()*a2) < a1)                351     if ((G4UniformRand()*a2) < a1)
403       tau = std::pow(taumin,G4UniformRand());  << 352       {
                                                   >> 353         tau = std::pow(taumin,G4UniformRand());
                                                   >> 354       }
404     else                                          355     else
405       tau = std::sqrt(1.0+G4UniformRand()*(tau << 
406     if (tau > taumax) tau = taumax; //prevent  << 
407     cdt1 = (1.0-tau)/(ek*tau);                 << 
408     //Incoherent scattering function           << 
409     S = 0.;                                    << 
410     for (size_t i=0;i<numberOfOscillators;i++) << 
411       {                                           356       {
412         ionEnergy = (*theTable)[i]->GetIonisat << 357         tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
413         if (photonEnergy0 > ionEnergy) //sum o << 358       }
414     {                                          << 359     cdt1 = (1.0-tau)/(ki*tau);
415       aux = photonEnergy0*(photonEnergy0-ionEn << 360     S=0.0;
416       hartreeFunc = (*theTable)[i]->GetHartree << 361     //Incoherent scattering function
417       oscStren = (*theTable)[i]->GetOscillator << 362     for (G4int i=0;i<nosc;i++){
418       pzomc = hartreeFunc*(aux-electron_mass_c << 363       ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
419         (electron_mass_c2*std::sqrt(2.0*aux+io << 364       if (photonEnergy0 > ionEnergy) //sum only on excitable levels
420       if (pzomc > 0)                           << 365         {
421         rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0 << 366     aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
                                                   >> 367     harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
                                                   >> 368     occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
                                                   >> 369     pzomc = harFunc*(aux-electron_mass_c2*ionEnergy)/
                                                   >> 370       (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
                                                   >> 371     if (pzomc > 0) 
                                                   >> 372       {
                                                   >> 373         rn[i] = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
422                (std::sqrt(0.5)+std::sqrt(2.0)*    374                (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
423       else                                     << 375       }
424         rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)- << 376     else
                                                   >> 377       {
                                                   >> 378         rn[i] = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
425            (std::sqrt(0.5)-std::sqrt(2.0)*pzom    379            (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
426       S += oscStren*rn[i];                     << 380       }
427       pac[i] = S;                              << 381     S = S + occupNb*rn[i];
428     }                                          << 382     pac[i] = S;
429         else                                   << 383         }
430     pac[i] = S-1e-6;                           << 384       else
431       }                                        << 385         {
                                                   >> 386     pac[i] = S-(1e-06);
                                                   >> 387         }
                                                   >> 388     }
432     //Rejection function                          389     //Rejection function
433     TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/ << 390     TST = S*(1.0+tau*(ki1+tau*(ki2+tau*ki3)))/(ki3*tau*(1.0+tau*tau));  
434   }while ((G4UniformRand()*s0) > TST);            391   }while ((G4UniformRand()*s0) > TST);
435                                                << 392       //Target electron shell
436       cosTheta = 1.0 - cdt1;                      393       cosTheta = 1.0 - cdt1;
437       G4double fpzmax=0.0,fpz=0.0;                394       G4double fpzmax=0.0,fpz=0.0;
438       G4double A=0.0;                             395       G4double A=0.0;
439       //Target electron shell                  << 
440       do                                          396       do
441   {                                               397   {
442     do                                            398     do
443       {                                           399       {
444         TST = S*G4UniformRand();               << 400         TST =S*G4UniformRand();
445         targetOscillator = numberOfOscillators << 401         iosc=nosc;
446         G4bool levelFound = false;             << 402         for (G4int i=0;i<nosc;i++){
447         for (size_t i=0;i<numberOfOscillators  << 403     if (pac[i]>TST) iosc = i;
                                                   >> 404     if (pac[i]>TST) break; 
                                                   >> 405         }
                                                   >> 406         A = G4UniformRand()*rn[iosc];
                                                   >> 407         harFunc = (*(hartreeFunction->find(Z)->second))[iosc]/fine_structure_const;
                                                   >> 408         occupNb = (G4int) (*(occupationNumber->find(Z)->second))[iosc];
                                                   >> 409         if (A < 0.5) {
                                                   >> 410     pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
                                                   >> 411       (std::sqrt(2.0)*harFunc);
                                                   >> 412         }
                                                   >> 413         else
448     {                                             414     {
449       if (pac[i]>TST)                          << 415       pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/
450         {                                      << 416         (std::sqrt(2.0)*harFunc);
451           targetOscillator = i;                << 
452           levelFound = true;                   << 
453         }                                      << 
454     }                                             417     }
455         A = G4UniformRand()*rn[targetOscillato << 
456         hartreeFunc = (*theTable)[targetOscill << 
457         oscStren = (*theTable)[targetOscillato << 
458         if (A < 0.5)                           << 
459     pzomc = (std::sqrt(0.5)-std::sqrt(0.5-G4Lo << 
460       (std::sqrt(2.0)*hartreeFunc);            << 
461         else                                   << 
462     pzomc = (std::sqrt(0.5-G4Log(2.0-2.0*A))-s << 
463       (std::sqrt(2.0)*hartreeFunc);            << 
464       } while (pzomc < -1);                       418       } while (pzomc < -1);
465                                                << 
466     // F(EP) rejection                            419     // F(EP) rejection
467     G4double XQC = 1.0+tau*(tau-2.0*cosTheta);    420     G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
468     G4double AF = std::sqrt(XQC)*(1.0+tau*(tau    421     G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
469     if (AF > 0)                                << 422     if (AF > 0) {
470       fpzmax = 1.0+AF*0.2;                        423       fpzmax = 1.0+AF*0.2;
                                                   >> 424     }
471     else                                          425     else
472       fpzmax = 1.0-AF*0.2;                     << 426       {
                                                   >> 427         fpzmax = 1.0-AF*0.2;
                                                   >> 428       }
473     fpz = 1.0+AF*std::max(std::min(pzomc,0.2),    429     fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
474   }while ((fpzmax*G4UniformRand())>fpz);          430   }while ((fpzmax*G4UniformRand())>fpz);
475                                                << 431   
476       //Energy of the scattered photon            432       //Energy of the scattered photon
477       G4double T = pzomc*pzomc;                   433       G4double T = pzomc*pzomc;
478       G4double b1 = 1.0-T*tau*tau;                434       G4double b1 = 1.0-T*tau*tau;
479       G4double b2 = 1.0-T*tau*cosTheta;           435       G4double b2 = 1.0-T*tau*cosTheta;
480       if (pzomc > 0.0)                            436       if (pzomc > 0.0)
481   epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2 << 437   {
                                                   >> 438     epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
                                                   >> 439   }
482       else                                        440       else
483   epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2 << 441   {
484     } //energy < 5 MeV                         << 442     epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
                                                   >> 443   }
                                                   >> 444     }
                                                   >> 445   
485                                                   446 
486   //Ok, the kinematics has been calculated.       447   //Ok, the kinematics has been calculated.
487   G4double sinTheta = std::sqrt(1-cosTheta*cos    448   G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
488   G4double phi = twopi * G4UniformRand() ;        449   G4double phi = twopi * G4UniformRand() ;
489   G4double dirx = sinTheta * std::cos(phi);       450   G4double dirx = sinTheta * std::cos(phi);
490   G4double diry = sinTheta * std::sin(phi);       451   G4double diry = sinTheta * std::sin(phi);
491   G4double dirz = cosTheta ;                      452   G4double dirz = cosTheta ;
492                                                   453 
493   // Update G4VParticleChange for the scattere    454   // Update G4VParticleChange for the scattered photon
494   G4ThreeVector photonDirection1(dirx,diry,dir    455   G4ThreeVector photonDirection1(dirx,diry,dirz);
495   photonDirection1.rotateUz(photonDirection0);    456   photonDirection1.rotateUz(photonDirection0);
496   fParticleChange->ProposeMomentumDirection(ph    457   fParticleChange->ProposeMomentumDirection(photonDirection1) ;
497                                                << 
498   G4double photonEnergy1 = epsilon * photonEne    458   G4double photonEnergy1 = epsilon * photonEnergy0;
499                                                   459 
500   if (photonEnergy1 > 0.)                         460   if (photonEnergy1 > 0.)
                                                   >> 461   {
501     fParticleChange->SetProposedKineticEnergy(    462     fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
                                                   >> 463   }
502   else                                            464   else
503   {                                               465   {
504     fParticleChange->SetProposedKineticEnergy(    466     fParticleChange->SetProposedKineticEnergy(0.) ;
505     fParticleChange->ProposeTrackStatus(fStopA    467     fParticleChange->ProposeTrackStatus(fStopAndKill);
506   }                                               468   }
507                                                << 469   
508   //Create scattered electron                     470   //Create scattered electron
509   G4double diffEnergy = photonEnergy0*(1-epsil    471   G4double diffEnergy = photonEnergy0*(1-epsilon);
510   ionEnergy = (*theTable)[targetOscillator]->G << 472   ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
511                                                << 473   G4double Q2 = 
512   G4double Q2 =                                << 
513     photonEnergy0*photonEnergy0+photonEnergy1*    474     photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
514   G4double cosThetaE = 0.; //scattering angle  << 475   G4double cosThetaE; //scattering angle for the electron
515                                                << 
516   if (Q2 > 1.0e-12)                               476   if (Q2 > 1.0e-12)
517     cosThetaE = (photonEnergy0-photonEnergy1*c << 477     {
                                                   >> 478       cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
                                                   >> 479     }
518   else                                            480   else
519     cosThetaE = 1.0;                           << 481     {
                                                   >> 482       cosThetaE = 1.0;
                                                   >> 483     }
520   G4double sinThetaE = std::sqrt(1-cosThetaE*c    484   G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
521                                                   485 
522   //Now, try to handle fluorescence            << 
523   //Notice: merged levels are indicated with Z << 
524   G4int shFlag = (*theTable)[targetOscillator] << 
525   G4int Z = (G4int) (*theTable)[targetOscillat << 
526                                                << 
527   //initialize here, then check photons create    486   //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
528   G4double bindingEnergy = 0.*eV;              << 487   std::vector<G4DynamicParticle*>* photonVector=0;
529   const G4AtomicShell* shell = 0;              << 
530                                                << 
531   //Real level                                 << 
532   if (Z > 0 && shFlag<30)                      << 
533     {                                          << 
534       shell = fTransitionManager->Shell(Z,shFl << 
535       bindingEnergy = shell->BindingEnergy();  << 
536     }                                          << 
537                                                   488 
                                                   >> 489   const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
                                                   >> 490   const G4AtomicShell* shell = transitionManager->Shell(Z,iosc);
                                                   >> 491   G4double bindingEnergy = shell->BindingEnergy();
                                                   >> 492   G4int shellId = shell->ShellId();
538   G4double ionEnergyInPenelopeDatabase = ionEn    493   G4double ionEnergyInPenelopeDatabase = ionEnergy;
539   //protection against energy non-conservation    494   //protection against energy non-conservation
540   ionEnergy = std::max(bindingEnergy,ionEnergy << 495   ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase);  
541                                                   496 
542   //subtract the excitation energy. If not emi    497   //subtract the excitation energy. If not emitted by fluorescence
543   //the ionization energy is deposited as loca    498   //the ionization energy is deposited as local energy deposition
544   G4double eKineticEnergy = diffEnergy - ionEn << 499   G4double eKineticEnergy = diffEnergy - ionEnergy; 
545   G4double localEnergyDeposit = ionEnergy;     << 500   G4double localEnergyDeposit = ionEnergy; 
546   G4double energyInFluorescence = 0.; //testin    501   G4double energyInFluorescence = 0.; //testing purposes only
547   G4double energyInAuger = 0; //testing purpos << 
548                                                   502 
549   if (eKineticEnergy < 0)                      << 503   if (eKineticEnergy < 0) 
550     {                                             504     {
551       //It means that there was some problem/m << 505       //It means that there was some problem/mismatch between the two databases. 
552       //Try to make it work                       506       //Try to make it work
553       //In this case available Energy (diffEne    507       //In this case available Energy (diffEnergy) < ionEnergy
554       //Full residual energy is deposited loca    508       //Full residual energy is deposited locally
555       localEnergyDeposit = diffEnergy;            509       localEnergyDeposit = diffEnergy;
556       eKineticEnergy = 0.0;                       510       eKineticEnergy = 0.0;
557     }                                             511     }
558                                                << 512      
559   //the local energy deposit is what remains:     513   //the local energy deposit is what remains: part of this may be spent for fluorescence.
560   //Notice: shell might be nullptr (invalid!)  << 514   if(DeexcitationFlag() && Z > 5) {
561   //Now, take care of fluorescence, if require << 
562   if (fAtomDeexcitation && shell)              << 
563     {                                          << 
564       G4int index = couple->GetIndex();        << 
565       if (fAtomDeexcitation->CheckDeexcitation << 
566   {                                            << 
567     size_t nBefore = fvect->size();            << 
568     fAtomDeexcitation->GenerateParticles(fvect << 
569     size_t nAfter = fvect->size();             << 
570                                                   515 
571     if (nAfter > nBefore) //actual production  << 516     const G4ProductionCutsTable* theCoupleTable=
572       {                                        << 517       G4ProductionCutsTable::GetProductionCutsTable();
573         for (size_t j=nBefore;j<nAfter;j++) // << 
574     {                                          << 
575       G4double itsEnergy = ((*fvect)[j])->GetK << 
576       if (itsEnergy < localEnergyDeposit) // v << 
577         {                                      << 
578           localEnergyDeposit -= itsEnergy;     << 
579           if (((*fvect)[j])->GetParticleDefini << 
580       energyInFluorescence += itsEnergy;       << 
581           else if (((*fvect)[j])->GetParticleD << 
582              G4Electron::Definition())         << 
583       energyInAuger += itsEnergy;              << 
584         }                                      << 
585       else //invalid secondary: takes more tha << 
586         {                                      << 
587           delete (*fvect)[j];                  << 
588           (*fvect)[j] = nullptr;               << 
589         }                                      << 
590     }                                          << 
591       }                                        << 
592                                                << 
593   }                                            << 
594     }                                          << 
595                                                   518 
596   //Always produce explicitly the electron     << 519     size_t index = couple->GetIndex();
597   G4DynamicParticle* electron = 0;             << 520     G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
                                                   >> 521     G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
                                                   >> 522 
                                                   >> 523     // Generation of fluorescence
                                                   >> 524     // Data in EADL are available only for Z > 5
                                                   >> 525     // Protection to avoid generating photons in the unphysical case of
                                                   >> 526     // shell binding energy > photon energy
                                                   >> 527     if (localEnergyDeposit > cutg || localEnergyDeposit > cute)
                                                   >> 528       { 
                                                   >> 529   G4DynamicParticle* aPhoton;
                                                   >> 530   deexcitationManager.SetCutForSecondaryPhotons(cutg);
                                                   >> 531   deexcitationManager.SetCutForAugerElectrons(cute);
                                                   >> 532       
                                                   >> 533   photonVector = deexcitationManager.GenerateParticles(Z,shellId);
                                                   >> 534   if(photonVector) 
                                                   >> 535     {
                                                   >> 536       size_t nPhotons = photonVector->size();
                                                   >> 537       for (size_t k=0; k<nPhotons; k++)
                                                   >> 538         {
                                                   >> 539     aPhoton = (*photonVector)[k];
                                                   >> 540     if (aPhoton)
                                                   >> 541       {
                                                   >> 542         G4double itsEnergy = aPhoton->GetKineticEnergy();
                                                   >> 543         if (itsEnergy <= localEnergyDeposit)
                                                   >> 544           {
                                                   >> 545       localEnergyDeposit -= itsEnergy;
                                                   >> 546       if (aPhoton->GetDefinition() == G4Gamma::Gamma()) 
                                                   >> 547         energyInFluorescence += itsEnergy;;
                                                   >> 548       fvect->push_back(aPhoton);        
                                                   >> 549           }
                                                   >> 550         else
                                                   >> 551           {
                                                   >> 552       delete aPhoton;
                                                   >> 553       (*photonVector)[k]=0;
                                                   >> 554           }
                                                   >> 555       }
                                                   >> 556         }
                                                   >> 557       delete photonVector;
                                                   >> 558     }
                                                   >> 559       }
                                                   >> 560   }
598                                                   561 
599   G4double xEl = sinThetaE * std::cos(phi+pi); << 562   //Produce explicitely the electron only if its proposed kinetic energy is 
600   G4double yEl = sinThetaE * std::sin(phi+pi); << 563   //above the cut, otherwise add local energy deposition
601   G4double zEl = cosThetaE;                    << 564   G4DynamicParticle* electron = 0;
602   G4ThreeVector eDirection(xEl,yEl,zEl); //ele << 565   //  if (eKineticEnergy > cutE) // VI: may be consistency problem if cut is applied here
603   eDirection.rotateUz(photonDirection0);       << 566   if (eKineticEnergy > 0.0)
604   electron = new G4DynamicParticle (G4Electron << 567     {
605             eDirection,eKineticEnergy) ;       << 568       G4double xEl = sinThetaE * std::cos(phi+pi); 
606   fvect->push_back(electron);                  << 569       G4double yEl = sinThetaE * std::sin(phi+pi);
                                                   >> 570       G4double zEl = cosThetaE;
                                                   >> 571       G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
                                                   >> 572       eDirection.rotateUz(photonDirection0);
                                                   >> 573       electron = new G4DynamicParticle (G4Electron::Electron(),
                                                   >> 574           eDirection,eKineticEnergy) ;
                                                   >> 575       fvect->push_back(electron);
                                                   >> 576     }
                                                   >> 577   else
                                                   >> 578     {
                                                   >> 579       localEnergyDeposit += eKineticEnergy;
                                                   >> 580     }
607                                                   581 
608   if (localEnergyDeposit < 0) //Should not be: << 582   if (localEnergyDeposit < 0)
609     {                                             583     {
610       G4Exception("G4PenelopeComptonModel::Sam << 584       G4cout << "WARNING-" 
611        "em2099",JustWarning,"WARNING: Negative << 585        << "G4PenelopeComptonModel::SampleSecondaries - Negative energy deposit"
                                                   >> 586        << G4endl;
612       localEnergyDeposit=0.;                      587       localEnergyDeposit=0.;
613     }                                             588     }
614   fParticleChange->ProposeLocalEnergyDeposit(l    589   fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
615                                                << 590   
616   G4double electronEnergy = 0.;                   591   G4double electronEnergy = 0.;
617   if (electron)                                << 592   if (verboseLevel > 1)
618     electronEnergy = eKineticEnergy;           << 
619   if (fVerboseLevel > 1)                       << 
620     {                                             593     {
621       G4cout << "-----------------------------    594       G4cout << "-----------------------------------------------------------" << G4endl;
622       G4cout << "Energy balance from G4Penelop    595       G4cout << "Energy balance from G4PenelopeCompton" << G4endl;
623       G4cout << "Incoming photon energy: " <<     596       G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl;
624       G4cout << "-----------------------------    597       G4cout << "-----------------------------------------------------------" << G4endl;
625       G4cout << "Scattered photon: " << photon    598       G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl;
                                                   >> 599       if (electron)
                                                   >> 600   electronEnergy = eKineticEnergy;
626       G4cout << "Scattered electron " << elect    601       G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl;
627       if (energyInFluorescence)                << 602       G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
628   G4cout << "Fluorescence x-rays: " << energyI << 
629       if (energyInAuger)                       << 
630   G4cout << "Auger electrons: " << energyInAug << 
631       G4cout << "Local energy deposit " << loc    603       G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
632       G4cout << "Total final state: " << (phot    604       G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+
633             localEnergyDeposit+energyInAuger)/ << 605             localEnergyDeposit)/keV << 
634   " keV" << G4endl;                               606   " keV" << G4endl;
635       G4cout << "-----------------------------    607       G4cout << "-----------------------------------------------------------" << G4endl;
636     }                                             608     }
637   if (fVerboseLevel > 0)                       << 609   if (verboseLevel > 0)
638     {                                             610     {
639       G4double energyDiff = std::fabs(photonEn    611       G4double energyDiff = std::fabs(photonEnergy1+
640               electronEnergy+energyInFluoresce    612               electronEnergy+energyInFluorescence+
641               localEnergyDeposit+energyInAuger << 613               localEnergyDeposit-photonEnergy0);
642       if (energyDiff > 0.05*keV)                  614       if (energyDiff > 0.05*keV)
643   G4cout << "Warning from G4PenelopeCompton: p << 615   G4cout << "Warning from G4PenelopeCompton: problem with energy conservation: " << 
644     (photonEnergy1+electronEnergy+energyInFluo << 616     (photonEnergy1+electronEnergy+energyInFluorescence+localEnergyDeposit)/keV << 
645     " keV (final) vs. " <<                     << 617     " keV (final) vs. " << 
646     photonEnergy0/keV << " keV (initial)" << G    618     photonEnergy0/keV << " keV (initial)" << G4endl;
647     }                                             619     }
648 }                                                 620 }
649                                                   621 
650 //....oooOO0OOooo........oooOO0OOooo........oo    622 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
651                                                   623 
652 G4double G4PenelopeComptonModel::DifferentialC << 624 void G4PenelopeComptonModel::ReadData()
653                   G4PenelopeOscillator* osc)   << 
654 {                                                 625 {
655   //                                           << 626   char* path = getenv("G4LEDATA");
656   // Penelope model v2008. Single differential << 627   if (!path)
657   // for photon Compton scattering by          << 628     {
658   // electrons in the given atomic oscillator, << 629       G4String excep = "G4PenelopeComptonModel - G4LEDATA environment variable not set!";
659   // scattering photon. This is in units of pi << 630       G4Exception(excep);
660   //                                           << 631     }
661   // D. Brusa et al., Nucl. Instrum. Meth. A 3 << 632   G4String pathString(path);
662   // The parametrization includes the J(p) dis << 633   G4String pathFile = pathString + "/penelope/compton-pen.dat";
663   // that are tabulated from Hartree-Fock calc << 634   std::ifstream file(pathFile);
664   // from F. Biggs et al., At. Data Nucl. Data << 635   
665   //                                           << 636   if (!file.is_open())
666   G4double ionEnergy = osc->GetIonisationEnerg << 
667   G4double harFunc = osc->GetHartreeFactor();  << 
668                                                << 
669   static const G4double k2 = std::sqrt(2.);    << 
670   static const G4double k1 = 1./k2;            << 
671                                                << 
672   if (energy < ionEnergy)                      << 
673     return 0;                                  << 
674                                                << 
675   //energy of the Compton line                 << 
676   G4double cdt1 = 1.0-cosTheta;                << 
677   G4double EOEC = 1.0+(energy/electron_mass_c2 << 
678   G4double ECOE = 1.0/EOEC;                    << 
679                                                << 
680   //Incoherent scattering function (analytical << 
681   G4double aux = energy*(energy-ionEnergy)*cdt << 
682   G4double Pzimax =                            << 
683     (aux - electron_mass_c2*ionEnergy)/(electr << 
684   G4double sia = 0.0;                          << 
685   G4double x = harFunc*Pzimax;                 << 
686   if (x > 0)                                   << 
687     sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x << 
688   else                                         << 
689     sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x));  << 
690                                                << 
691   //1st order correction, integral of Pz times << 
692   //Calculated approximately using a free-elec << 
693   G4double pf = 3.0/(4.0*harFunc);             << 
694   if (std::fabs(Pzimax) < pf)                  << 
695     {                                             637     {
696       G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE* << 638       G4String excep = "G4PenelopeComptonModel - data file " + pathFile + " not found!";
697       G4double p2 = Pzimax*Pzimax;             << 639       G4Exception(excep);
698       G4double dspz = std::sqrt(QCOE2)*        << 
699   (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc     << 
700   *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf));        << 
701       sia += std::max(dspz,-1.0*sia);          << 
702     }                                             640     }
703                                                   641 
704   G4double XKN = EOEC+ECOE-1.0+cosTheta*cosThe << 642   G4int k1,test,test1;
                                                   >> 643   G4double a1,a2;
                                                   >> 644   G4int Z=1,nLevels=0;
705                                                   645 
706   //Differential cross section (per electron,  << 646   if (!ionizationEnergy || !hartreeFunction || !occupationNumber)
707   G4double diffCS = ECOE*ECOE*XKN*sia;         << 647     {
                                                   >> 648       G4String excep = "G4PenelopeComptonModel: problem with reading data from file";
                                                   >> 649       G4Exception(excep);
                                                   >> 650     }
708                                                   651 
709   return diffCS;                               << 652   do{
                                                   >> 653     G4double harOfElectronsBelowThreshold = 0;
                                                   >> 654     G4int nbOfElectronsBelowThreshold = 0;
                                                   >> 655     G4DataVector* occVector = new G4DataVector;
                                                   >> 656     G4DataVector* harVector = new G4DataVector;
                                                   >> 657     G4DataVector* bindingEVector = new G4DataVector;
                                                   >> 658     file >> Z >> nLevels;
                                                   >> 659     for (G4int h=0;h<nLevels;h++)
                                                   >> 660       {
                                                   >> 661   file >> k1 >> a1 >> a2;
                                                   >> 662   //Make explicit unit of measurements for ionisation energy, which is MeV
                                                   >> 663         a1 *= MeV; 
                                                   >> 664   if (a1 > 15*eV)
                                                   >> 665     {
                                                   >> 666       occVector->push_back((G4double) k1);
                                                   >> 667       bindingEVector->push_back(a1);
                                                   >> 668       harVector->push_back(a2);
                                                   >> 669     }
                                                   >> 670   else
                                                   >> 671     {
                                                   >> 672       nbOfElectronsBelowThreshold += k1;
                                                   >> 673       harOfElectronsBelowThreshold += k1*a2;
                                                   >> 674     }
                                                   >> 675       }
                                                   >> 676     //Add the "final" level
                                                   >> 677     if (nbOfElectronsBelowThreshold)
                                                   >> 678       {
                                                   >> 679   occVector->push_back(nbOfElectronsBelowThreshold);
                                                   >> 680   bindingEVector->push_back(0*eV);
                                                   >> 681   G4double averageHartree = 
                                                   >> 682     harOfElectronsBelowThreshold/((G4double) nbOfElectronsBelowThreshold);
                                                   >> 683   harVector->push_back(averageHartree);
                                                   >> 684       }
                                                   >> 685     //Ok, done for element Z
                                                   >> 686     occupationNumber->insert(std::make_pair(Z,occVector));
                                                   >> 687     ionizationEnergy->insert(std::make_pair(Z,bindingEVector));
                                                   >> 688     hartreeFunction->insert(std::make_pair(Z,harVector));
                                                   >> 689     file >> test >> test1; //-1 -1 close the data for each Z
                                                   >> 690     if (test > 0) {
                                                   >> 691       G4String excep = "G4PenelopeComptonModel - data file corrupted!";
                                                   >> 692       G4Exception(excep);
                                                   >> 693     }
                                                   >> 694   }while (test != -2); //the very last Z is closed with -2 instead of -1
                                                   >> 695   file.close();
                                                   >> 696   if (verboseLevel > 2)
                                                   >> 697     {
                                                   >> 698       G4cout << "Data from G4PenelopeComptonModel read " << G4endl;
                                                   >> 699     }
710 }                                                 700 }
711                                                   701 
712 //....oooOO0OOooo........oooOO0OOooo........oo    702 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
713                                                   703 
714 G4double G4PenelopeComptonModel::OscillatorTot << 704 G4double G4PenelopeComptonModel::DifferentialCrossSection(G4double cosTheta)
715 {                                                 705 {
716   //Total cross section (integrated) for the g << 706   //
717   //pi*classic_electr_radius^2                 << 707   // Penelope model for the Compton scattering differential cross section 
718                                                << 708   // dSigma/dOmega.
719   //Integrate differential cross section for e << 709   // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
720   G4double stre = osc->GetOscillatorStrength() << 710   // The parametrization includes the J(p) distribution profiles for the atomic shells, 
721                                                << 711   // that are tabulated from Hartree-Fock calculations 
722   // here one uses the  using the 20-point     << 712   // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
723   // Gauss quadrature method with an adaptive  << 713   //
724   const G4int npoints=10;                      << 714   const G4double k2 = std::sqrt(2.0);
725   const G4int ncallsmax=20000;                 << 715   const G4double k1 = std::sqrt(0.5);
726   const G4int nst=256;                         << 716   const G4double k12 = 0.5;
727   static const G4double Abscissas[10] = {7.652 << 717   G4double cdt1 = 1.0-cosTheta;
728           5.1086700195082710e-01,6.36053680726 << 718   G4double energy = energyForIntegration;
729           8.3911697182221882e-01,9.12234428251 << 719   G4int Z = ZForIntegration;
730           9.9312859918509492e-01};             << 720   //energy of Compton line;
731   static const G4double Weights[10] = {1.52753 << 721   G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1; 
732         1.3168863844917663e-01,1.1819453196151 << 722   G4double ECOE = 1.0/EOEC;
733         8.3276741576704749e-02,6.2672048334109 << 723   //Incoherent scattering function (analytical profile)
734         1.7614007139152118e-02};               << 724   G4double sia = 0.0;
735                                                << 725   G4int nosc = occupationNumber->find(Z)->second->size();
736   G4double MaxError = 1e-5;                    << 726   for (G4int i=0;i<nosc;i++){
737   //Error control                              << 727     G4double ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
738   G4double Ctol = std::min(std::max(MaxError,1 << 728     //Sum only of those shells for which E>Eion
739   G4double Ptol = 0.01*Ctol;                   << 729     if (energy > ionEnergy)
740   G4double Err=1e35;                           << 
741                                                << 
742   //Gauss integration from -1 to 1             << 
743   G4double LowPoint = -1.0;                    << 
744   G4double HighPoint = 1.0;                    << 
745                                                << 
746   G4double h=HighPoint-LowPoint;               << 
747   G4double sumga=0.0;                          << 
748   G4double a=0.5*(HighPoint-LowPoint);         << 
749   G4double b=0.5*(HighPoint+LowPoint);         << 
750   G4double c=a*Abscissas[0];                   << 
751   G4double d= Weights[0]*                      << 
752     (DifferentialCrossSection(b+c,energy,osc)+ << 
753   for (G4int i=2;i<=npoints;i++)               << 
754     {                                          << 
755       c=a*Abscissas[i-1];                      << 
756       d += Weights[i-1]*                       << 
757   (DifferentialCrossSection(b+c,energy,osc)+Di << 
758     }                                          << 
759   G4int icall = 2*npoints;                     << 
760   G4int LH=1;                                  << 
761   G4double S[nst],x[nst],sn[nst],xrn[nst];     << 
762   S[0]=d*a;                                    << 
763   x[0]=LowPoint;                               << 
764                                                << 
765   G4bool loopAgain = true;                     << 
766                                                << 
767   //Adaptive bipartition scheme                << 
768   do{                                          << 
769     G4double h0=h;                             << 
770     h=0.5*h; //bipartition                     << 
771     G4double sumr=0;                           << 
772     G4int LHN=0;                               << 
773     G4double si,xa,xb,xc;                      << 
774     for (G4int i=1;i<=LH;i++){                 << 
775       si=S[i-1];                               << 
776       xa=x[i-1];                               << 
777       xb=xa+h;                                 << 
778       xc=xa+h0;                                << 
779       a=0.5*(xb-xa);                           << 
780       b=0.5*(xb+xa);                           << 
781       c=a*Abscissas[0];                        << 
782       G4double dLocal = Weights[0]*            << 
783   (DifferentialCrossSection(b+c,energy,osc)+Di << 
784                                                << 
785       for (G4int j=1;j<npoints;j++)            << 
786   {                                            << 
787     c=a*Abscissas[j];                          << 
788     dLocal += Weights[j]*                      << 
789       (DifferentialCrossSection(b+c,energy,osc << 
790   }                                            << 
791       G4double s1=dLocal*a;                    << 
792       a=0.5*(xc-xb);                           << 
793       b=0.5*(xc+xb);                           << 
794       c=a*Abscissas[0];                        << 
795       dLocal=Weights[0]*                       << 
796   (DifferentialCrossSection(b+c,energy,osc)+Di << 
797                                                << 
798       for (G4int j=1;j<npoints;j++)            << 
799   {                                            << 
800     c=a*Abscissas[j];                          << 
801     dLocal += Weights[j]*                      << 
802       (DifferentialCrossSection(b+c,energy,osc << 
803   }                                            << 
804       G4double s2=dLocal*a;                    << 
805       icall=icall+4*npoints;                   << 
806       G4double s12=s1+s2;                      << 
807       if (std::abs(s12-si)<std::max(Ptol*std:: << 
808   sumga += s12;                                << 
809       else                                     << 
810   {                                            << 
811     sumr += s12;                               << 
812     LHN += 2;                                  << 
813     sn[LHN-1]=s2;                              << 
814     xrn[LHN-1]=xb;                             << 
815     sn[LHN-2]=s1;                              << 
816     xrn[LHN-2]=xa;                             << 
817   }                                            << 
818                                                << 
819       if (icall>ncallsmax || LHN>nst)          << 
820   {                                            << 
821     G4cout << "G4PenelopeComptonModel: " << G4 << 
822     G4cout << "LowPoint: " << LowPoint << ", H << 
823     G4cout << "Tolerance: " << MaxError << G4e << 
824     G4cout << "Calls: " << icall << ", Integra << 
825     G4cout << "Number of open subintervals: "  << 
826     G4cout << "WARNING: the required accuracy  << 
827     loopAgain = false;                         << 
828   }                                            << 
829     }                                          << 
830     Err=std::abs(sumr)/std::max(std::abs(sumr+ << 
831     if (Err < Ctol || LHN == 0)                << 
832       loopAgain = false; //end of cycle        << 
833     LH=LHN;                                    << 
834     for (G4int i=0;i<LH;i++)                   << 
835       {                                           730       {
836   S[i]=sn[i];                                  << 731         G4double aux = energy * (energy-ionEnergy)*cdt1;
837   x[i]=xrn[i];                                 << 732   G4double Pzimax = 
                                                   >> 733     (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
                                                   >> 734   G4double harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
                                                   >> 735   G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
                                                   >> 736   G4double x = harFunc*Pzimax;
                                                   >> 737   G4double siap = 0;
                                                   >> 738   if (x > 0) 
                                                   >> 739     {
                                                   >> 740       siap = 1.0-0.5*std::exp(k12-(k1+k2*x)*(k1+k2*x));
                                                   >> 741     }
                                                   >> 742   else
                                                   >> 743     {
                                                   >> 744       siap = 0.5*std::exp(k12-(k1-k2*x)*(k1-k2*x));
                                                   >> 745     }
                                                   >> 746   sia = sia + occupNb*siap; //sum of all contributions;
838       }                                           747       }
839   }while(Ctol < 1.0 && loopAgain);             << 748   }
840                                                << 749   G4double XKN = EOEC+ECOE-1+cosTheta*cosTheta;
841   G4double xs = stre*sumga;                    << 750   G4double diffCS = pi*classic_electr_radius*classic_electr_radius*ECOE*ECOE*XKN*sia;
842                                                << 751   return diffCS;
843   return xs;                                   << 
844 }                                                 752 }
845                                                   753 
846 //....oooOO0OOooo........oooOO0OOooo........oo    754 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
847                                                   755 
848 G4double G4PenelopeComptonModel::KleinNishinaC << 756 void G4PenelopeComptonModel::ActivateAuger(G4bool augerbool)
849                const G4Material* material)     << 
850 {                                                 757 {
851   // use Klein-Nishina formula                 << 758   if (!DeexcitationFlag() && augerbool)
852   // total cross section in units of pi*classi << 759     {
853   G4double cs = 0;                             << 760       G4cout << "WARNING - G4PenelopeComptonModel" << G4endl;
854                                                << 761       G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
855   G4double ek =energy/electron_mass_c2;        << 762       G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
856   G4double eks = ek*ek;                        << 763     }
857   G4double ek2 = 1.0+ek+ek;                    << 764   deexcitationManager.ActivateAugerElectronProduction(augerbool);
858   G4double ek1 = eks-ek2-1.0;                  << 765   if (verboseLevel > 1)
859                                                << 766     G4cout << "Auger production set to " << augerbool << G4endl;
860   G4double t0 = 1.0/ek2;                       << 
861   G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*G4Lo << 
862                                                << 
863   G4PenelopeOscillatorTable* theTable = fOscMa << 
864                                                << 
865   for (size_t i=0;i<theTable->size();i++)      << 
866     {                                          << 
867       G4PenelopeOscillator* theOsc = (*theTabl << 
868       G4double ionEnergy = theOsc->GetIonisati << 
869       G4double tau=(energy-ionEnergy)/energy;  << 
870       if (tau > t0)                            << 
871   {                                            << 
872     G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1 << 
873     G4double stre = theOsc->GetOscillatorStren << 
874                                                << 
875     cs += stre*(csu-csl);                      << 
876   }                                            << 
877     }                                          << 
878   cs /= (ek*eks);                              << 
879                                                << 
880   return cs;                                   << 
881                                                << 
882 }                                                 767 }
883                                                   768 
884 //....oooOO0OOooo........oooOO0OOooo........oo << 769 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
885                                                << 
886 void G4PenelopeComptonModel::SetParticle(const << 
887 {                                              << 
888   if(!fParticle) {                             << 
889     fParticle = p;                             << 
890   }                                            << 
891 }                                              << 
892                                                   770