Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopeGammaConversionModel.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/G4PenelopeGammaConversionModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4PenelopeGammaConversionModel.cc (Version 10.0.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 // $Id: G4PenelopeGammaConversionModel.cc 79186 2014-02-20 09:20:02Z gcosmo $
 26 //                                                 27 //
 27 // Author: Luciano Pandola                         28 // Author: Luciano Pandola
 28 //                                                 29 //
 29 // History:                                        30 // History:
 30 // --------                                        31 // --------
 31 // 13 Jan 2010   L Pandola    First implementa     32 // 13 Jan 2010   L Pandola    First implementation (updated to Penelope08)
 32 // 24 May 2011   L Pandola    Renamed (make v2     33 // 24 May 2011   L Pandola    Renamed (make v2008 as default Penelope)
 33 // 18 Sep 2013   L Pandola    Migration to MT  <<  34 // 18 Sep 2013   L Pandola    Migration to MT paradigm. Only master model deals with 
 34 //                             data and create     35 //                             data and creates shared tables
 35 //                                                 36 //
 36                                                    37 
 37 #include "G4PenelopeGammaConversionModel.hh"       38 #include "G4PenelopeGammaConversionModel.hh"
 38 #include "G4PhysicalConstants.hh"                  39 #include "G4PhysicalConstants.hh"
 39 #include "G4SystemOfUnits.hh"                      40 #include "G4SystemOfUnits.hh"
 40 #include "G4ParticleDefinition.hh"                 41 #include "G4ParticleDefinition.hh"
 41 #include "G4MaterialCutsCouple.hh"                 42 #include "G4MaterialCutsCouple.hh"
 42 #include "G4ProductionCutsTable.hh"                43 #include "G4ProductionCutsTable.hh"
 43 #include "G4DynamicParticle.hh"                    44 #include "G4DynamicParticle.hh"
 44 #include "G4Element.hh"                            45 #include "G4Element.hh"
 45 #include "G4Gamma.hh"                              46 #include "G4Gamma.hh"
 46 #include "G4Electron.hh"                           47 #include "G4Electron.hh"
 47 #include "G4Positron.hh"                           48 #include "G4Positron.hh"
 48 #include "G4PhysicsFreeVector.hh"                  49 #include "G4PhysicsFreeVector.hh"
 49 #include "G4MaterialCutsCouple.hh"                 50 #include "G4MaterialCutsCouple.hh"
 50 #include "G4AutoLock.hh"                           51 #include "G4AutoLock.hh"
 51 #include "G4Exp.hh"                            << 
 52                                                << 
 53 //....oooOO0OOooo........oooOO0OOooo........oo << 
 54 const G4int G4PenelopeGammaConversionModel::fM << 
 55 G4PhysicsFreeVector* G4PenelopeGammaConversion << 
 56 G4double G4PenelopeGammaConversionModel::fAtom << 
 57                      1.2281e+02,7.3167e+01,6.9 << 
 58                      6.4696e+01,6.1228e+01,5.7 << 
 59                      5.0787e+01,4.7851e+01,4.6 << 
 60                      4.4503e+01,4.3815e+01,4.3 << 
 61                      4.1586e+01,4.0953e+01,4.0 << 
 62                      3.9756e+01,3.9144e+01,3.8 << 
 63                      3.7174e+01,3.6663e+01,3.5 << 
 64                      3.4688e+01,3.4197e+01,3.3 << 
 65                      3.3068e+01,3.2740e+01,3.2 << 
 66                      3.1884e+01,3.1622e+01,3.1 << 
 67                      3.0950e+01,3.0758e+01,3.0 << 
 68                      3.0097e+01,2.9832e+01,2.9 << 
 69                      2.9247e+01,2.9085e+01,2.8 << 
 70                      2.8580e+01,2.8442e+01,2.8 << 
 71                      2.7973e+01,2.7819e+01,2.7 << 
 72                      2.7285e+01,2.7093e+01,2.6 << 
 73                      2.6516e+01,2.6304e+01,2.6 << 
 74                      2.5730e+01,2.5577e+01,2.5 << 
 75                      2.5100e+01,2.4941e+01,2.4 << 
 76                      2.4506e+01,2.4391e+01,2.4 << 
 77                      2.4039e+01,2.3922e+01,2.3 << 
 78                      2.3621e+01,2.3523e+01,2.3 << 
 79                      2.3238e+01,2.3139e+01,2.3 << 
 80                      2.2833e+01,2.2694e+01,2.2 << 
 81                      2.2446e+01,2.2358e+01,2.2 << 
 82                                                    52 
 83 //....oooOO0OOooo........oooOO0OOooo........oo     53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 84                                                    54 
 85 G4PenelopeGammaConversionModel::G4PenelopeGamm     55 G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition* part,
 86                      const G4String& nam)          56                      const G4String& nam)
 87   :G4VEmModel(nam),fParticleChange(nullptr),fP <<  57   :G4VEmModel(nam),fParticleChange(0),fParticle(0),
 88    fEffectiveCharge(nullptr),fMaterialInvScree <<  58    logAtomicCrossSection(0),
 89    fScreeningFunction(nullptr),fIsInitialised( <<  59    fEffectiveCharge(0),fMaterialInvScreeningRadius(0),
 90 {                                              <<  60    fScreeningFunction(0),isInitialised(false),fLocalTable(false)
                                                   >>  61 
                                                   >>  62 {  
 91   fIntrinsicLowEnergyLimit = 2.0*electron_mass     63   fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
 92   fIntrinsicHighEnergyLimit = 100.0*GeV;           64   fIntrinsicHighEnergyLimit = 100.0*GeV;
 93   fSmallEnergy = 1.1*MeV;                          65   fSmallEnergy = 1.1*MeV;
                                                   >>  66   
                                                   >>  67   InitializeScreeningRadii();
 94                                                    68 
 95   if (part)                                        69   if (part)
 96     SetParticle(part);                             70     SetParticle(part);
 97                                                    71 
 98   //  SetLowEnergyLimit(fIntrinsicLowEnergyLim     72   //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
 99   SetHighEnergyLimit(fIntrinsicHighEnergyLimit     73   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
100   //                                               74   //
101   fVerboseLevel= 0;                            <<  75   verboseLevel= 0;
102   // Verbosity scale:                              76   // Verbosity scale:
103   // 0 = nothing                               <<  77   // 0 = nothing 
104   // 1 = warning for energy non-conservation   <<  78   // 1 = warning for energy non-conservation 
105   // 2 = details of energy budget                  79   // 2 = details of energy budget
106   // 3 = calculation of cross sections, file o     80   // 3 = calculation of cross sections, file openings, sampling of atoms
107   // 4 = entering in methods                       81   // 4 = entering in methods
108 }                                                  82 }
109                                                    83 
110 //....oooOO0OOooo........oooOO0OOooo........oo     84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111                                                    85 
112 G4PenelopeGammaConversionModel::~G4PenelopeGam     86 G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel()
113 {                                                  87 {
114   //Delete shared tables, they exist only in t     88   //Delete shared tables, they exist only in the master model
115   if (IsMaster() || fLocalTable)                   89   if (IsMaster() || fLocalTable)
116     {                                          <<  90     {      
117       for(G4int i=0; i<=fMaxZ; ++i)            <<  91       if (logAtomicCrossSection)
118   {                                                92   {
119     if(fLogAtomicCrossSection[i]) {            <<  93     /*
120       delete fLogAtomicCrossSection[i];        <<  94       std::map <G4int,G4PhysicsFreeVector*>::iterator i;
121       fLogAtomicCrossSection[i] = nullptr;     <<  95       for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
122     }                                          <<  96       if (i->second) delete i->second;
                                                   >>  97     */
                                                   >>  98     delete logAtomicCrossSection;
123   }                                                99   }
124       if (fEffectiveCharge)                       100       if (fEffectiveCharge)
125   delete fEffectiveCharge;                        101   delete fEffectiveCharge;
126       if (fMaterialInvScreeningRadius)            102       if (fMaterialInvScreeningRadius)
127   delete fMaterialInvScreeningRadius;             103   delete fMaterialInvScreeningRadius;
128       if (fScreeningFunction)                     104       if (fScreeningFunction)
129   delete fScreeningFunction;                      105   delete fScreeningFunction;
130     }                                             106     }
131 }                                                 107 }
132                                                   108 
                                                   >> 109 
133 //....oooOO0OOooo........oooOO0OOooo........oo    110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134                                                   111 
135 void G4PenelopeGammaConversionModel::Initialis    112 void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition* part,
136             const G4DataVector&)                  113             const G4DataVector&)
137 {                                                 114 {
138   if (fVerboseLevel > 3)                       << 115   if (verboseLevel > 3)
139     G4cout << "Calling  G4PenelopeGammaConvers    116     G4cout << "Calling  G4PenelopeGammaConversionModel::Initialise()" << G4endl;
140                                                   117 
141   SetParticle(part);                              118   SetParticle(part);
142                                                   119 
143   //Only the master model creates/fills/destro    120   //Only the master model creates/fills/destroys the tables
144   if (IsMaster() && part == fParticle)            121   if (IsMaster() && part == fParticle)
145     {                                             122     {
                                                   >> 123       // logAtomicCrossSection is created only once, since it is never cleared
                                                   >> 124       if (!logAtomicCrossSection)
                                                   >> 125   logAtomicCrossSection =  new std::map<G4int,G4PhysicsFreeVector*>;
                                                   >> 126 
146       //delete old material data...               127       //delete old material data...
147       if (fEffectiveCharge)                       128       if (fEffectiveCharge)
148   {                                               129   {
149     delete fEffectiveCharge;                      130     delete fEffectiveCharge;
150     fEffectiveCharge = nullptr;                << 131     fEffectiveCharge = 0;
151   }                                               132   }
152       if (fMaterialInvScreeningRadius)            133       if (fMaterialInvScreeningRadius)
153   {                                               134   {
154     delete fMaterialInvScreeningRadius;           135     delete fMaterialInvScreeningRadius;
155     fMaterialInvScreeningRadius = nullptr;     << 136     fMaterialInvScreeningRadius = 0;
156   }                                               137   }
157       if (fScreeningFunction)                     138       if (fScreeningFunction)
158   {                                               139   {
159     delete fScreeningFunction;                    140     delete fScreeningFunction;
160     fScreeningFunction = nullptr;              << 141     fScreeningFunction = 0;
161   }                                            << 142   }      
162       //and create new ones                       143       //and create new ones
163       fEffectiveCharge = new std::map<const G4    144       fEffectiveCharge = new std::map<const G4Material*,G4double>;
164       fMaterialInvScreeningRadius = new std::m    145       fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
165       fScreeningFunction = new std::map<const     146       fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
166                                                << 147       
167       G4ProductionCutsTable* theCoupleTable =  << 148       G4ProductionCutsTable* theCoupleTable = 
168   G4ProductionCutsTable::GetProductionCutsTabl    149   G4ProductionCutsTable::GetProductionCutsTable();
169                                                   150 
170       for (G4int i=0;i<(G4int)theCoupleTable-> << 151       for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
171   {                                               152   {
172     const G4Material* material =               << 153     const G4Material* material = 
173       theCoupleTable->GetMaterialCutsCouple(i)    154       theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
174     const G4ElementVector* theElementVector =     155     const G4ElementVector* theElementVector = material->GetElementVector();
175                                                << 156     
176     for (std::size_t j=0;j<material->GetNumber << 157     for (size_t j=0;j<material->GetNumberOfElements();j++)
177       {                                           158       {
178         G4int iZ = theElementVector->at(j)->Ge << 159         G4int iZ = (G4int) theElementVector->at(j)->GetZ();
179         //read data files only in the master      160         //read data files only in the master
180         if (iZ <= fMaxZ &&  !fLogAtomicCrossSe << 161         if (!logAtomicCrossSection->count(iZ))
181     ReadDataFile(iZ);                             162     ReadDataFile(iZ);
182       }                                           163       }
183                                                   164 
184     //check if material data are available        165     //check if material data are available
185     if (!fEffectiveCharge->count(material))       166     if (!fEffectiveCharge->count(material))
186       InitializeScreeningFunctions(material);  << 167       InitializeScreeningFunctions(material); 
187   }                                               168   }
188       if (fVerboseLevel > 0) {                 << 169 
                                                   >> 170 
                                                   >> 171       if (verboseLevel > 0) { 
189   G4cout << "Penelope Gamma Conversion model v    172   G4cout << "Penelope Gamma Conversion model v2008 is initialized " << G4endl
190          << "Energy range: "                      173          << "Energy range: "
191          << LowEnergyLimit() / MeV << " MeV -     174          << LowEnergyLimit() / MeV << " MeV - "
192          << HighEnergyLimit() / GeV << " GeV"     175          << HighEnergyLimit() / GeV << " GeV"
193          << G4endl;                               176          << G4endl;
194       }                                           177       }
                                                   >> 178 
195     }                                             179     }
196   if(fIsInitialised) return;                   << 180 
                                                   >> 181 
                                                   >> 182   if(isInitialised) return;
197   fParticleChange = GetParticleChangeForGamma(    183   fParticleChange = GetParticleChangeForGamma();
198   fIsInitialised = true;                       << 184   isInitialised = true;
199 }                                                 185 }
200                                                   186 
201 //....oooOO0OOooo........oooOO0OOooo........oo    187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202                                                   188 
203 void G4PenelopeGammaConversionModel::Initialis    189 void G4PenelopeGammaConversionModel::InitialiseLocal(const G4ParticleDefinition* part,
204                  G4VEmModel *masterModel)         190                  G4VEmModel *masterModel)
205 {                                                 191 {
206   if (fVerboseLevel > 3)                       << 192   if (verboseLevel > 3)
207     G4cout << "Calling  G4PenelopeGammaConvers    193     G4cout << "Calling  G4PenelopeGammaConversionModel::InitialiseLocal()" << G4endl;
                                                   >> 194 
208   //                                              195   //
209   //Check that particle matches: one might hav << 196   //Check that particle matches: one might have multiple master models (e.g. 
210   //for e+ and e-).                               197   //for e+ and e-).
211   //                                              198   //
212   if (part == fParticle)                          199   if (part == fParticle)
213     {                                             200     {
214       //Get the const table pointers from the     201       //Get the const table pointers from the master to the workers
215       const G4PenelopeGammaConversionModel* th << 202       const G4PenelopeGammaConversionModel* theModel = 
216   static_cast<G4PenelopeGammaConversionModel*>    203   static_cast<G4PenelopeGammaConversionModel*> (masterModel);
217                                                << 204       
218       //Copy pointers to the data tables          205       //Copy pointers to the data tables
219       fEffectiveCharge = theModel->fEffectiveC    206       fEffectiveCharge = theModel->fEffectiveCharge;
220       fMaterialInvScreeningRadius = theModel->    207       fMaterialInvScreeningRadius = theModel->fMaterialInvScreeningRadius;
221       fScreeningFunction = theModel->fScreenin << 208       fScreeningFunction = theModel->fScreeningFunction;
222       for(G4int i=0; i<=fMaxZ; ++i)            << 209       logAtomicCrossSection = theModel->logAtomicCrossSection;
223   fLogAtomicCrossSection[i] = theModel->fLogAt << 210       
224                                                << 
225       //Same verbosity for all workers, as the    211       //Same verbosity for all workers, as the master
226       fVerboseLevel = theModel->fVerboseLevel; << 212       verboseLevel = theModel->verboseLevel;
227     }                                          << 213     }      
228                                                   214 
229   return;                                         215   return;
230 }                                                 216 }
231                                                   217 
232 //....oooOO0OOooo........oooOO0OOooo........oo    218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 namespace { G4Mutex  PenelopeGammaConversionMo    219 namespace { G4Mutex  PenelopeGammaConversionModelMutex = G4MUTEX_INITIALIZER; }
234                                                   220 
235 G4double G4PenelopeGammaConversionModel::Compu    221 G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom(
236                     const G4ParticleDefinition    222                     const G4ParticleDefinition*,
237                     G4double energy,              223                     G4double energy,
238                     G4double Z, G4double,         224                     G4double Z, G4double,
239                     G4double, G4double)           225                     G4double, G4double)
240 {                                                 226 {
241   //                                              227   //
242   // Penelope model v2008.                        228   // Penelope model v2008.
243   // Cross section (including triplet producti << 229   // Cross section (including triplet production) read from database and managed 
244   // through the G4CrossSectionHandler utility    230   // through the G4CrossSectionHandler utility. Cross section data are from
245   // M.J. Berger and J.H. Hubbel (XCOM), Repor    231   // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
246   //                                              232   //
247                                                << 233   
248   if (energy < fIntrinsicLowEnergyLimit)          234   if (energy < fIntrinsicLowEnergyLimit)
249     return 0;                                     235     return 0;
250                                                   236 
251   G4int iZ = G4int(Z);                         << 237   G4int iZ = (G4int) Z;
252                                                << 238  
253   if (!fLogAtomicCrossSection[iZ])             << 239   //Either Initialize() was not called, or we are in a slave and InitializeLocal() was 
                                                   >> 240   //not invoked
                                                   >> 241   if (!logAtomicCrossSection)
                                                   >> 242     {
                                                   >> 243       //create a **thread-local** version of the table. Used only for G4EmCalculator and 
                                                   >> 244       //Unit Tests
                                                   >> 245       fLocalTable = true;
                                                   >> 246       logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
                                                   >> 247     }
                                                   >> 248   //now it should be ok
                                                   >> 249   if (!logAtomicCrossSection->count(iZ))
254      {                                            250      {
255        //If we are here, it means that Initial << 251        //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 
256        //not filled up. This can happen in a U    252        //not filled up. This can happen in a UnitTest or via G4EmCalculator
257        if (fVerboseLevel > 0)                  << 253        if (verboseLevel > 0)
258    {                                              254    {
259      //Issue a G4Exception (warning) only in v    255      //Issue a G4Exception (warning) only in verbose mode
260      G4ExceptionDescription ed;                   256      G4ExceptionDescription ed;
261      ed << "Unable to retrieve the cross secti    257      ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
262      ed << "This can happen only in Unit Tests << 258      ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;   
263      G4Exception("G4PenelopeGammaConversionMod    259      G4Exception("G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom()",
264            "em2018",JustWarning,ed);              260            "em2018",JustWarning,ed);
265    }                                              261    }
266        //protect file reading via autolock        262        //protect file reading via autolock
267        G4AutoLock lock(&PenelopeGammaConversio    263        G4AutoLock lock(&PenelopeGammaConversionModelMutex);
268        ReadDataFile(iZ);                       << 264        ReadDataFile(iZ);            
269        lock.unlock();                             265        lock.unlock();
270        fLocalTable = true;                     << 
271      }                                            266      }
                                                   >> 267 
272   G4double cs = 0;                                268   G4double cs = 0;
273   G4double logene = G4Log(energy);             << 269   G4double logene = std::log(energy);
274   G4PhysicsFreeVector* theVec = fLogAtomicCros << 270 
                                                   >> 271   G4PhysicsFreeVector* theVec = logAtomicCrossSection->find(iZ)->second;
                                                   >> 272 
275   G4double logXS = theVec->Value(logene);         273   G4double logXS = theVec->Value(logene);
276   cs = G4Exp(logXS);                           << 274   cs = std::exp(logXS);
277                                                   275 
278   if (fVerboseLevel > 2)                       << 276   if (verboseLevel > 2)
279     G4cout << "Gamma conversion cross section  << 277     G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z << 
280       " = " << cs/barn << " barn" << G4endl;      278       " = " << cs/barn << " barn" << G4endl;
281   return cs;                                      279   return cs;
282 }                                                 280 }
283                                                   281 
284 //....oooOO0OOooo........oooOO0OOooo........oo    282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
285                                                   283 
286 void                                           << 284 void 
287 G4PenelopeGammaConversionModel::SampleSecondar    285 G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
288               const G4MaterialCutsCouple* coup    286               const G4MaterialCutsCouple* couple,
289               const G4DynamicParticle* aDynami    287               const G4DynamicParticle* aDynamicGamma,
290               G4double,                           288               G4double,
291               G4double)                           289               G4double)
292 {                                                 290 {
293   //                                              291   //
294   // Penelope model v2008.                        292   // Penelope model v2008.
295   // Final state is sampled according to the B    293   // Final state is sampled according to the Bethe-Heitler model with Coulomb
296   // corrections, according to the semi-empiri    294   // corrections, according to the semi-empirical model of
297   //  J. Baro' et al., Radiat. Phys. Chem. 44     295   //  J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
298   //                                              296   //
299   // The model uses the high energy Coulomb co << 297   // The model uses the high energy Coulomb correction from 
300   //  H. Davies et al., Phys. Rev. 93 (1954) 7    298   //  H. Davies et al., Phys. Rev. 93 (1954) 788
301   // and atomic screening radii tabulated from << 299   // and atomic screening radii tabulated from 
302   //  J.H. Hubbel et al., J. Phys. Chem. Ref.     300   //  J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
303   // for Z= 1 to 92.                           << 301   // for Z= 1 to 92. 
304   //                                              302   //
305   if (fVerboseLevel > 3)                       << 303   if (verboseLevel > 3)
306     G4cout << "Calling SamplingSecondaries() o    304     G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl;
307                                                   305 
308   G4double photonEnergy = aDynamicGamma->GetKi    306   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
309                                                   307 
310   // Always kill primary                          308   // Always kill primary
311   fParticleChange->ProposeTrackStatus(fStopAnd    309   fParticleChange->ProposeTrackStatus(fStopAndKill);
312   fParticleChange->SetProposedKineticEnergy(0.    310   fParticleChange->SetProposedKineticEnergy(0.);
313                                                   311 
314   if (photonEnergy <= fIntrinsicLowEnergyLimit    312   if (photonEnergy <= fIntrinsicLowEnergyLimit)
315     {                                             313     {
316       fParticleChange->ProposeLocalEnergyDepos    314       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
317       return ;                                    315       return ;
318     }                                             316     }
319                                                   317 
320   G4ParticleMomentum photonDirection = aDynami    318   G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
321   const G4Material* mat = couple->GetMaterial(    319   const G4Material* mat = couple->GetMaterial();
322                                                   320 
323   //Either Initialize() was not called, or we  << 321   //Either Initialize() was not called, or we are in a slave and InitializeLocal() was 
324   //not invoked                                   322   //not invoked
325   if (!fEffectiveCharge)                          323   if (!fEffectiveCharge)
326     {                                             324     {
327       //create a **thread-local** version of t << 325       //create a **thread-local** version of the table. Used only for G4EmCalculator and 
328       //Unit Tests                                326       //Unit Tests
329       fLocalTable = true;                      << 327       fLocalTable = true;    
330       fEffectiveCharge = new std::map<const G4    328       fEffectiveCharge = new std::map<const G4Material*,G4double>;
331       fMaterialInvScreeningRadius = new std::m    329       fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
332       fScreeningFunction = new std::map<const     330       fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
333     }                                             331     }
334                                                   332 
335   if (!fEffectiveCharge->count(mat))              333   if (!fEffectiveCharge->count(mat))
336     {                                             334     {
337       //If we are here, it means that Initiali << 335       //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 
338       //not filled up. This can happen in a Un    336       //not filled up. This can happen in a UnitTest or via G4EmCalculator
339       if (fVerboseLevel > 0)                   << 337       if (verboseLevel > 0)
340   {                                               338   {
341     //Issue a G4Exception (warning) only in ve    339     //Issue a G4Exception (warning) only in verbose mode
342     G4ExceptionDescription ed;                    340     G4ExceptionDescription ed;
343     ed << "Unable to allocate the EffectiveCha << 341     ed << "Unable to allocate the EffectiveCharge data for " << 
344       mat->GetName() << G4endl;                   342       mat->GetName() << G4endl;
345     ed << "This can happen only in Unit Tests" << 343     ed << "This can happen only in Unit Tests" << G4endl;   
346     G4Exception("G4PenelopeGammaConversionMode    344     G4Exception("G4PenelopeGammaConversionModel::SampleSecondaries()",
347           "em2019",JustWarning,ed);               345           "em2019",JustWarning,ed);
348   }                                               346   }
349       //protect file reading via autolock         347       //protect file reading via autolock
350       G4AutoLock lock(&PenelopeGammaConversion    348       G4AutoLock lock(&PenelopeGammaConversionModelMutex);
351       InitializeScreeningFunctions(mat);          349       InitializeScreeningFunctions(mat);
352       lock.unlock();                              350       lock.unlock();
353     }                                             351     }
354                                                   352 
355   // eps is the fraction of the photon energy     353   // eps is the fraction of the photon energy assigned to e- (including rest mass)
356   G4double eps = 0;                               354   G4double eps = 0;
357   G4double eki = electron_mass_c2/photonEnergy    355   G4double eki = electron_mass_c2/photonEnergy;
358                                                   356 
359   //Do it fast for photon energy < 1.1 MeV (cl    357   //Do it fast for photon energy < 1.1 MeV (close to threshold)
360   if (photonEnergy < fSmallEnergy)                358   if (photonEnergy < fSmallEnergy)
361     eps = eki + (1.0-2.0*eki)*G4UniformRand();    359     eps = eki + (1.0-2.0*eki)*G4UniformRand();
362   else                                            360   else
363     {                                             361     {
364       //Complete calculation                      362       //Complete calculation
365       G4double effC = fEffectiveCharge->find(m    363       G4double effC = fEffectiveCharge->find(mat)->second;
366       G4double alz = effC*fine_structure_const    364       G4double alz = effC*fine_structure_const;
367       G4double T = std::sqrt(2.0*eki);            365       G4double T = std::sqrt(2.0*eki);
368       G4double F00=(-1.774-1.210e1*alz+1.118e1    366       G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
369          +(8.523+7.326e1*alz-4.441e1*alz*alz)*    367          +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
370          -(1.352e1+1.211e2*alz-9.641e1*alz*alz    368          -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
371   +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T    369   +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
372                                                << 370      
373       G4double F0b = fScreeningFunction->find(    371       G4double F0b = fScreeningFunction->find(mat)->second.second;
374       G4double g0 = F0b + F00;                    372       G4double g0 = F0b + F00;
375       G4double invRad = fMaterialInvScreeningR    373       G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
376       G4double bmin = 4.0*eki/invRad;             374       G4double bmin = 4.0*eki/invRad;
377       std::pair<G4double,G4double> scree =  Ge    375       std::pair<G4double,G4double> scree =  GetScreeningFunctions(bmin);
378       G4double g1 = scree.first;                  376       G4double g1 = scree.first;
379       G4double g2 = scree.second;                 377       G4double g2 = scree.second;
380       G4double g1min = g1+g0;                     378       G4double g1min = g1+g0;
381       G4double g2min = g2+g0;                     379       G4double g2min = g2+g0;
382       G4double xr = 0.5-eki;                      380       G4double xr = 0.5-eki;
383       G4double a1 = 2.*g1min*xr*xr/3.;         << 381       G4double a1 = 2.*g1min*xr*xr/3.;      
384       G4double p1 = a1/(a1+g2min);                382       G4double p1 = a1/(a1+g2min);
385                                                   383 
386       G4bool loopAgain = false;                   384       G4bool loopAgain = false;
387       //Random sampling of eps                    385       //Random sampling of eps
388       do{                                         386       do{
389   loopAgain = false;                              387   loopAgain = false;
390   if (G4UniformRand() <= p1)                      388   if (G4UniformRand() <= p1)
391     {                                             389     {
392       G4double  ru2m1 = 2.0*G4UniformRand()-1.    390       G4double  ru2m1 = 2.0*G4UniformRand()-1.0;
393       if (ru2m1 < 0)                              391       if (ru2m1 < 0)
394         eps = 0.5-xr*std::pow(std::abs(ru2m1),    392         eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
395       else                                        393       else
396         eps = 0.5+xr*std::pow(ru2m1,1./3.);       394         eps = 0.5+xr*std::pow(ru2m1,1./3.);
397       G4double B = eki/(invRad*eps*(1.0-eps));    395       G4double B = eki/(invRad*eps*(1.0-eps));
398       scree =  GetScreeningFunctions(B);          396       scree =  GetScreeningFunctions(B);
399       g1 = scree.first;                           397       g1 = scree.first;
400       g1 = std::max(g1+g0,0.);                    398       g1 = std::max(g1+g0,0.);
401       if (G4UniformRand()*g1min > g1)          << 399       if (G4UniformRand()*g1min > g1) 
402         loopAgain = true;                         400         loopAgain = true;
403     }                                             401     }
404   else                                            402   else
405     {                                             403     {
406       eps = eki+2.0*xr*G4UniformRand();           404       eps = eki+2.0*xr*G4UniformRand();
407       G4double B = eki/(invRad*eps*(1.0-eps));    405       G4double B = eki/(invRad*eps*(1.0-eps));
408       scree =  GetScreeningFunctions(B);          406       scree =  GetScreeningFunctions(B);
409       g2 = scree.second;                          407       g2 = scree.second;
410       g2 = std::max(g2+g0,0.);                    408       g2 = std::max(g2+g0,0.);
411       if (G4UniformRand()*g2min > g2)             409       if (G4UniformRand()*g2min > g2)
412         loopAgain = true;                         410         loopAgain = true;
413     }                                          << 411     }      
414       }while(loopAgain);                          412       }while(loopAgain);
                                                   >> 413       
415     }                                             414     }
416   if (fVerboseLevel > 4)                       << 415   if (verboseLevel > 4)
417     G4cout << "Sampled eps = " << eps << G4end    416     G4cout << "Sampled eps = " << eps << G4endl;
418                                                   417 
419   G4double electronTotEnergy = eps*photonEnerg    418   G4double electronTotEnergy = eps*photonEnergy;
420   G4double positronTotEnergy = (1.0-eps)*photo    419   G4double positronTotEnergy = (1.0-eps)*photonEnergy;
421                                                << 420   
422   // Scattered electron (positron) angles. ( Z    421   // Scattered electron (positron) angles. ( Z - axis along the parent photon)
423                                                   422 
424   //electron kinematics                           423   //electron kinematics
425   G4double electronKineEnergy = std::max(0.,el << 424   G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; 
426   G4double costheta_el = G4UniformRand()*2.0-1    425   G4double costheta_el = G4UniformRand()*2.0-1.0;
427   G4double kk = std::sqrt(electronKineEnergy*(    426   G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
428   costheta_el = (costheta_el*electronTotEnergy    427   costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
429   G4double phi_el  = twopi * G4UniformRand() ;    428   G4double phi_el  = twopi * G4UniformRand() ;
430   G4double dirX_el = std::sqrt(1.-costheta_el*    429   G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
431   G4double dirY_el = std::sqrt(1.-costheta_el*    430   G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
432   G4double dirZ_el = costheta_el;                 431   G4double dirZ_el = costheta_el;
433                                                   432 
434   //positron kinematics                           433   //positron kinematics
435   G4double positronKineEnergy = std::max(0.,po    434   G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
436   G4double costheta_po = G4UniformRand()*2.0-1    435   G4double costheta_po = G4UniformRand()*2.0-1.0;
437   kk = std::sqrt(positronKineEnergy*(positronK    436   kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
438   costheta_po = (costheta_po*positronTotEnergy    437   costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
439   G4double phi_po  = twopi * G4UniformRand() ;    438   G4double phi_po  = twopi * G4UniformRand() ;
440   G4double dirX_po = std::sqrt(1.-costheta_po*    439   G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
441   G4double dirY_po = std::sqrt(1.-costheta_po*    440   G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
442   G4double dirZ_po = costheta_po;                 441   G4double dirZ_po = costheta_po;
443                                                   442 
444   // Kinematics of the created pair:              443   // Kinematics of the created pair:
445   // the electron and positron are assumed to     444   // the electron and positron are assumed to have a symetric angular
446   // distribution with respect to the Z axis a    445   // distribution with respect to the Z axis along the parent photon
447   G4double localEnergyDeposit = 0. ;              446   G4double localEnergyDeposit = 0. ;
448                                                << 447  
449   if (electronKineEnergy > 0.0)                   448   if (electronKineEnergy > 0.0)
450     {                                             449     {
451       G4ThreeVector electronDirection ( dirX_e    450       G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
452       electronDirection.rotateUz(photonDirecti    451       electronDirection.rotateUz(photonDirection);
453       G4DynamicParticle* electron = new G4Dyna    452       G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
454                  electronDirection,               453                  electronDirection,
455                  electronKineEnergy);             454                  electronKineEnergy);
456       fvect->push_back(electron);                 455       fvect->push_back(electron);
457     }                                             456     }
458   else                                            457   else
459     {                                             458     {
460       localEnergyDeposit += electronKineEnergy    459       localEnergyDeposit += electronKineEnergy;
461       electronKineEnergy = 0;                     460       electronKineEnergy = 0;
462     }                                             461     }
463                                                   462 
464   //Generate the positron. Real particle in an    463   //Generate the positron. Real particle in any case, because it will annihilate. If below
465   //threshold, produce it at rest                 464   //threshold, produce it at rest
466   if (positronKineEnergy < 0.0)                   465   if (positronKineEnergy < 0.0)
467     {                                             466     {
468       localEnergyDeposit += positronKineEnergy    467       localEnergyDeposit += positronKineEnergy;
469       positronKineEnergy = 0; //produce it at     468       positronKineEnergy = 0; //produce it at rest
470     }                                             469     }
471   G4ThreeVector positronDirection(dirX_po,dirY    470   G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
472   positronDirection.rotateUz(photonDirection);    471   positronDirection.rotateUz(photonDirection);
473   G4DynamicParticle* positron = new G4DynamicP    472   G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(),
474                   positronDirection, positronK    473                   positronDirection, positronKineEnergy);
475   fvect->push_back(positron);                     474   fvect->push_back(positron);
476                                                   475 
477   //Add rest of energy to the local energy dep    476   //Add rest of energy to the local energy deposit
478   fParticleChange->ProposeLocalEnergyDeposit(l    477   fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
479                                                << 478   
480   if (fVerboseLevel > 1)                       << 479   if (verboseLevel > 1)
481     {                                             480     {
482       G4cout << "-----------------------------    481       G4cout << "-----------------------------------------------------------" << G4endl;
483       G4cout << "Energy balance from G4Penelop    482       G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl;
484       G4cout << "Incoming photon energy: " <<     483       G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
485       G4cout << "-----------------------------    484       G4cout << "-----------------------------------------------------------" << G4endl;
486       if (electronKineEnergy)                     485       if (electronKineEnergy)
487   G4cout << "Electron (explicitly produced) "  << 486   G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV" 
488          << G4endl;                               487          << G4endl;
489       if (positronKineEnergy)                     488       if (positronKineEnergy)
490   G4cout << "Positron (not at rest) " << posit    489   G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
491       G4cout << "Rest masses of e+/- " << 2.0*    490       G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
492       if (localEnergyDeposit)                     491       if (localEnergyDeposit)
493   G4cout << "Local energy deposit " << localEn    492   G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
494       G4cout << "Total final state: " << (elec    493       G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
495             localEnergyDeposit+2.0*electron_ma    494             localEnergyDeposit+2.0*electron_mass_c2)/keV <<
496         " keV" << G4endl;                         495         " keV" << G4endl;
497       G4cout << "-----------------------------    496       G4cout << "-----------------------------------------------------------" << G4endl;
498     }                                             497     }
499  if (fVerboseLevel > 0)                        << 498  if (verboseLevel > 0)
500     {                                             499     {
501       G4double energyDiff = std::fabs(electron    500       G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
502               localEnergyDeposit+2.0*electron_    501               localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
503       if (energyDiff > 0.05*keV)                  502       if (energyDiff > 0.05*keV)
504   G4cout << "Warning from G4PenelopeGammaConve << 503   G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: " 
505          << (electronKineEnergy+positronKineEn    504          << (electronKineEnergy+positronKineEnergy+
506        localEnergyDeposit+2.0*electron_mass_c2 << 505        localEnergyDeposit+2.0*electron_mass_c2)/keV 
507          << " keV (final) vs. " << photonEnerg    506          << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
508     }                                          << 507     } 
509 }                                                 508 }
510                                                   509 
511 //....oooOO0OOooo........oooOO0OOooo........oo    510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
512                                                   511 
513 void G4PenelopeGammaConversionModel::ReadDataF    512 void G4PenelopeGammaConversionModel::ReadDataFile(const G4int Z)
514 {                                                 513 {
515   if (!IsMaster())                             << 514   if (!IsMaster())    
516       //Should not be here!                       515       //Should not be here!
517     G4Exception("G4PenelopeGammaConversionMode    516     G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
518     "em0100",FatalException,"Worker thread in  << 517     "em0100",FatalException,"Worker thread in this method");    
519                                                   518 
520   if (fVerboseLevel > 2)                       << 519  if (verboseLevel > 2)
521     {                                             520     {
522       G4cout << "G4PenelopeGammaConversionMode    521       G4cout << "G4PenelopeGammaConversionModel::ReadDataFile()" << G4endl;
523       G4cout << "Going to read Gamma Conversio    522       G4cout << "Going to read Gamma Conversion data files for Z=" << Z << G4endl;
524     }                                             523     }
525                                                << 524  
526     const char* path = G4FindDataDir("G4LEDATA << 525   char* path = getenv("G4LEDATA");
527     if(!path)                                  << 526   if (!path)
528     {                                             527     {
529       G4String excep =                         << 528       G4String excep = 
530   "G4PenelopeGammaConversionModel - G4LEDATA e    529   "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
531       G4Exception("G4PenelopeGammaConversionMo    530       G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
532       "em0006",FatalException,excep);             531       "em0006",FatalException,excep);
533       return;                                     532       return;
534     }                                             533     }
535                                                << 534  
536   /*                                              535   /*
537     Read the cross section file                   536     Read the cross section file
538   */                                              537   */
539   std::ostringstream ost;                         538   std::ostringstream ost;
540   if (Z>9)                                        539   if (Z>9)
541     ost << path << "/penelope/pairproduction/p    540     ost << path << "/penelope/pairproduction/pdgpp" << Z << ".p08";
542   else                                            541   else
543     ost << path << "/penelope/pairproduction/p    542     ost << path << "/penelope/pairproduction/pdgpp0" << Z << ".p08";
544   std::ifstream file(ost.str().c_str());          543   std::ifstream file(ost.str().c_str());
545   if (!file.is_open())                            544   if (!file.is_open())
546     {                                             545     {
547       G4String excep = "G4PenelopeGammaConvers << 546       G4String excep = "G4PenelopeGammaConversionModel - data file " + 
548   G4String(ost.str()) + " not found!";            547   G4String(ost.str()) + " not found!";
549       G4Exception("G4PenelopeGammaConversionMo    548       G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
550       "em0003",FatalException,excep);             549       "em0003",FatalException,excep);
551     }                                             550     }
552                                                   551 
553   //I have to know in advance how many points     552   //I have to know in advance how many points are in the data list
554   //to initialize the G4PhysicsFreeVector()       553   //to initialize the G4PhysicsFreeVector()
555   std::size_t ndata=0;                         << 554   size_t ndata=0;
556   G4String line;                                  555   G4String line;
557   while( getline(file, line) )                    556   while( getline(file, line) )
558     ndata++;                                      557     ndata++;
559   ndata -= 1; //remove one header line            558   ndata -= 1; //remove one header line
                                                   >> 559   //G4cout << "Found: " << ndata << " lines" << G4endl;
560                                                   560 
561   file.clear();                                   561   file.clear();
562   file.close();                                   562   file.close();
563   file.open(ost.str().c_str());                   563   file.open(ost.str().c_str());
564   G4int readZ =0;                                 564   G4int readZ =0;
565   file >> readZ;                               << 565   file >> readZ; 
566                                                   566 
567   if (fVerboseLevel > 3)                       << 567   if (verboseLevel > 3)
568     G4cout << "Element Z=" << Z << G4endl;        568     G4cout << "Element Z=" << Z << G4endl;
569                                                   569 
570   //check the right file is opened.               570   //check the right file is opened.
571   if (readZ != Z)                                 571   if (readZ != Z)
572     {                                             572     {
573       G4ExceptionDescription ed;                  573       G4ExceptionDescription ed;
574       ed << "Corrupted data file for Z=" << Z  << 574       ed << "Corrupted data file for Z=" << Z << G4endl;      
575       G4Exception("G4PenelopeGammaConversionMo    575       G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
576       "em0005",FatalException,ed);                576       "em0005",FatalException,ed);
577     }                                             577     }
578                                                   578 
579   fLogAtomicCrossSection[Z] = new G4PhysicsFre << 579   G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(ndata);
580   G4double ene=0,xs=0;                            580   G4double ene=0,xs=0;
581   for (std::size_t i=0;i<ndata;++i)            << 581   for (size_t i=0;i<ndata;i++)
582     {                                             582     {
583       file >> ene >> xs;                          583       file >> ene >> xs;
584       //dimensional quantities                    584       //dimensional quantities
585       ene *= eV;                                  585       ene *= eV;
586       xs *= barn;                                 586       xs *= barn;
587       if (xs < 1e-40*cm2) //protection against    587       if (xs < 1e-40*cm2) //protection against log(0)
588   xs = 1e-40*cm2;                                 588   xs = 1e-40*cm2;
589       fLogAtomicCrossSection[Z]->PutValue(i,G4 << 589       theVec->PutValue(i,std::log(ene),std::log(xs));      
590     }                                             590     }
591   file.close();                                   591   file.close();
592                                                   592 
                                                   >> 593   if (!logAtomicCrossSection)
                                                   >> 594     {
                                                   >> 595       G4ExceptionDescription ed;
                                                   >> 596       ed << "Problem with allocation of logAtomicCrossSection data table " << G4endl;
                                                   >> 597       G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
                                                   >> 598       "em2020",FatalException,ed);
                                                   >> 599       delete theVec;
                                                   >> 600       return;
                                                   >> 601     }
                                                   >> 602   logAtomicCrossSection->insert(std::make_pair(Z,theVec));
                                                   >> 603 
593   return;                                         604   return;
                                                   >> 605 
                                                   >> 606 }
                                                   >> 607 
                                                   >> 608 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 609 
                                                   >> 610 void G4PenelopeGammaConversionModel::InitializeScreeningRadii()
                                                   >> 611 {
                                                   >> 612   G4double temp[99] = {1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,6.4696e+01,
                                                   >> 613            6.1228e+01,5.7524e+01,5.4033e+01,5.0787e+01,4.7851e+01,4.6373e+01,
                                                   >> 614            4.5401e+01,4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,4.1586e+01,
                                                   >> 615            4.0953e+01,4.0524e+01,4.0256e+01,3.9756e+01,3.9144e+01,3.8462e+01,
                                                   >> 616            3.7778e+01,3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,3.4688e+01,
                                                   >> 617            3.4197e+01,3.3786e+01,3.3422e+01,3.3068e+01,3.2740e+01,3.2438e+01,
                                                   >> 618            3.2143e+01,3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,3.0950e+01,
                                                   >> 619            3.0758e+01,3.0561e+01,3.0285e+01,3.0097e+01,2.9832e+01,2.9581e+01,
                                                   >> 620            2.9411e+01,2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,2.8580e+01,
                                                   >> 621            2.8442e+01,2.8312e+01,2.8139e+01,2.7973e+01,2.7819e+01,2.7675e+01,
                                                   >> 622            2.7496e+01,2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,2.6516e+01,
                                                   >> 623            2.6304e+01,2.6108e+01,2.5929e+01,2.5730e+01,2.5577e+01,2.5403e+01,
                                                   >> 624            2.5245e+01,2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,2.4506e+01,
                                                   >> 625            2.4391e+01,2.4262e+01,2.4145e+01,2.4039e+01,2.3922e+01,2.3813e+01,
                                                   >> 626            2.3712e+01,2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,2.3238e+01,
                                                   >> 627            2.3139e+01,2.3048e+01,2.2967e+01,2.2833e+01,2.2694e+01,2.2624e+01,
                                                   >> 628            2.2545e+01,2.2446e+01,2.2358e+01,2.2264e+01};
                                                   >> 629 
                                                   >> 630   //copy temporary vector in class data member
                                                   >> 631   for (G4int i=0;i<99;i++)
                                                   >> 632     fAtomicScreeningRadius[i] = temp[i];
594 }                                                 633 }
595                                                   634 
596 //....oooOO0OOooo........oooOO0OOooo........oo    635 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
597                                                   636 
598 void G4PenelopeGammaConversionModel::Initializ    637 void G4PenelopeGammaConversionModel::InitializeScreeningFunctions(const G4Material* material)
599 {                                                 638 {
                                                   >> 639   /*
                                                   >> 640   if (!IsMaster())    
                                                   >> 641       //Should not be here!
                                                   >> 642     G4Exception("G4PenelopeGammaConversionModel::InitializeScreeningFunctions()",
                                                   >> 643     "em01001",FatalException,"Worker thread in this method");    
                                                   >> 644   */
                                                   >> 645 
600   // This is subroutine GPPa0 of Penelope         646   // This is subroutine GPPa0 of Penelope
601   //                                              647   //
602   // 1) calculate the effective Z for the purp    648   // 1) calculate the effective Z for the purpose
603   //                                              649   //
604   G4double zeff = 0;                              650   G4double zeff = 0;
605   G4int intZ = 0;                                 651   G4int intZ = 0;
606   G4int nElements = (G4int)material->GetNumber << 652   G4int nElements = material->GetNumberOfElements();
607   const G4ElementVector* elementVector = mater    653   const G4ElementVector* elementVector = material->GetElementVector();
608                                                   654 
609   //avoid calculations if only one building el    655   //avoid calculations if only one building element!
610   if (nElements == 1)                             656   if (nElements == 1)
611     {                                             657     {
612       zeff = (*elementVector)[0]->GetZ();         658       zeff = (*elementVector)[0]->GetZ();
613       intZ = (G4int) zeff;                        659       intZ = (G4int) zeff;
614     }                                             660     }
615   else // many elements...let's do the calcula    661   else // many elements...let's do the calculation
616     {                                             662     {
617       const G4double* fractionVector = materia    663       const G4double* fractionVector = material->GetVecNbOfAtomsPerVolume();
618                                                << 664       
619       G4double atot = 0;                          665       G4double atot = 0;
620       for (G4int i=0;i<nElements;i++)             666       for (G4int i=0;i<nElements;i++)
621   {                                               667   {
622     G4double Zelement = (*elementVector)[i]->G    668     G4double Zelement = (*elementVector)[i]->GetZ();
623     G4double Aelement = (*elementVector)[i]->G    669     G4double Aelement = (*elementVector)[i]->GetAtomicMassAmu();
624     atot += Aelement*fractionVector[i];           670     atot += Aelement*fractionVector[i];
625     zeff += Zelement*Aelement*fractionVector[i    671     zeff += Zelement*Aelement*fractionVector[i]; //average with the number of nuclei
626   }                                               672   }
627       atot /= material->GetTotNbOfAtomsPerVolu    673       atot /= material->GetTotNbOfAtomsPerVolume();
628       zeff /= (material->GetTotNbOfAtomsPerVol    674       zeff /= (material->GetTotNbOfAtomsPerVolume()*atot);
629                                                << 675       
630       intZ = (G4int) (zeff+0.25);                 676       intZ = (G4int) (zeff+0.25);
631       if (intZ <= 0)                              677       if (intZ <= 0)
632   intZ = 1;                                       678   intZ = 1;
633       if (intZ > fMaxZ)                        << 679       if (intZ > 99)
634   intZ = fMaxZ;                                << 680   intZ = 99;
635     }                                             681     }
636                                                   682 
637   if (fEffectiveCharge)                           683   if (fEffectiveCharge)
638     fEffectiveCharge->insert(std::make_pair(ma    684     fEffectiveCharge->insert(std::make_pair(material,zeff));
639                                                   685 
640   //                                              686   //
641   // 2) Calculate Coulomb Correction              687   // 2) Calculate Coulomb Correction
642   //                                              688   //
643   G4double alz = fine_structure_const*zeff;       689   G4double alz = fine_structure_const*zeff;
644   G4double alzSquared = alz*alz;                  690   G4double alzSquared = alz*alz;
645   G4double fc =  alzSquared*(0.202059-alzSquar    691   G4double fc =  alzSquared*(0.202059-alzSquared*
646            (0.03693-alzSquared*                   692            (0.03693-alzSquared*
647             (0.00835-alzSquared*(0.00201-alzSq    693             (0.00835-alzSquared*(0.00201-alzSquared*
648                (0.00049-alzSquared*               694                (0.00049-alzSquared*
649                 (0.00012-alzSquared*0.00003)))    695                 (0.00012-alzSquared*0.00003)))))
650            +1.0/(alzSquared+1.0));                696            +1.0/(alzSquared+1.0));
651   //                                              697   //
652   // 3) Screening functions and low-energy cor    698   // 3) Screening functions and low-energy corrections
653   //                                              699   //
654   G4double matRadius = 2.0/ fAtomicScreeningRa << 700   G4double matRadius = 2.0/ fAtomicScreeningRadius[intZ-1];
655   if (fMaterialInvScreeningRadius)                701   if (fMaterialInvScreeningRadius)
656     fMaterialInvScreeningRadius->insert(std::m    702     fMaterialInvScreeningRadius->insert(std::make_pair(material,matRadius));
657                                                   703 
658   std::pair<G4double,G4double> myPair(0,0);       704   std::pair<G4double,G4double> myPair(0,0);
659   G4double f0a = 4.0*G4Log(fAtomicScreeningRad << 705   G4double f0a = 4.0*std::log(fAtomicScreeningRadius[intZ-1]);
660   G4double f0b = f0a - 4.0*fc;                    706   G4double f0b = f0a - 4.0*fc;
661   myPair.first = f0a;                             707   myPair.first = f0a;
662   myPair.second = f0b;                            708   myPair.second = f0b;
663                                                   709 
664   if (fScreeningFunction)                         710   if (fScreeningFunction)
665     fScreeningFunction->insert(std::make_pair(    711     fScreeningFunction->insert(std::make_pair(material,myPair));
666                                                   712 
667   if (fVerboseLevel > 2)                       << 713   if (verboseLevel > 2)
668     {                                             714     {
669       G4cout << "Average Z for material " << m << 715       G4cout << "Average Z for material " << material->GetName() << " = " << 
670   zeff << G4endl;                                 716   zeff << G4endl;
671       G4cout << "Effective radius for material << 717       G4cout << "Effective radius for material " << material->GetName() << " = " << 
672   fAtomicScreeningRadius[intZ] << " m_e*c/hbar << 718   fAtomicScreeningRadius[intZ-1] << " m_e*c/hbar --> BCB = " << 
673   matRadius << G4endl;                            719   matRadius << G4endl;
674       G4cout << "Screening parameters F0 for m << 720       G4cout << "Screening parameters F0 for material " << material->GetName() << " = " << 
675   f0a << "," << f0b << G4endl;                    721   f0a << "," << f0b << G4endl;
676     }                                             722     }
677   return;                                         723   return;
678 }                                                 724 }
679                                                   725 
680 //....oooOO0OOooo........oooOO0OOooo........oo    726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
681                                                   727 
682 std::pair<G4double,G4double>                   << 728 std::pair<G4double,G4double> 
683 G4PenelopeGammaConversionModel::GetScreeningFu    729 G4PenelopeGammaConversionModel::GetScreeningFunctions(G4double B)
684 {                                                 730 {
685   // This is subroutine SCHIFF of Penelope        731   // This is subroutine SCHIFF of Penelope
686   //                                              732   //
687   // Screening Functions F1(B) and F2(B) in th << 733   // Screening Functions F1(B) and F2(B) in the Bethe-Heitler differential cross 
688   // section for pair production                  734   // section for pair production
689   //                                              735   //
690   std::pair<G4double,G4double> result(0.,0.);     736   std::pair<G4double,G4double> result(0.,0.);
691   G4double BSquared = B*B;                        737   G4double BSquared = B*B;
692   G4double f1 = 2.0-2.0*G4Log(1.0+BSquared);   << 738   G4double f1 = 2.0-2.0*std::log(1.0+BSquared);
693   G4double f2 = f1 - 6.66666666e-1; // (-2/3)     739   G4double f2 = f1 - 6.66666666e-1; // (-2/3)
694   if (B < 1.0e-10)                                740   if (B < 1.0e-10)
695     f1 = f1-twopi*B;                              741     f1 = f1-twopi*B;
696   else                                            742   else
697     {                                             743     {
698       G4double a0 = 4.0*B*std::atan(1./B);        744       G4double a0 = 4.0*B*std::atan(1./B);
699       f1 = f1 - a0;                               745       f1 = f1 - a0;
700       f2 += 2.0*BSquared*(4.0-a0-3.0*G4Log((1. << 746       f2 += 2.0*BSquared*(4.0-a0-3.0*std::log((1.0+BSquared)/BSquared));
701     }                                             747     }
702   G4double g1 = 0.5*(3.0*f1-f2);                  748   G4double g1 = 0.5*(3.0*f1-f2);
703   G4double g2 = 0.25*(3.0*f1+f2);                 749   G4double g2 = 0.25*(3.0*f1+f2);
704                                                   750 
705   result.first = g1;                              751   result.first = g1;
706   result.second = g2;                             752   result.second = g2;
707                                                   753 
708   return result;                                  754   return result;
709 }                                                 755 }
710                                                   756 
711 //....oooOO0OOooo........oooOO0OOooo........oo    757 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
712                                                   758 
713 void G4PenelopeGammaConversionModel::SetPartic    759 void G4PenelopeGammaConversionModel::SetParticle(const G4ParticleDefinition* p)
714 {                                                 760 {
715   if(!fParticle) {                                761   if(!fParticle) {
716     fParticle = p;                             << 762     fParticle = p;  
717   }                                               763   }
718 }                                                 764 }
                                                   >> 765 
719                                                   766