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.4)


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