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 9.6)


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