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 10.7.p2)


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