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


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