Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungFS.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/G4PenelopeBremsstrahlungFS.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungFS.cc (Version 10.4.p3)


  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$
 26 //                                                 27 //
 27 // Author: Luciano Pandola                         28 // Author: Luciano Pandola
 28 //                                                 29 //
 29 // History:                                        30 // History:
 30 // --------                                        31 // --------
 31 // 23 Nov 2010   L Pandola    First complete i     32 // 23 Nov 2010   L Pandola    First complete implementation
 32 // 02 May 2011   L.Pandola    Remove dependenc     33 // 02 May 2011   L.Pandola    Remove dependency on CLHEP::HepMatrix
 33 // 24 May 2011   L.Pandola    Renamed (make v2     34 // 24 May 2011   L.Pandola    Renamed (make v2008 as default Penelope)
 34 // 03 Oct 2013   L.Pandola    Migration to MT      35 // 03 Oct 2013   L.Pandola    Migration to MT
 35 // 30 Oct 2013   L.Pandola    Use G4Cache to a     36 // 30 Oct 2013   L.Pandola    Use G4Cache to avoid new/delete of the
 36 //                             data vector on      37 //                             data vector on the fly in SampleGammaEnergy()
 37 //                                                 38 //
 38                                                    39 
 39 //....oooOO0OOooo........oooOO0OOooo........oo     40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 40 #include "G4PenelopeBremsstrahlungFS.hh"           41 #include "G4PenelopeBremsstrahlungFS.hh"
 41 #include "G4PhysicsFreeVector.hh"                  42 #include "G4PhysicsFreeVector.hh"
 42 #include "G4PhysicsTable.hh"                       43 #include "G4PhysicsTable.hh"
 43 #include "G4Material.hh"                           44 #include "G4Material.hh"
 44 #include "Randomize.hh"                            45 #include "Randomize.hh"
 45 #include "G4AutoDelete.hh"                         46 #include "G4AutoDelete.hh"
 46 #include "G4Exp.hh"                                47 #include "G4Exp.hh"
 47 #include "G4PhysicalConstants.hh"                  48 #include "G4PhysicalConstants.hh"
 48 #include "G4SystemOfUnits.hh"                      49 #include "G4SystemOfUnits.hh"
 49                                                    50 
 50 //....oooOO0OOooo........oooOO0OOooo........oo     51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 51                                                    52 
 52 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstr     53 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstrahlungFS(G4int verbosity) :
 53   fReducedXSTable(nullptr),fEffectiveZSq(nullp <<  54   theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
 54   fPBcut(nullptr),fVerbosity(verbosity)        <<  55   thePBcut(0),fVerbosity(verbosity)
 55 {                                                  56 {
 56   fCache.Put(0);                                   57   fCache.Put(0);
 57   G4double tempvector[fNBinsX] =               <<  58   G4double tempvector[nBinsX] =
 58     {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15     59     {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0,
 59     0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6     60     0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0,
 60     0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0     61     0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0,
 61     0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e     62     0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
 62                                                    63 
 63   for (std::size_t ix=0;ix<fNBinsX;++ix)       <<  64   for (size_t ix=0;ix<nBinsX;ix++)
 64     theXGrid[ix] = tempvector[ix];                 65     theXGrid[ix] = tempvector[ix];
 65                                                    66 
 66   for (std::size_t i=0;i<fNBinsE;++i)          <<  67   for (size_t i=0;i<nBinsE;i++)
 67     theEGrid[i] = 0.;                              68     theEGrid[i] = 0.;
 68                                                    69 
 69   fElementData = new std::map<G4int,G4DataVect <<  70   theElementData = new std::map<G4int,G4DataVector*>;
 70 }                                                  71 }
 71                                                    72 
 72 //....oooOO0OOooo........oooOO0OOooo........oo     73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 73                                                    74 
 74 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsst     75 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsstrahlungFS()
 75 {                                                  76 {
 76   ClearTables();                                   77   ClearTables();
 77                                                    78 
 78   //The G4Physics*Vector pointers contained in     79   //The G4Physics*Vector pointers contained in the fCache are automatically deleted by
 79   //the G4AutoDelete so there is no need to ta     80   //the G4AutoDelete so there is no need to take care of them manually
 80                                                    81 
 81   //Clear manually fElementData                <<  82   //Clear manually theElementData
 82   if (fElementData)                            <<  83   if (theElementData)
 83     {                                              84     {
 84       for (auto& item : (*fElementData))       <<  85       for (auto& item : (*theElementData)) 
 85   delete item.second;                              86   delete item.second;
 86       delete fElementData;                     <<  87       delete theElementData;
 87       fElementData = nullptr;                  <<  88       theElementData = nullptr;
 88     }                                              89     }
                                                   >>  90 
 89 }                                                  91 }
 90                                                    92 
 91 //....oooOO0OOooo........oooOO0OOooo........oo     93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
 92                                                    94 
 93                                                    95 
 94 void G4PenelopeBremsstrahlungFS::ClearTables(G     96 void G4PenelopeBremsstrahlungFS::ClearTables(G4bool isMaster)
 95 {                                                  97 {
 96   //Just to check                                  98   //Just to check
 97   if (!isMaster)                                   99   if (!isMaster)
 98     G4Exception("G4PenelopeBremsstrahlungFS::C    100     G4Exception("G4PenelopeBremsstrahlungFS::ClearTables()",
 99     "em0100",FatalException,"Worker thread in     101     "em0100",FatalException,"Worker thread in this method");
100                                                   102 
101   if (fReducedXSTable)                         << 103 
                                                   >> 104   if (theReducedXSTable)
102     {                                             105     {
103       for (auto& item : (*fReducedXSTable))    << 106       for (auto& item : (*theReducedXSTable))
104   {                                               107   {
105     G4PhysicsTable* tab = item.second;         << 108     //G4PhysicsTable* tab = item.second;
106     tab->clearAndDestroy();                    << 109     //tab->clearAndDestroy();
107     delete tab;                                << 110     delete item.second;
108   }                                               111   }
109       fReducedXSTable->clear();                << 112       delete theReducedXSTable;
110       delete fReducedXSTable;                  << 113       theReducedXSTable = nullptr;
111       fReducedXSTable = nullptr;               << 
112     }                                             114     }
113                                                   115 
114   if (fSamplingTable)                          << 116   if (theSamplingTable)
115     {                                             117     {
116       for (auto& item : (*fSamplingTable))     << 118       for (auto& item : (*theSamplingTable))
117   {                                               119   {
118     G4PhysicsTable* tab = item.second;         << 120     //G4PhysicsTable* tab = item.second;
119     tab->clearAndDestroy();                    << 121     // tab->clearAndDestroy();
120           delete tab;                          << 122           delete item.second;
121   }                                               123   }
122       fSamplingTable->clear();                 << 124       delete theSamplingTable;
123       delete fSamplingTable;                   << 125       theSamplingTable = nullptr;
124       fSamplingTable = nullptr;                << 
125     }                                             126     }
126   if (fPBcut)                                  << 127   if (thePBcut)
127     {                                             128     {
128       /*                                          129       /*
129   std::map< std::pair<const G4Material*,G4doub    130   std::map< std::pair<const G4Material*,G4double> ,G4PhysicsFreeVector*>::iterator kk;
130   for (kk=fPBcut->begin(); kk != fPBcut->end() << 131   for (kk=thePBcut->begin(); kk != thePBcut->end(); kk++)
131   delete kk->second;                              132   delete kk->second;
132       */                                          133       */
133       delete fPBcut;                           << 134       delete thePBcut;
134       fPBcut = nullptr;                        << 135       thePBcut = nullptr;
135     }                                             136     }
136                                                   137 
137   if (fEffectiveZSq)                           << 138 
                                                   >> 139   if (theEffectiveZSq)
138     {                                             140     {
139       delete fEffectiveZSq;                    << 141       delete theEffectiveZSq;
140       fEffectiveZSq = nullptr;                 << 142       theEffectiveZSq = nullptr;
141     }                                             143     }
142                                                   144 
143  return;                                          145  return;
144 }                                                 146 }
145                                                   147 
146 //....oooOO0OOooo........oooOO0OOooo........oo    148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147                                                   149 
148 G4double G4PenelopeBremsstrahlungFS::GetEffect    150 G4double G4PenelopeBremsstrahlungFS::GetEffectiveZSquared(const G4Material* material) const
149 {                                                 151 {
150   if (!fEffectiveZSq)                          << 152   if (!theEffectiveZSq)
151     {                                             153     {
152       G4ExceptionDescription ed;                  154       G4ExceptionDescription ed;
153       ed << "The container for the <Z^2> value    155       ed << "The container for the <Z^2> values is not initialized" << G4endl;
154       G4Exception("G4PenelopeBremsstrahlungFS:    156       G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
155       "em2007",FatalException,ed);                157       "em2007",FatalException,ed);
156       return 0;                                   158       return 0;
157     }                                             159     }
158   //found in the table: return it                 160   //found in the table: return it
159   if (fEffectiveZSq->count(material))          << 161   if (theEffectiveZSq->count(material))
160     return fEffectiveZSq->find(material)->seco << 162     return theEffectiveZSq->find(material)->second;
161   else                                            163   else
162     {                                             164     {
163       G4ExceptionDescription ed;                  165       G4ExceptionDescription ed;
164       ed << "The value of  <Z^2> is not proper    166       ed << "The value of  <Z^2> is not properly set for material " <<
165   material->GetName() << G4endl;                  167   material->GetName() << G4endl;
166       //requires running of BuildScaledXSTable    168       //requires running of BuildScaledXSTable()
167       G4Exception("G4PenelopeBremsstrahlungFS:    169       G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
168       "em2008",FatalException,ed);                170       "em2008",FatalException,ed);
169     }                                             171     }
170   return 0;                                       172   return 0;
171 }                                                 173 }
172                                                   174 
173 //....oooOO0OOooo........oooOO0OOooo........oo    175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
174                                                   176 
175 void G4PenelopeBremsstrahlungFS::BuildScaledXS    177 void G4PenelopeBremsstrahlungFS::BuildScaledXSTable(const G4Material* material,
176                 G4double cut,G4bool isMaster)     178                 G4double cut,G4bool isMaster)
177 {                                                 179 {
178   //Corresponds to subroutines EBRaW and EBRaR    180   //Corresponds to subroutines EBRaW and EBRaR of PENELOPE
179   /*                                              181   /*
180     This method generates the table of the sca    182     This method generates the table of the scaled energy-loss cross section from
181     bremsstrahlung emission for the given mate    183     bremsstrahlung emission for the given material. Original data are read from
182     file. The table is normalized according to    184     file. The table is normalized according to the Berger-Seltzer cross section.
183   */                                              185   */
184                                                   186 
185   //Just to check                                 187   //Just to check
186   if (!isMaster)                                  188   if (!isMaster)
187     G4Exception("G4PenelopeBremsstrahlungFS::B    189     G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
188     "em0100",FatalException,"Worker thread in     190     "em0100",FatalException,"Worker thread in this method");
189                                                   191 
190   if (fVerbosity > 2)                             192   if (fVerbosity > 2)
191     {                                             193     {
192       G4cout << "Entering in G4PenelopeBremsst    194       G4cout << "Entering in G4PenelopeBremsstrahlungFS::BuildScaledXSTable for " <<
193   material->GetName() << G4endl;                  195   material->GetName() << G4endl;
194       G4cout << "Threshold = " << cut/keV << "    196       G4cout << "Threshold = " << cut/keV << " keV, isMaster= " << isMaster <<
195   G4endl;                                         197   G4endl;
196     }                                             198     }
197                                                   199 
198   //This method should be accessed by the mast    200   //This method should be accessed by the master only
199   if (!fSamplingTable)                         << 201   if (!theSamplingTable)
200     fSamplingTable =                           << 202     theSamplingTable =
201       new std::map< std::pair<const G4Material    203       new std::map< std::pair<const G4Material*,G4double> , G4PhysicsTable*>;
202   if (!fPBcut)                                 << 204   if (!thePBcut)
203     fPBcut =                                   << 205     thePBcut =
204       new std::map< std::pair<const G4Material    206       new std::map< std::pair<const G4Material*,G4double> , G4PhysicsFreeVector* >;
205                                                   207 
206   //check if the container exists (if not, cre    208   //check if the container exists (if not, create it)
207   if (!fReducedXSTable)                        << 209   if (!theReducedXSTable)
208     fReducedXSTable = new std::map< std::pair< << 210     theReducedXSTable = new std::map< std::pair<const G4Material*,G4double> ,
209     G4PhysicsTable*>;                             211     G4PhysicsTable*>;
210   if (!fEffectiveZSq)                          << 212   if (!theEffectiveZSq)
211     fEffectiveZSq = new std::map<const G4Mater << 213     theEffectiveZSq = new std::map<const G4Material*,G4double>;
212                                                   214 
213   //****************************************** << 215 
                                                   >> 216 
                                                   >> 217  //*********************************************************************
214   //Determine the equivalent atomic number <Z^    218   //Determine the equivalent atomic number <Z^2>
215   //******************************************    219   //*********************************************************************
216   std::vector<G4double> *StechiometricFactors     220   std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
217   std::size_t nElements = material->GetNumberO << 221   G4int nElements = material->GetNumberOfElements();
218   const G4ElementVector* elementVector = mater    222   const G4ElementVector* elementVector = material->GetElementVector();
219   const G4double* fractionVector = material->G    223   const G4double* fractionVector = material->GetFractionVector();
220   for (std::size_t i=0;i<nElements;i++)        << 224   for (G4int i=0;i<nElements;i++)
221     {                                             225     {
222       G4double fraction = fractionVector[i];      226       G4double fraction = fractionVector[i];
223       G4double atomicWeigth = (*elementVector)    227       G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
224       StechiometricFactors->push_back(fraction    228       StechiometricFactors->push_back(fraction/atomicWeigth);
225     }                                             229     }
226   //Find max                                      230   //Find max
227   G4double MaxStechiometricFactor = 0.;           231   G4double MaxStechiometricFactor = 0.;
228   for (std::size_t i=0;i<nElements;i++)        << 232   for (G4int i=0;i<nElements;i++)
229     {                                             233     {
230       if ((*StechiometricFactors)[i] > MaxStec    234       if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
231         MaxStechiometricFactor = (*Stechiometr    235         MaxStechiometricFactor = (*StechiometricFactors)[i];
232     }                                             236     }
233   //Normalize                                     237   //Normalize
234   for (std::size_t i=0;i<nElements;i++)        << 238   for (G4int i=0;i<nElements;i++)
235     (*StechiometricFactors)[i] /=  MaxStechiom    239     (*StechiometricFactors)[i] /=  MaxStechiometricFactor;
236                                                   240 
237   G4double sumz2 = 0;                             241   G4double sumz2 = 0;
238   G4double sums = 0;                              242   G4double sums = 0;
239   for (std::size_t i=0;i<nElements;i++)        << 243   for (G4int i=0;i<nElements;i++)
240     {                                             244     {
241       G4double Z = (*elementVector)[i]->GetZ()    245       G4double Z = (*elementVector)[i]->GetZ();
242       sumz2 += (*StechiometricFactors)[i]*Z*Z;    246       sumz2 += (*StechiometricFactors)[i]*Z*Z;
243       sums  += (*StechiometricFactors)[i];        247       sums  += (*StechiometricFactors)[i];
244     }                                             248     }
245   G4double ZBR2 = sumz2/sums;                     249   G4double ZBR2 = sumz2/sums;
246                                                   250 
247   fEffectiveZSq->insert(std::make_pair(materia << 251   theEffectiveZSq->insert(std::make_pair(material,ZBR2));
                                                   >> 252 
248                                                   253 
249   //******************************************    254   //*********************************************************************
250   // loop on elements and read data files         255   // loop on elements and read data files
251   //******************************************    256   //*********************************************************************
252   G4DataVector* tempData = new G4DataVector(fN << 257   G4DataVector* tempData = new G4DataVector(nBinsE);
253   G4DataVector* tempMatrix = new G4DataVector( << 258   G4DataVector* tempMatrix = new G4DataVector(nBinsE*nBinsX,0.);
254                                                   259 
255   for (std::size_t iel=0;iel<nElements;iel++)  << 260   for (G4int iel=0;iel<nElements;iel++)
256     {                                             261     {
257       G4double Z = (*elementVector)[iel]->GetZ    262       G4double Z = (*elementVector)[iel]->GetZ();
258       G4int iZ = (G4int) Z;                       263       G4int iZ = (G4int) Z;
259       G4double wgt = (*StechiometricFactors)[i    264       G4double wgt = (*StechiometricFactors)[iel]*Z*Z/ZBR2;
260                                                << 265       //
                                                   >> 266 
261       //the element is not already loaded         267       //the element is not already loaded
262       if (!fElementData->count(iZ))            << 268       if (!theElementData->count(iZ))
263   {                                               269   {
264     ReadDataFile(iZ);                             270     ReadDataFile(iZ);
265     if (!fElementData->count(iZ))              << 271     if (!theElementData->count(iZ))
266       {                                           272       {
267         G4ExceptionDescription ed;                273         G4ExceptionDescription ed;
268         ed << "Error in G4PenelopeBremsstrahlu    274         ed << "Error in G4PenelopeBremsstrahlungFS::BuildScaledXSTable" << G4endl;
269         ed << "Unable to retrieve data for ele    275         ed << "Unable to retrieve data for element " << iZ << G4endl;
270         G4Exception("G4PenelopeBremsstrahlungF    276         G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
271         "em2009",FatalException,ed);              277         "em2009",FatalException,ed);
272       }                                           278       }
273   }                                               279   }
274                                                   280 
275       G4DataVector* atomData = fElementData->f << 281       G4DataVector* atomData = theElementData->find(iZ)->second;
276                                                   282 
277       for (std::size_t ie=0;ie<fNBinsE;++ie)   << 283       for (size_t ie=0;ie<nBinsE;ie++)
278   {                                               284   {
279     (*tempData)[ie] += wgt*(*atomData)[ie*(fNB << 285     (*tempData)[ie] += wgt*(*atomData)[ie*(nBinsX+1)+nBinsX]; //last column contains total XS
280     for (std::size_t ix=0;ix<fNBinsX;++ix)     << 286     for (size_t ix=0;ix<nBinsX;ix++)
281       (*tempMatrix)[ie*fNBinsX+ix] += wgt*(*at << 287       (*tempMatrix)[ie*nBinsX+ix] += wgt*(*atomData)[ie*(nBinsX+1)+ix];
282   }                                               288   }
283     }                                             289     }
284                                                   290 
285   //******************************************    291   //*********************************************************************
286   // the total energy loss spectrum is re-norm    292   // the total energy loss spectrum is re-normalized to reproduce the total
287   // scaled cross section of Berger and Seltze    293   // scaled cross section of Berger and Seltzer
288   //******************************************    294   //*********************************************************************
289   for (std::size_t ie=0;ie<fNBinsE;++ie)       << 295   for (size_t ie=0;ie<nBinsE;ie++)
290     {                                             296     {
291       //for each energy, calculate integral of    297       //for each energy, calculate integral of dSigma/dx over dx
292       G4double* tempData2 = new G4double[fNBin << 298       G4double* tempData2 = new G4double[nBinsX];
293       for (std::size_t ix=0;ix<fNBinsX;++ix)   << 299       for (size_t ix=0;ix<nBinsX;ix++)
294   tempData2[ix] = (*tempMatrix)[ie*fNBinsX+ix] << 300   tempData2[ix] = (*tempMatrix)[ie*nBinsX+ix];
295       G4double rsum = GetMomentumIntegral(temp    301       G4double rsum = GetMomentumIntegral(tempData2,1.0,0);
296       delete[] tempData2;                         302       delete[] tempData2;
297       G4double fact = millibarn*(theEGrid[ie]+    303       G4double fact = millibarn*(theEGrid[ie]+electron_mass_c2)*(1./fine_structure_const)/
298   (classic_electr_radius*classic_electr_radius    304   (classic_electr_radius*classic_electr_radius*(theEGrid[ie]+2.0*electron_mass_c2));
299       G4double fnorm = (*tempData)[ie]/(rsum*f    305       G4double fnorm = (*tempData)[ie]/(rsum*fact);
300       G4double TST = 100.*std::fabs(fnorm-1.0)    306       G4double TST = 100.*std::fabs(fnorm-1.0);
301       if (TST > 1.0)                              307       if (TST > 1.0)
302   {                                               308   {
303     G4ExceptionDescription ed;                    309     G4ExceptionDescription ed;
304     ed << "G4PenelopeBremsstrahlungFS. Corrupt    310     ed << "G4PenelopeBremsstrahlungFS. Corrupted data files?" << G4endl;
305     G4cout << "TST= " << TST << "; fnorm = " <    311     G4cout << "TST= " << TST << "; fnorm = " << fnorm << G4endl;
306     G4cout << "rsum = " << rsum << G4endl;        312     G4cout << "rsum = " << rsum << G4endl;
307     G4cout << "fact = " << fact << G4endl;        313     G4cout << "fact = " << fact << G4endl;
308     G4cout << ie << " " << theEGrid[ie]/keV <<    314     G4cout << ie << " " << theEGrid[ie]/keV << " " << (*tempData)[ie]/barn << G4endl;
309     G4Exception("G4PenelopeBremsstrahlungFS::B    315     G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
310           "em2010",FatalException,ed);            316           "em2010",FatalException,ed);
311   }                                               317   }
312       for (std::size_t ix=0;ix<fNBinsX;++ix)   << 318       for (size_t ix=0;ix<nBinsX;ix++)
313   (*tempMatrix)[ie*fNBinsX+ix] *= fnorm;       << 319   (*tempMatrix)[ie*nBinsX+ix] *= fnorm;
314     }                                             320     }
315                                                   321 
316   //******************************************    322   //*********************************************************************
317   // create and fill the tables                   323   // create and fill the tables
318   //******************************************    324   //*********************************************************************
319   G4PhysicsTable* thePhysicsTable = new G4Phys    325   G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
320   // the table will contain 32 G4PhysicsFreeVe    326   // the table will contain 32 G4PhysicsFreeVectors with different
321   // values of x. Each of the G4PhysicsFreeVec    327   // values of x. Each of the G4PhysicsFreeVectors has a profile of
322   // log(XS) vs. log(E)                           328   // log(XS) vs. log(E)
323                                                   329 
324   //reserve space of the vectors. Everything i    330   //reserve space of the vectors. Everything is log-log
325   //I add one extra "fake" point at low energy    331   //I add one extra "fake" point at low energy, since the Penelope
326   //table starts at 1 keV                         332   //table starts at 1 keV
327   for (std::size_t i=0;i<fNBinsX;++i)          << 333   for (size_t i=0;i<nBinsX;i++)
328     thePhysicsTable->push_back(new G4PhysicsFr << 334     thePhysicsTable->push_back(new G4PhysicsFreeVector(nBinsE+1));
329                                                   335 
330   for (std::size_t ix=0;ix<fNBinsX;++ix)       << 336   for (size_t ix=0;ix<nBinsX;ix++)
331     {                                             337     {
332       G4PhysicsFreeVector* theVec =               338       G4PhysicsFreeVector* theVec =
333   (G4PhysicsFreeVector*) ((*thePhysicsTable)[i    339   (G4PhysicsFreeVector*) ((*thePhysicsTable)[ix]);
334       for (std::size_t ie=0;ie<fNBinsE;++ie)   << 340       for (size_t ie=0;ie<nBinsE;ie++)
335   {                                               341   {
336     G4double logene = G4Log(theEGrid[ie]);     << 342     G4double logene = std::log(theEGrid[ie]);
337     G4double aValue = (*tempMatrix)[ie*fNBinsX << 343     G4double aValue = (*tempMatrix)[ie*nBinsX+ix];
338     if (aValue < 1e-20*millibarn) //protection    344     if (aValue < 1e-20*millibarn) //protection against log(0)
339       aValue = 1e-20*millibarn;                   345       aValue = 1e-20*millibarn;
340     theVec->PutValues(ie+1,logene,G4Log(aValue << 346     theVec->PutValue(ie+1,logene,std::log(aValue));
341   }                                               347   }
342       //Add fake point at 1 eV using an extrap    348       //Add fake point at 1 eV using an extrapolation with the derivative
343       //at the first valid point (Penelope app    349       //at the first valid point (Penelope approach)
344       G4double derivative = ((*theVec)[2]-(*th    350       G4double derivative = ((*theVec)[2]-(*theVec)[1])/(theVec->Energy(2) - theVec->Energy(1));
345       G4double log1eV = G4Log(1*eV);           << 351       G4double log1eV = std::log(1*eV);
346       G4double val1eV = (*theVec)[1]+derivativ    352       G4double val1eV = (*theVec)[1]+derivative*(log1eV-theVec->Energy(1));
347       //fake point at very low energy             353       //fake point at very low energy
348       theVec->PutValues(0,log1eV,val1eV);      << 354       theVec->PutValue(0,log1eV,val1eV);
349     }                                             355     }
350   std::pair<const G4Material*,G4double> theKey    356   std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
351   fReducedXSTable->insert(std::make_pair(theKe << 357   theReducedXSTable->insert(std::make_pair(theKey,thePhysicsTable));
352                                                   358 
353   delete StechiometricFactors;                    359   delete StechiometricFactors;
354   delete tempData;                                360   delete tempData;
355   delete tempMatrix;                              361   delete tempMatrix;
356                                                   362 
357   //Do here also the initialization of the ene    363   //Do here also the initialization of the energy sampling
358   if (!(fSamplingTable->count(theKey)))        << 364   if (!(theSamplingTable->count(theKey)))
359     InitializeEnergySampling(material,cut);       365     InitializeEnergySampling(material,cut);
360                                                   366 
361   return;                                         367   return;
362 }                                                 368 }
363                                                   369 
364 //....oooOO0OOooo........oooOO0OOooo........oo    370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
365                                                   371 
366 void G4PenelopeBremsstrahlungFS::ReadDataFile(    372 void G4PenelopeBremsstrahlungFS::ReadDataFile(G4int Z)
367 {                                                 373 {
368   const char* path = G4FindDataDir("G4LEDATA") << 374 
                                                   >> 375   char* path = getenv("G4LEDATA");
369   if (!path)                                      376   if (!path)
370     {                                             377     {
371       G4String excep = "G4PenelopeBremsstrahlu    378       G4String excep = "G4PenelopeBremsstrahlungFS - G4LEDATA environment variable not set!";
372       G4Exception("G4PenelopeBremsstrahlungFS:    379       G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
373       "em0006",FatalException,excep);             380       "em0006",FatalException,excep);
374       return;                                     381       return;
375     }                                             382     }
376   /*                                              383   /*
377     Read the cross section file                   384     Read the cross section file
378   */                                              385   */
379   std::ostringstream ost;                         386   std::ostringstream ost;
380   if (Z>9)                                        387   if (Z>9)
381     ost << path << "/penelope/bremsstrahlung/p    388     ost << path << "/penelope/bremsstrahlung/pdebr" << Z << ".p08";
382   else                                            389   else
383     ost << path << "/penelope/bremsstrahlung/p    390     ost << path << "/penelope/bremsstrahlung/pdebr0" << Z << ".p08";
384   std::ifstream file(ost.str().c_str());          391   std::ifstream file(ost.str().c_str());
385   if (!file.is_open())                            392   if (!file.is_open())
386     {                                             393     {
387       G4String excep = "G4PenelopeBremsstrahlu    394       G4String excep = "G4PenelopeBremsstrahlungFS - data file " +
388   G4String(ost.str()) + " not found!";            395   G4String(ost.str()) + " not found!";
389       G4Exception("G4PenelopeBremsstrahlungFS:    396       G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
390       "em0003",FatalException,excep);             397       "em0003",FatalException,excep);
391       return;                                     398       return;
392     }                                             399     }
393                                                   400 
394   G4int readZ =0;                                 401   G4int readZ =0;
395   file >> readZ;                                  402   file >> readZ;
396                                                   403 
397   //check the right file is opened.               404   //check the right file is opened.
398   if (readZ != Z)                                 405   if (readZ != Z)
399     {                                             406     {
400       G4ExceptionDescription ed;                  407       G4ExceptionDescription ed;
401       ed << "Corrupted data file for Z=" << Z     408       ed << "Corrupted data file for Z=" << Z << G4endl;
402       G4Exception("G4PenelopeBremsstrahlungFS:    409       G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
403       "em0005",FatalException,ed);                410       "em0005",FatalException,ed);
404       return;                                     411       return;
405     }                                             412     }
406                                                   413 
407   G4DataVector* theMatrix = new G4DataVector(f << 414   G4DataVector* theMatrix = new G4DataVector(nBinsE*(nBinsX+1),0.); //initialized with zeros
408                                                   415 
409   for (std::size_t ie=0;ie<fNBinsE;++ie)       << 416   for (size_t ie=0;ie<nBinsE;ie++)
410     {                                             417     {
411       G4double myDouble = 0;                      418       G4double myDouble = 0;
412       file >> myDouble; //energy (eV)             419       file >> myDouble; //energy (eV)
413       if (!theEGrid[ie]) //fill only the first    420       if (!theEGrid[ie]) //fill only the first time
414   theEGrid[ie] = myDouble*eV;                     421   theEGrid[ie] = myDouble*eV;
415       //                                          422       //
416       for (std::size_t ix=0;ix<fNBinsX;++ix)   << 423       for (size_t ix=0;ix<nBinsX;ix++)
417   {                                               424   {
418     file >> myDouble;                             425     file >> myDouble;
419     (*theMatrix)[ie*(fNBinsX+1)+ix] = myDouble << 426     (*theMatrix)[ie*(nBinsX+1)+ix] = myDouble*millibarn;
420   }                                               427   }
421       file >> myDouble; //total cross section     428       file >> myDouble; //total cross section
422       (*theMatrix)[ie*(fNBinsX+1)+fNBinsX] = m << 429       (*theMatrix)[ie*(nBinsX+1)+nBinsX] = myDouble*millibarn;
423     }                                             430     }
424                                                   431 
425   if (fElementData)                            << 432   if (theElementData)
426     fElementData->insert(std::make_pair(Z,theM << 433     theElementData->insert(std::make_pair(Z,theMatrix));
427   else                                            434   else
428     delete theMatrix;                             435     delete theMatrix;
429   file.close();                                   436   file.close();
430   return;                                         437   return;
431 }                                                 438 }
432                                                << 
433 //....oooOO0OOooo........oooOO0OOooo........oo    439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
434                                                   440 
435 G4double G4PenelopeBremsstrahlungFS::GetMoment    441 G4double G4PenelopeBremsstrahlungFS::GetMomentumIntegral(G4double* y,
436                  G4double xup,G4int momOrder)     442                  G4double xup,G4int momOrder) const
437 //x is always the gridX                           443 //x is always the gridX
438 {                                                 444 {
439   //Corresponds to the function RLMOM of Penel    445   //Corresponds to the function RLMOM of Penelope
440   //This method performs the calculation of th    446   //This method performs the calculation of the integral of (x^momOrder)*y over the interval
441   //from x[0] to xup, obtained by linear inter    447   //from x[0] to xup, obtained by linear interpolation on a table of y.
442   //The independent variable is assumed to tak    448   //The independent variable is assumed to take positive values only.
443   //                                              449   //
444   std::size_t size = fNBinsX;                  << 450   size_t size = nBinsX;
445   const G4double eps = 1e-35;                     451   const G4double eps = 1e-35;
446                                                   452 
447   //Check that the call is valid                  453   //Check that the call is valid
448   if (momOrder<-1 || size<2 || theXGrid[0]<0)     454   if (momOrder<-1 || size<2 || theXGrid[0]<0)
449     {                                             455     {
450       G4Exception("G4PenelopeBremsstrahlungFS:    456       G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
451       "em2011",FatalException,"Invalid call");    457       "em2011",FatalException,"Invalid call");
452     }                                             458     }
453                                                   459 
454   for (std::size_t i=1;i<size;++i)             << 460   for (size_t i=1;i<size;i++)
455     {                                             461     {
456       if (theXGrid[i]<0 || theXGrid[i]<theXGri    462       if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1])
457   {                                               463   {
458     G4ExceptionDescription ed;                    464     G4ExceptionDescription ed;
459     ed << "Invalid call for bin " << i << G4en    465     ed << "Invalid call for bin " << i << G4endl;
460     G4Exception("G4PenelopeBremsstrahlungFS::G    466     G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
461       "em2012",FatalException,ed);                467       "em2012",FatalException,ed);
462   }                                               468   }
463     }                                             469     }
464                                                   470 
465   //Compute the integral                          471   //Compute the integral
466   G4double result = 0;                            472   G4double result = 0;
467   if (xup < theXGrid[0])                          473   if (xup < theXGrid[0])
468     return result;                                474     return result;
469   G4bool loopAgain = true;                     << 475   bool loopAgain = true;
470   G4double xt = std::min(xup,theXGrid[size-1])    476   G4double xt = std::min(xup,theXGrid[size-1]);
471   G4double xtc = 0;                               477   G4double xtc = 0;
472   for (std::size_t i=0;i<size-1;++i)           << 478   for (size_t i=0;i<size-1;i++)
473     {                                             479     {
474       G4double x1 = std::max(theXGrid[i],eps);    480       G4double x1 = std::max(theXGrid[i],eps);
475       G4double y1 = y[i];                         481       G4double y1 = y[i];
476       G4double x2 = std::max(theXGrid[i+1],eps    482       G4double x2 = std::max(theXGrid[i+1],eps);
477       G4double y2 = y[i+1];                       483       G4double y2 = y[i+1];
478       if (xt < x2)                                484       if (xt < x2)
479   {                                               485   {
480     xtc = xt;                                     486     xtc = xt;
481     loopAgain = false;                            487     loopAgain = false;
482   }                                               488   }
483       else                                        489       else
484   xtc = x2;                                       490   xtc = x2;
485       G4double dx = x2-x1;                        491       G4double dx = x2-x1;
486       G4double dy = y2-y1;                        492       G4double dy = y2-y1;
487       G4double ds = 0;                            493       G4double ds = 0;
488       if (std::fabs(dx)>1e-14*std::fabs(dy))      494       if (std::fabs(dx)>1e-14*std::fabs(dy))
489   {                                               495   {
490     G4double b=dy/dx;                             496     G4double b=dy/dx;
491     G4double a=y1-b*x1;                           497     G4double a=y1-b*x1;
492     if (momOrder == -1)                           498     if (momOrder == -1)
493       ds = a*G4Log(xtc/x1)+b*(xtc-x1);         << 499       ds = a*std::log(xtc/x1)+b*(xtc-x1);
494     else if (momOrder == 0) //speed it up, not    500     else if (momOrder == 0) //speed it up, not using pow()
495       ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1);    501       ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1);
496     else                                          502     else
497       ds = a*(std::pow(xtc,momOrder+1)-std::po    503       ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((G4double) (momOrder + 1))
498         + b*(std::pow(xtc,momOrder+2)-std::pow    504         + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((G4double) (momOrder + 2));
499   }                                               505   }
500       else                                        506       else
501   ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOr    507   ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOrder);
502       result += ds;                               508       result += ds;
503       if (!loopAgain)                             509       if (!loopAgain)
504   return result;                                  510   return result;
505     }                                             511     }
506   return result;                                  512   return result;
507 }                                                 513 }
508                                                   514 
509 //....oooOO0OOooo........oooOO0OOooo........oo    515 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
510                                                   516 
511 const G4PhysicsTable* G4PenelopeBremsstrahlung    517 const G4PhysicsTable* G4PenelopeBremsstrahlungFS::GetScaledXSTable(const G4Material* mat,
512                    const G4double cut) const      518                    const G4double cut) const
513 {                                                 519 {
514   //check if it already contains the entry        520   //check if it already contains the entry
515   std::pair<const G4Material*,G4double> theKey    521   std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
516                                                   522 
517   if (!(fReducedXSTable->count(theKey)))       << 523   if (!(theReducedXSTable->count(theKey)))
518     {                                             524     {
519       G4Exception("G4PenelopeBremsstrahlungFS:    525       G4Exception("G4PenelopeBremsstrahlungFS::GetScaledXSTable()",
520       "em2013",FatalException,"Unable to retri    526       "em2013",FatalException,"Unable to retrieve the cross section table");
521     }                                             527     }
522                                                   528 
523   return fReducedXSTable->find(theKey)->second << 529   return theReducedXSTable->find(theKey)->second;
524 }                                                 530 }
525                                                   531 
526 //....oooOO0OOooo........oooOO0OOooo........oo    532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
527                                                   533 
528 void G4PenelopeBremsstrahlungFS::InitializeEne    534 void G4PenelopeBremsstrahlungFS::InitializeEnergySampling(const G4Material* material,
529                   G4double cut)                   535                   G4double cut)
530 {                                                 536 {
531   if (fVerbosity > 2)                             537   if (fVerbosity > 2)
532     G4cout << "Entering in G4PenelopeBremsstra    538     G4cout << "Entering in G4PenelopeBremsstrahlungFS::InitializeEnergySampling() for " <<
533       material->GetName() << G4endl;              539       material->GetName() << G4endl;
534                                                   540 
535   //This method should be accessed by the mast    541   //This method should be accessed by the master only
536   std::pair<const G4Material*,G4double> theKey    542   std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
537                                                   543 
538   G4PhysicsTable* thePhysicsTable = new G4Phys    544   G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
539   // the table will contain 57 G4PhysicsFreeVe    545   // the table will contain 57 G4PhysicsFreeVectors with different
540   // values of E.                                 546   // values of E.
541   G4PhysicsFreeVector* thePBvec = new G4Physic << 547   G4PhysicsFreeVector* thePBvec = new G4PhysicsFreeVector(nBinsE);
542                                                   548 
543   //I reserve space of the vectors.               549   //I reserve space of the vectors.
544   for (std::size_t i=0;i<fNBinsE;++i)          << 550   for (size_t i=0;i<nBinsE;i++)
545     thePhysicsTable->push_back(new G4PhysicsFr << 551     thePhysicsTable->push_back(new G4PhysicsFreeVector(nBinsX));
546                                                   552 
547   //Retrieve the table. Must already exist at     553   //Retrieve the table. Must already exist at this point, because this
548   //method is invoked by GetScaledXSTable()       554   //method is invoked by GetScaledXSTable()
549   if (!(fReducedXSTable->count(theKey)))       << 555   if (!(theReducedXSTable->count(theKey)))
550     G4Exception("G4PenelopeBremsstrahlungFS::I    556     G4Exception("G4PenelopeBremsstrahlungFS::InitializeEnergySampling()",
551     "em2013",FatalException,"Unable to retriev    557     "em2013",FatalException,"Unable to retrieve the cross section table");
552   G4PhysicsTable* theTableReduced = fReducedXS << 558   G4PhysicsTable* theTableReduced = theReducedXSTable->find(theKey)->second;
553                                                   559 
554   for (std::size_t ie=0;ie<fNBinsE;++ie)       << 560   for (size_t ie=0;ie<nBinsE;ie++)
555     {                                             561     {
556       G4PhysicsFreeVector* theVec =               562       G4PhysicsFreeVector* theVec =
557   (G4PhysicsFreeVector*) ((*thePhysicsTable)[i    563   (G4PhysicsFreeVector*) ((*thePhysicsTable)[ie]);
558       //Fill the table                            564       //Fill the table
559       G4double value = 0; //first value           565       G4double value = 0; //first value
560       theVec->PutValues(0,theXGrid[0],value);  << 566       theVec->PutValue(0,theXGrid[0],value);
561       for (std::size_t ix=1;ix<fNBinsX;++ix)   << 567       for (size_t ix=1;ix<nBinsX;ix++)
562   {                                               568   {
563     //Here calculate the cumulative distributi    569     //Here calculate the cumulative distribution
564     // int_{0}^{x} dSigma(x',E)/dx' (1/x') dx'    570     // int_{0}^{x} dSigma(x',E)/dx' (1/x') dx'
565     G4PhysicsFreeVector* v1 = (G4PhysicsFreeVe    571     G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableReduced)[ix-1];
566     G4PhysicsFreeVector* v2 = (G4PhysicsFreeVe    572     G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableReduced)[ix];
567                                                   573 
568     G4double x1=std::max(theXGrid[ix-1],1.0e-3    574     G4double x1=std::max(theXGrid[ix-1],1.0e-35);
569     //Remember: the table fReducedXSTable has  << 575     //Remember: the table theReducedXSTable has a fake first point in energy
570     //so, it contains one more bin than fNBins << 576     //so, it contains one more bin than nBinsE.
571     G4double y1=G4Exp((*v1)[ie+1]);               577     G4double y1=G4Exp((*v1)[ie+1]);
572     G4double x2=std::max(theXGrid[ix],1.0e-35)    578     G4double x2=std::max(theXGrid[ix],1.0e-35);
573     G4double y2=G4Exp((*v2)[ie+1]);               579     G4double y2=G4Exp((*v2)[ie+1]);
574     G4double B = (y2-y1)/(x2-x1);                 580     G4double B = (y2-y1)/(x2-x1);
575     G4double A = y1-B*x1;                         581     G4double A = y1-B*x1;
576     G4double dS = A*G4Log(x2/x1)+B*(x2-x1);    << 582     G4double dS = A*std::log(x2/x1)+B*(x2-x1);
577     value += dS;                                  583     value += dS;
578     theVec->PutValues(ix,theXGrid[ix],value);  << 584     theVec->PutValue(ix,theXGrid[ix],value);
579   }                                               585   }
                                                   >> 586 
580       //fill the PB vector                        587       //fill the PB vector
581       G4double xc = cut/theEGrid[ie];             588       G4double xc = cut/theEGrid[ie];
582       //Fill a temp data vector                   589       //Fill a temp data vector
583       G4double* tempData = new G4double[fNBins << 590       G4double* tempData = new G4double[nBinsX];
584       for (std::size_t ix=0;ix<fNBinsX;++ix)   << 591       for (size_t ix=0;ix<nBinsX;ix++)
585   {                                               592   {
586     G4PhysicsFreeVector* vv = (G4PhysicsFreeVe    593     G4PhysicsFreeVector* vv = (G4PhysicsFreeVector*) (*theTableReduced)[ix];
587     tempData[ix] = G4Exp((*vv)[ie+1]);            594     tempData[ix] = G4Exp((*vv)[ie+1]);
588   }                                               595   }
589       G4double pbval = (xc<=1) ?                  596       G4double pbval = (xc<=1) ?
590   GetMomentumIntegral(tempData,xc,-1) :           597   GetMomentumIntegral(tempData,xc,-1) :
591   GetMomentumIntegral(tempData,1.0,-1);           598   GetMomentumIntegral(tempData,1.0,-1);
592       thePBvec->PutValues(ie,theEGrid[ie],pbva << 599       thePBvec->PutValue(ie,theEGrid[ie],pbval);
593       delete[] tempData;                          600       delete[] tempData;
594     }                                             601     }
595                                                   602 
596   fSamplingTable->insert(std::make_pair(theKey << 603   theSamplingTable->insert(std::make_pair(theKey,thePhysicsTable));
597   fPBcut->insert(std::make_pair(theKey,thePBve << 604   thePBcut->insert(std::make_pair(theKey,thePBvec));
598   return;                                         605   return;
599 }                                                 606 }
600                                                   607 
601 //....oooOO0OOooo........oooOO0OOooo........oo    608 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
602                                                   609 
603 G4double G4PenelopeBremsstrahlungFS::SampleGam    610 G4double G4PenelopeBremsstrahlungFS::SampleGammaEnergy(G4double energy,const G4Material* mat,
604                    const G4double cut) const      611                    const G4double cut) const
605 {                                                 612 {
606   std::pair<const G4Material*,G4double> theKey    613   std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
607   if (!(fSamplingTable->count(theKey)) || !(fP << 614   if (!(theSamplingTable->count(theKey)) || !(thePBcut->count(theKey)))
608     {                                             615     {
609       G4ExceptionDescription ed;                  616       G4ExceptionDescription ed;
610       ed << "Unable to retrieve the SamplingTa    617       ed << "Unable to retrieve the SamplingTable: " <<
611   fSamplingTable->count(theKey) << " " <<      << 618   theSamplingTable->count(theKey) << " " <<
612   fPBcut->count(theKey) << G4endl;             << 619   thePBcut->count(theKey) << G4endl;
613       G4Exception("G4PenelopeBremsstrahlungFS:    620       G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
614       "em2014",FatalException,ed);                621       "em2014",FatalException,ed);
615     }                                             622     }
616   const G4PhysicsTable* theTableInte = fSampli << 623 
617   const G4PhysicsTable* theTableRed = fReduced << 624 
                                                   >> 625   const G4PhysicsTable* theTableInte = theSamplingTable->find(theKey)->second;
                                                   >> 626   const G4PhysicsTable* theTableRed = theReducedXSTable->find(theKey)->second;
618                                                   627 
619   //Find the energy bin using bi-partition        628   //Find the energy bin using bi-partition
620   std::size_t eBin = 0;                        << 629   size_t eBin = 0;
621   G4bool firstOrLastBin = false;                  630   G4bool firstOrLastBin = false;
622                                                   631 
623   if (energy < theEGrid[0]) //below first bin     632   if (energy < theEGrid[0]) //below first bin
624     {                                             633     {
625       eBin = 0;                                   634       eBin = 0;
626       firstOrLastBin = true;                      635       firstOrLastBin = true;
627     }                                             636     }
628   else if (energy > theEGrid[fNBinsE-1]) //aft << 637   else if (energy > theEGrid[nBinsE-1]) //after last bin
629     {                                             638     {
630       eBin = fNBinsE-1;                        << 639       eBin = nBinsE-1;
631       firstOrLastBin = true;                      640       firstOrLastBin = true;
632     }                                             641     }
633   else                                            642   else
634     {                                             643     {
635       std::size_t i=0;                         << 644       size_t i=0;
636       std::size_t j=fNBinsE-1;                 << 645       size_t j=nBinsE-1;
637       while ((j-i)>1)                             646       while ((j-i)>1)
638   {                                               647   {
639     std::size_t k = (i+j)/2;                   << 648     size_t k = (i+j)/2;
640     if (energy > theEGrid[k])                     649     if (energy > theEGrid[k])
641       i = k;                                      650       i = k;
642     else                                          651     else
643       j = k;                                      652       j = k;
644   }                                               653   }
645       eBin = i;                                   654       eBin = i;
646     }                                             655     }
647                                                   656 
648   //Get the appropriate physics vector            657   //Get the appropriate physics vector
649   const G4PhysicsFreeVector* theVec1 = (G4Phys    658   const G4PhysicsFreeVector* theVec1 = (G4PhysicsFreeVector*) (*theTableInte)[eBin];
650                                                   659 
651   //Use a "temporary" vector which contains th    660   //Use a "temporary" vector which contains the linear interpolation of the x spectra
652   //in energy. The temporary vector is thread-    661   //in energy. The temporary vector is thread-local, so that there is no conflict.
653   //This is achieved via G4Cache. The theTempV    662   //This is achieved via G4Cache. The theTempVect is allocated only once per thread
654   //(member variable), but it is overwritten a    663   //(member variable), but it is overwritten at every call of this method
655   //(because the interpolation factors change!    664   //(because the interpolation factors change!)
656   G4PhysicsFreeVector* theTempVec = fCache.Get    665   G4PhysicsFreeVector* theTempVec = fCache.Get();
657   if (!theTempVec) //First time this thread ge    666   if (!theTempVec) //First time this thread gets the cache
658     {                                             667     {
659       theTempVec = new G4PhysicsFreeVector(fNB << 668       theTempVec = new G4PhysicsFreeVector(nBinsX);   
660       fCache.Put(theTempVec);                     669       fCache.Put(theTempVec);
661       // The G4AutoDelete takes care here to c    670       // The G4AutoDelete takes care here to clean up the vectors
662       G4AutoDelete::Register(theTempVec);         671       G4AutoDelete::Register(theTempVec);
663       if (fVerbosity > 4)                         672       if (fVerbosity > 4)
664   G4cout << "Creating new instance of G4Physic    673   G4cout << "Creating new instance of G4PhysicsFreeVector() on the worker" << G4endl;
665     }                                             674     }
666                                                   675 
667   //theTempVect is allocated only once (member    676   //theTempVect is allocated only once (member variable), but it is overwritten at
668   //every call of this method (because the int    677   //every call of this method (because the interpolation factors change!)
669   if (!firstOrLastBin)                            678   if (!firstOrLastBin)
670     {                                             679     {
671       const G4PhysicsFreeVector* theVec2 = (G4    680       const G4PhysicsFreeVector* theVec2 = (G4PhysicsFreeVector*) (*theTableInte)[eBin+1];
672       for (std::size_t iloop=0;iloop<fNBinsX;+ << 681       for (size_t iloop=0;iloop<nBinsX;iloop++)
673   {                                               682   {
674     G4double val = (*theVec1)[iloop]+(((*theVe    683     G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))*
675       (energy-theEGrid[eBin])/(theEGrid[eBin+1    684       (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
676     theTempVec->PutValues(iloop,theXGrid[iloop << 685     theTempVec->PutValue(iloop,theXGrid[iloop],val);
677   }                                               686   }
678     }                                             687     }
679   else //first or last bin, no interpolation      688   else //first or last bin, no interpolation
680     {                                             689     {
681       for (std::size_t iloop=0;iloop<fNBinsX;+ << 690       for (size_t iloop=0;iloop<nBinsX;iloop++)
682   theTempVec->PutValues(iloop,theXGrid[iloop], << 691   theTempVec->PutValue(iloop,theXGrid[iloop],(*theVec1)[iloop]);
683     }                                             692     }
684                                                   693 
685   //Start the game                                694   //Start the game
686   G4double pbcut = (*(fPBcut->find(theKey)->se << 695   G4double pbcut = (*(thePBcut->find(theKey)->second))[eBin];
687                                                   696 
688   if (!firstOrLastBin) //linear interpolation     697   if (!firstOrLastBin) //linear interpolation on pbcut as well
689     {                                             698     {
690       pbcut = (*(fPBcut->find(theKey)->second) << 699       pbcut = (*(thePBcut->find(theKey)->second))[eBin] +
691   ((*(fPBcut->find(theKey)->second))[eBin+1]-( << 700   ((*(thePBcut->find(theKey)->second))[eBin+1]-(*(thePBcut->find(theKey)->second))[eBin])*
692   (energy-theEGrid[eBin])/(theEGrid[eBin+1]-th    701   (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
693     }                                             702     }
694                                                   703 
695   G4double pCumulative = (*theTempVec)[fNBinsX << 704   G4double pCumulative = (*theTempVec)[nBinsX-1]; //last value
696                                                   705 
697   G4double eGamma = 0;                            706   G4double eGamma = 0;
698   G4int nIterations = 0;                          707   G4int nIterations = 0;
699   do                                              708   do
700     {                                             709     {
701       G4double pt = pbcut + G4UniformRand()*(p    710       G4double pt = pbcut + G4UniformRand()*(pCumulative - pbcut);
702       nIterations++;                              711       nIterations++;
703                                                   712 
704       //find where it is                          713       //find where it is
705       std::size_t ibin = 0;                    << 714       size_t ibin = 0;
706       if (pt < (*theTempVec)[0])                  715       if (pt < (*theTempVec)[0])
707   ibin = 0;                                       716   ibin = 0;
708       else if (pt > (*theTempVec)[fNBinsX-1])  << 717       else if (pt > (*theTempVec)[nBinsX-1])
709   {                                               718   {
710     //We observed problems due to numerical ro    719     //We observed problems due to numerical rounding here (STT).
711     //delta here is a tiny positive number        720     //delta here is a tiny positive number
712     G4double delta = pt-(*theTempVec)[fNBinsX- << 721     G4double delta = pt-(*theTempVec)[nBinsX-1];
713     if (delta < pt*1e-10) // very small! Numer    722     if (delta < pt*1e-10) // very small! Numerical rounding only
714       {                                           723       {
715         ibin = fNBinsX-2;                      << 724         ibin = nBinsX-2;
716         G4ExceptionDescription ed;                725         G4ExceptionDescription ed;
717         ed << "Found that (pt > (*theTempVec)[ << 726         ed << "Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
718     " , (*theTempVec)[fNBinsX-1] = " << (*theT << 727     " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[nBinsX-1] << " and delta = " <<
719     (pt-(*theTempVec)[fNBinsX-1]) << G4endl;   << 728     (pt-(*theTempVec)[nBinsX-1]) << G4endl;
720         ed << "Possible symptom of problem wit    729         ed << "Possible symptom of problem with numerical precision" << G4endl;
721         G4Exception("G4PenelopeBremsstrahlungF    730         G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
722         "em2015",JustWarning,ed);                 731         "em2015",JustWarning,ed);
723       }                                           732       }
724     else //real problem                           733     else //real problem
725       {                                           734       {
726         G4ExceptionDescription ed;                735         G4ExceptionDescription ed;
727         ed << "Crash at (pt > (*theTempVec)[fN << 736         ed << "Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
728     " , (*theTempVec)[fNBinsX-1]=" << (*theTem << 737     " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[nBinsX-1] << " and nBinsX = " <<
729     fNBinsX << G4endl;                         << 738     nBinsX << G4endl;
730         ed << "Material: " << mat->GetName() <    739         ed << "Material: " << mat->GetName() << ", energy = " << energy/keV << " keV" <<
731     G4endl;                                       740     G4endl;
732         G4Exception("G4PenelopeBremsstrahlungF    741         G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
733         "em2015",FatalException,ed);              742         "em2015",FatalException,ed);
734       }                                           743       }
735   }                                               744   }
736       else                                        745       else
737   {                                               746   {
738     std::size_t i=0;                           << 747     size_t i=0;
739     std::size_t j=fNBinsX-1;                   << 748     size_t j=nBinsX-1;
740     while ((j-i)>1)                               749     while ((j-i)>1)
741       {                                           750       {
742         std::size_t k = (i+j)/2;               << 751         size_t k = (i+j)/2;
743         if (pt > (*theTempVec)[k])                752         if (pt > (*theTempVec)[k])
744     i = k;                                        753     i = k;
745         else                                      754         else
746     j = k;                                        755     j = k;
747       }                                           756       }
748     ibin = i;                                     757     ibin = i;
749   }                                               758   }
750                                                   759 
751       G4double w1 = theXGrid[ibin];               760       G4double w1 = theXGrid[ibin];
752       G4double w2 = theXGrid[ibin+1];             761       G4double w2 = theXGrid[ibin+1];
753                                                   762 
754       const G4PhysicsFreeVector* v1 = (G4Physi    763       const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableRed)[ibin];
755       const G4PhysicsFreeVector* v2 = (G4Physi    764       const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableRed)[ibin+1];
756       //Remember: the table fReducedXSTable ha << 765       //Remember: the table theReducedXSTable has a fake first point in energy
757       //so, it contains one more bin than fNBi << 766       //so, it contains one more bin than nBinsE.
758       G4double pdf1 = G4Exp((*v1)[eBin+1]);       767       G4double pdf1 = G4Exp((*v1)[eBin+1]);
759       G4double pdf2 = G4Exp((*v2)[eBin+1]);       768       G4double pdf2 = G4Exp((*v2)[eBin+1]);
760       G4double deltaW = w2-w1;                    769       G4double deltaW = w2-w1;
761       G4double dpdfb = pdf2-pdf1;                 770       G4double dpdfb = pdf2-pdf1;
762       G4double B = dpdfb/deltaW;                  771       G4double B = dpdfb/deltaW;
763       G4double A = pdf1-B*w1;                     772       G4double A = pdf1-B*w1;
764       //I already made an interpolation in ene    773       //I already made an interpolation in energy, so I can use the actual value for the
765       //calculation of the wbcut, instead of t    774       //calculation of the wbcut, instead of the grid values (except for the last bin)
766       G4double wbcut  = (cut < energy) ? cut/e    775       G4double wbcut  = (cut < energy) ? cut/energy : 1.0;
767       if (firstOrLastBin) //this is an particu    776       if (firstOrLastBin) //this is an particular case: no interpolation available
768   wbcut  = (cut < theEGrid[eBin]) ? cut/theEGr    777   wbcut  = (cut < theEGrid[eBin]) ? cut/theEGrid[eBin] : 1.0;
769                                                   778 
770       if (w1 < wbcut)                             779       if (w1 < wbcut)
771   w1 = wbcut;                                     780   w1 = wbcut;
772       if (w2 < w1)                                781       if (w2 < w1)
773   {                                               782   {
774     //This configuration can happen if initial    783     //This configuration can happen if initially wbcut > w2 > w1. Due to the previous
775     //statement, (w1 = wbcut), it becomes wbcu    784     //statement, (w1 = wbcut), it becomes wbcut = w1 > w2. In this case, it is not a
776     //real problem. It becomes a problem if w2    785     //real problem. It becomes a problem if w2 < w1 before the w1 = wbcut statement. Issue
777     //a warning only in this specific case.       786     //a warning only in this specific case.
778     if (w2 > wbcut)                               787     if (w2 > wbcut)
779       {                                           788       {
780         G4ExceptionDescription ed;                789         G4ExceptionDescription ed;
781         ed << "Warning in G4PenelopeBremsstrah    790         ed << "Warning in G4PenelopeBremsstrahlungFS::SampleX()" << G4endl;
782         ed << "Conflicting end-point values: w    791         ed << "Conflicting end-point values: w1=" << w1 << "; w2 = " << w2 << G4endl;
783         ed << "wbcut = " << wbcut << " energy=    792         ed << "wbcut = " << wbcut << " energy= " << energy/keV << " keV" << G4endl;
784         ed << "cut = " << cut/keV << " keV" <<    793         ed << "cut = " << cut/keV << " keV" << G4endl;
785         G4Exception("G4PenelopeBremsstrahlungF    794         G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()","em2015",
786         JustWarning,ed);                          795         JustWarning,ed);
787       }                                           796       }
788     return w1*energy;                             797     return w1*energy;
789   }                                               798   }
790                                                   799 
791       G4double pmax = std::max(A+B*w1,A+B*w2);    800       G4double pmax = std::max(A+B*w1,A+B*w2);
792       G4bool loopAgain = false;                   801       G4bool loopAgain = false;
793       do                                          802       do
794   {                                               803   {
795     loopAgain = false;                            804     loopAgain = false;
796     eGamma = w1* std::pow((w2/w1),G4UniformRan    805     eGamma = w1* std::pow((w2/w1),G4UniformRand());
797     if  (G4UniformRand()*pmax > (A+B*eGamma))     806     if  (G4UniformRand()*pmax > (A+B*eGamma))
798       loopAgain = true;                           807       loopAgain = true;
799   }while(loopAgain);                              808   }while(loopAgain);
800       eGamma *= energy;                           809       eGamma *= energy;
801       if (nIterations > 100) //protection agai    810       if (nIterations > 100) //protection against infinite loops
802   return eGamma;                                  811   return eGamma;
803     }while(eGamma < cut); //repeat if sampled     812     }while(eGamma < cut); //repeat if sampled sub-cut!
804                                                   813 
805   return eGamma;                                  814   return eGamma;
806 }                                                 815 }
807                                                   816