Geant4 Cross Reference

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


  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: G4PenelopeRayleighModel.cc,v 1.4 2010-12-15 10:26:41 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 // 03 Dec 2009   L Pandola    First implementa     33 // 03 Dec 2009   L Pandola    First implementation
 32 // 25 May 2011   L.Pandola    Renamed (make v2     34 // 25 May 2011   L.Pandola    Renamed (make v2008 as default Penelope)
 33 // 19 Sep 2013   L.Pandola    Migration to MT  << 
 34 //                                                 35 //
 35                                                    36 
 36 #include "G4PenelopeRayleighModel.hh"              37 #include "G4PenelopeRayleighModel.hh"
 37 #include "G4PhysicalConstants.hh"              << 
 38 #include "G4SystemOfUnits.hh"                  << 
 39 #include "G4PenelopeSamplingData.hh"               38 #include "G4PenelopeSamplingData.hh"
 40 #include "G4ParticleDefinition.hh"                 39 #include "G4ParticleDefinition.hh"
 41 #include "G4MaterialCutsCouple.hh"                 40 #include "G4MaterialCutsCouple.hh"
 42 #include "G4ProductionCutsTable.hh"                41 #include "G4ProductionCutsTable.hh"
 43 #include "G4DynamicParticle.hh"                    42 #include "G4DynamicParticle.hh"
 44 #include "G4PhysicsTable.hh"                       43 #include "G4PhysicsTable.hh"
 45 #include "G4ElementTable.hh"                       44 #include "G4ElementTable.hh"
 46 #include "G4Element.hh"                            45 #include "G4Element.hh"
 47 #include "G4PhysicsFreeVector.hh"                  46 #include "G4PhysicsFreeVector.hh"
 48 #include "G4AutoLock.hh"                       << 
 49 #include "G4Exp.hh"                            << 
 50                                                    47 
 51 //....oooOO0OOooo........oooOO0OOooo........oo << 
 52                                                << 
 53 const G4int G4PenelopeRayleighModel::fMaxZ;    << 
 54 G4PhysicsFreeVector* G4PenelopeRayleighModel:: << 
 55 G4PhysicsFreeVector* G4PenelopeRayleighModel:: << 
 56                                                    48 
 57 //....oooOO0OOooo........oooOO0OOooo........oo     49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 58                                                    50 
 59 G4PenelopeRayleighModel::G4PenelopeRayleighMod <<  51 
                                                   >>  52 G4PenelopeRayleighModel::G4PenelopeRayleighModel(const G4ParticleDefinition*,
 60              const G4String& nam)                  53              const G4String& nam)
 61   :G4VEmModel(nam),fParticleChange(nullptr),fP <<  54   :G4VEmModel(nam),isInitialised(false),logAtomicCrossSection(0),   
 62    fLogFormFactorTable(nullptr),fPMaxTable(nul <<  55    atomicFormFactor(0),logFormFactorTable(0),pMaxTable(0),samplingTable(0)
 63    fIsInitialised(false),fLocalTable(false)    << 
 64 {                                                  56 {
 65   fIntrinsicLowEnergyLimit = 100.0*eV;             57   fIntrinsicLowEnergyLimit = 100.0*eV;
 66   fIntrinsicHighEnergyLimit = 100.0*GeV;           58   fIntrinsicHighEnergyLimit = 100.0*GeV;
 67   //  SetLowEnergyLimit(fIntrinsicLowEnergyLim     59   //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
 68   SetHighEnergyLimit(fIntrinsicHighEnergyLimit     60   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
 69                                                << 
 70   if (part)                                    << 
 71     SetParticle(part);                         << 
 72                                                << 
 73   //                                               61   //
 74   fVerboseLevel= 0;                            <<  62   verboseLevel= 0;
 75   // Verbosity scale:                              63   // Verbosity scale:
 76   // 0 = nothing                               <<  64   // 0 = nothing 
 77   // 1 = warning for energy non-conservation   <<  65   // 1 = warning for energy non-conservation 
 78   // 2 = details of energy budget                  66   // 2 = details of energy budget
 79   // 3 = calculation of cross sections, file o     67   // 3 = calculation of cross sections, file openings, sampling of atoms
 80   // 4 = entering in methods                   <<  68   // 4 = entering in methods 
 81                                                    69 
 82   //build the energy grid. It is the same for      70   //build the energy grid. It is the same for all materials
 83   G4double logenergy = G4Log(fIntrinsicLowEner <<  71   G4double logenergy = std::log(fIntrinsicLowEnergyLimit/2.);
 84   G4double logmaxenergy = G4Log(1.5*fIntrinsic <<  72   G4double logmaxenergy = std::log(1.5*fIntrinsicHighEnergyLimit);
 85   //finer grid below 160 keV                       73   //finer grid below 160 keV
 86   G4double logtransitionenergy = G4Log(160*keV <<  74   G4double logtransitionenergy = std::log(160*keV); 
 87   G4double logfactor1 = G4Log(10.)/250.;       <<  75   G4double logfactor1 = std::log(10.)/250.;
 88   G4double logfactor2 = logfactor1*10;             76   G4double logfactor2 = logfactor1*10;
 89   fLogEnergyGridPMax.push_back(logenergy);     <<  77   logEnergyGridPMax.push_back(logenergy);
 90   do{                                              78   do{
 91     if (logenergy < logtransitionenergy)           79     if (logenergy < logtransitionenergy)
 92       logenergy += logfactor1;                     80       logenergy += logfactor1;
 93     else                                           81     else
 94       logenergy += logfactor2;                     82       logenergy += logfactor2;
 95     fLogEnergyGridPMax.push_back(logenergy);   <<  83     logEnergyGridPMax.push_back(logenergy);      
 96   }while (logenergy < logmaxenergy);               84   }while (logenergy < logmaxenergy);
 97 }                                                  85 }
 98                                                    86 
 99 //....oooOO0OOooo........oooOO0OOooo........oo     87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100                                                    88 
101 G4PenelopeRayleighModel::~G4PenelopeRayleighMo     89 G4PenelopeRayleighModel::~G4PenelopeRayleighModel()
102 {                                                  90 {
103   if (IsMaster() || fLocalTable)               <<  91   std::map <const G4int,G4PhysicsFreeVector*>::iterator i;
                                                   >>  92   if (logAtomicCrossSection)
104     {                                              93     {
105                                                <<  94       for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
106       for(G4int i=0; i<=fMaxZ; ++i)            <<  95   if (i->second) delete i->second;
107   {                                            <<  96       delete logAtomicCrossSection;
108     if(fLogAtomicCrossSection[i])              <<  97      }
109       {                                        <<  98 
110         delete fLogAtomicCrossSection[i];      <<  99    if (atomicFormFactor)
111         fLogAtomicCrossSection[i] = nullptr;   << 100      {
112       }                                        << 101        for (i=atomicFormFactor->begin();i != atomicFormFactor->end();i++)
113     if(fAtomicFormFactor[i])                   << 102    if (i->second) delete i->second;
114       {                                        << 103        delete atomicFormFactor;
115         delete fAtomicFormFactor[i];           << 104      }
116         fAtomicFormFactor[i] = nullptr;        << 105 
117       }                                        << 106   ClearTables();
118   }                                            << 
119       ClearTables();                           << 
120     }                                          << 
121 }                                                 107 }
122                                                   108 
123 //....oooOO0OOooo........oooOO0OOooo........oo    109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 void G4PenelopeRayleighModel::ClearTables()       110 void G4PenelopeRayleighModel::ClearTables()
125 {                                                 111 {
126    if (fLogFormFactorTable)                    << 112    std::map <const G4Material*,G4PhysicsFreeVector*>::iterator i;
                                                   >> 113  
                                                   >> 114    if (logFormFactorTable)
127      {                                            115      {
128        for (auto& item : (*fLogFormFactorTable << 116        for (i=logFormFactorTable->begin(); i != logFormFactorTable->end(); i++)
129    if (item.second) delete item.second;        << 117    if (i->second) delete i->second;
130        delete fLogFormFactorTable;             << 118        delete logFormFactorTable;
131        fLogFormFactorTable = nullptr; //zero e << 119        logFormFactorTable = 0; //zero explicitely
132      }                                            120      }
133    if (fPMaxTable)                             << 121 
                                                   >> 122    if (pMaxTable)
134      {                                            123      {
135        for (auto& item : (*fPMaxTable))        << 124        for (i=pMaxTable->begin(); i != pMaxTable->end(); i++)
136    if (item.second) delete item.second;        << 125    if (i->second) delete i->second;
137        delete fPMaxTable;                      << 126        delete pMaxTable;
138        fPMaxTable = nullptr; //zero explicitly << 127        pMaxTable = 0; //zero explicitely
139      }                                            128      }
140    if (fSamplingTable)                         << 129 
                                                   >> 130    std::map<const G4Material*,G4PenelopeSamplingData*>::iterator ii;
                                                   >> 131    if (samplingTable)
141      {                                            132      {
142        for (auto& item : (*fSamplingTable))    << 133        for (ii=samplingTable->begin(); ii != samplingTable->end(); ii++)
143    if (item.second) delete item.second;        << 134    if (ii->second) delete ii->second;
144        delete fSamplingTable;                  << 135        delete samplingTable;
145        fSamplingTable = nullptr; //zero explic << 136        samplingTable = 0; //zero explicitely
146      }                                         << 137      }     
                                                   >> 138 
147    return;                                        139    return;
148 }                                                 140 }
149                                                   141 
150 //....oooOO0OOooo........oooOO0OOooo........oo    142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151                                                   143 
152 void G4PenelopeRayleighModel::Initialise(const << 144 void G4PenelopeRayleighModel::Initialise(const G4ParticleDefinition* ,
153            const G4DataVector& )                  145            const G4DataVector& )
154 {                                                 146 {
155   if (fVerboseLevel > 3)                       << 147   if (verboseLevel > 3)
156     G4cout << "Calling G4PenelopeRayleighModel    148     G4cout << "Calling G4PenelopeRayleighModel::Initialise()" << G4endl;
157                                                   149 
158   SetParticle(part);                           << 150   //clear tables depending on materials, not the atomic ones
159                                                << 151   ClearTables();
160   //Only the master model creates/fills/destro << 152   
161   if (IsMaster() && part == fParticle)         << 153   //create new tables
162     {                                          << 154   //
163       //clear tables depending on materials, n << 155   // logAtomicCrossSection and atomicFormFactor are created only once,
164       ClearTables();                           << 156   // since they are never cleared
165                                                << 157   if (!logAtomicCrossSection)
166       if (fVerboseLevel > 3)                   << 158     logAtomicCrossSection = new std::map<const G4int,G4PhysicsFreeVector*>;
167   G4cout << "Calling G4PenelopeRayleighModel:: << 159   if (!atomicFormFactor)
168                                                << 160     atomicFormFactor = new std::map<const G4int,G4PhysicsFreeVector*>;
169       //create new tables                      << 161 
170       if (!fLogFormFactorTable)                << 162   if (!logFormFactorTable)
171   fLogFormFactorTable = new std::map<const G4M << 163     logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
172       if (!fPMaxTable)                         << 164   if (!pMaxTable)
173   fPMaxTable = new std::map<const G4Material*, << 165     pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
174       if (!fSamplingTable)                     << 166   if (!samplingTable)
175   fSamplingTable = new std::map<const G4Materi << 167     samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
176                                                << 168 
177       G4ProductionCutsTable* theCoupleTable =  << 169 
178   G4ProductionCutsTable::GetProductionCutsTabl << 170   if (verboseLevel > 0) {
179                                                << 171     G4cout << "Penelope Rayleigh model v2008 is initialized " << G4endl
180       for (G4int i=0;i<(G4int)theCoupleTable-> << 172      << "Energy range: "
181   {                                            << 173      << LowEnergyLimit() / keV << " keV - "
182     const G4Material* material =               << 174      << HighEnergyLimit() / GeV << " GeV"
183       theCoupleTable->GetMaterialCutsCouple(i) << 175      << G4endl;
184     const G4ElementVector* theElementVector =  << 176   }
185                                                << 
186     for (std::size_t j=0;j<material->GetNumber << 
187       {                                        << 
188         G4int iZ = theElementVector->at(j)->Ge << 
189         //read data files only in the master   << 
190         if (!fLogAtomicCrossSection[iZ])       << 
191     ReadDataFile(iZ);                          << 
192       }                                        << 
193                                                << 
194     //1) If the table has not been built for t << 
195     if (!fLogFormFactorTable->count(material)) << 
196       BuildFormFactorTable(material);          << 
197                                                << 
198     //2) retrieve or build the sampling table  << 
199     if (!(fSamplingTable->count(material)))    << 
200       InitializeSamplingAlgorithm(material);   << 
201                                                << 
202     //3) retrieve or build the pMax data       << 
203     if (!fPMaxTable->count(material))          << 
204       GetPMaxTable(material);                  << 
205   }                                            << 
206                                                << 
207       if (fVerboseLevel > 1) {                 << 
208   G4cout << "Penelope Rayleigh model v2008 is  << 
209          << "Energy range: "                   << 
210          << LowEnergyLimit() / keV << " keV -  << 
211          << HighEnergyLimit() / GeV << " GeV"  << 
212          << G4endl;                            << 
213       }                                        << 
214     }                                          << 
215                                                   177 
216   if(fIsInitialised) return;                   << 178   if(isInitialised) return;
217   fParticleChange = GetParticleChangeForGamma(    179   fParticleChange = GetParticleChangeForGamma();
218   fIsInitialised = true;                       << 180   isInitialised = true;
219 }                                                 181 }
220                                                   182 
221 //....oooOO0OOooo........oooOO0OOooo........oo    183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222                                                   184 
223 void G4PenelopeRayleighModel::InitialiseLocal( << 
224                  G4VEmModel *masterModel)      << 
225 {                                              << 
226   if (fVerboseLevel > 3)                       << 
227     G4cout << "Calling  G4PenelopeRayleighMode << 
228   //                                           << 
229   //Check that particle matches: one might hav << 
230   //for e+ and e-).                            << 
231   //                                           << 
232   if (part == fParticle)                       << 
233     {                                          << 
234       //Get the const table pointers from the  << 
235       const G4PenelopeRayleighModel* theModel  << 
236   static_cast<G4PenelopeRayleighModel*> (maste << 
237                                                << 
238       //Copy pointers to the data tables       << 
239       for(G4int i=0; i<=fMaxZ; ++i)            << 
240   {                                            << 
241     fLogAtomicCrossSection[i] = theModel->fLog << 
242     fAtomicFormFactor[i] = theModel->fAtomicFo << 
243   }                                            << 
244       fLogFormFactorTable = theModel->fLogForm << 
245       fPMaxTable = theModel->fPMaxTable;       << 
246       fSamplingTable = theModel->fSamplingTabl << 
247                                                << 
248       //copy the G4DataVector with the grid    << 
249       fLogQSquareGrid = theModel->fLogQSquareG << 
250                                                << 
251       //Same verbosity for all workers, as the << 
252       fVerboseLevel = theModel->fVerboseLevel; << 
253     }                                          << 
254                                                << 
255   return;                                      << 
256 }                                              << 
257                                                << 
258                                                << 
259 //....oooOO0OOooo........oooOO0OOooo........oo << 
260 namespace { G4Mutex  PenelopeRayleighModelMute << 
261 G4double G4PenelopeRayleighModel::ComputeCross    185 G4double G4PenelopeRayleighModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
262                    G4double energy,               186                    G4double energy,
263                    G4double Z,                    187                    G4double Z,
264                    G4double,                      188                    G4double,
265                    G4double,                      189                    G4double,
266                    G4double)                      190                    G4double)
267 {                                                 191 {
268   // Cross section of Rayleigh scattering in P << 192   // Cross section of Rayleigh scattering in Penelope v2008 is calculated by the EPDL97 
269   // tabulation, Cuellen et al. (1997), with n << 193   // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel 
270   // et al. J. Phys. Chem. Ref. Data 4 (1975)     194   // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615.
271                                                   195 
272    if (fVerboseLevel > 3)                      << 196    if (verboseLevel > 3)
273     G4cout << "Calling CrossSectionPerAtom() o    197     G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" << G4endl;
                                                   >> 198  
                                                   >> 199    G4int iZ = (G4int) Z;
274                                                   200 
275    G4int iZ = G4int(Z);                        << 201    //read data files
276                                                << 202    if (!logAtomicCrossSection->count(iZ))
277    if (!fLogAtomicCrossSection[iZ])            << 203      ReadDataFile(iZ);
                                                   >> 204    //now it should be ok
                                                   >> 205    if (!logAtomicCrossSection->count(iZ))
278      {                                            206      {
279        //If we are here, it means that Initial << 207        G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
280        //not filled up. This can happen in a U << 208        "em2040",FatalException,"Unable to load the cross section table");
281        if (fVerboseLevel > 0)                  << 
282   {                                            << 
283     //Issue a G4Exception (warning) only in ve << 
284     G4ExceptionDescription ed;                 << 
285     ed << "Unable to retrieve the cross sectio << 
286     ed << "This can happen only in Unit Tests  << 
287     G4Exception("G4PenelopeRayleighModel::Comp << 
288           "em2040",JustWarning,ed);            << 
289   }                                            << 
290        //protect file reading via autolock     << 
291        G4AutoLock lock(&PenelopeRayleighModelM << 
292        ReadDataFile(iZ);                       << 
293        lock.unlock();                          << 
294      }                                            209      }
295                                                   210 
296    G4double cross = 0;                            211    G4double cross = 0;
297    G4PhysicsFreeVector* atom = fLogAtomicCross << 212 
                                                   >> 213    G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second;
298    if (!atom)                                     214    if (!atom)
299      {                                            215      {
300        G4ExceptionDescription ed;                 216        G4ExceptionDescription ed;
301        ed << "Unable to find Z=" << iZ << " in    217        ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl;
302        G4Exception("G4PenelopeRayleighModel::C    218        G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
303        "em2041",FatalException,ed);               219        "em2041",FatalException,ed);
304        return 0;                                  220        return 0;
305      }                                            221      }
306    G4double logene = G4Log(energy);            << 222    G4double logene = std::log(energy);
307    G4double logXS = atom->Value(logene);          223    G4double logXS = atom->Value(logene);
308    cross = G4Exp(logXS);                       << 224    cross = std::exp(logXS);
309                                                   225 
310    if (fVerboseLevel > 2)                      << 226    if (verboseLevel > 2)
311      {                                         << 227     G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z << 
312        G4cout << "Rayleigh cross section at "  << 228       " = " << cross/barn << " barn" << G4endl;
313    " = " << cross/barn << " barn" << G4endl;   << 229     return cross;
314      }                                         << 
315    return cross;                               << 
316 }                                                 230 }
317                                                   231 
318                                                   232 
319 //....oooOO0OOooo........oooOO0OOooo........oo    233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
320 void G4PenelopeRayleighModel::BuildFormFactorT    234 void G4PenelopeRayleighModel::BuildFormFactorTable(const G4Material* material)
321 {                                                 235 {
                                                   >> 236  
322   /*                                              237   /*
323     1) get composition and equivalent molecula    238     1) get composition and equivalent molecular density
324   */                                              239   */
325   std::size_t nElements = material->GetNumberO << 240   
                                                   >> 241   G4int nElements = material->GetNumberOfElements();
326   const G4ElementVector* elementVector = mater    242   const G4ElementVector* elementVector = material->GetElementVector();
327   const G4double* fractionVector = material->G    243   const G4double* fractionVector = material->GetFractionVector();
328                                                   244 
329   std::vector<G4double> *StechiometricFactors     245   std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
330   for (std::size_t i=0;i<nElements;++i)        << 246   for (G4int i=0;i<nElements;i++)
331     {                                             247     {
332       G4double fraction = fractionVector[i];      248       G4double fraction = fractionVector[i];
333       G4double atomicWeigth = (*elementVector)    249       G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
334       StechiometricFactors->push_back(fraction    250       StechiometricFactors->push_back(fraction/atomicWeigth);
335     }                                             251     }
336   //Find max                                      252   //Find max
337   G4double MaxStechiometricFactor = 0.;           253   G4double MaxStechiometricFactor = 0.;
338   for (std::size_t i=0;i<nElements;++i)        << 254   for (G4int i=0;i<nElements;i++)
339     {                                             255     {
340       if ((*StechiometricFactors)[i] > MaxStec    256       if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
341         MaxStechiometricFactor = (*Stechiometr    257         MaxStechiometricFactor = (*StechiometricFactors)[i];
342     }                                             258     }
343   if (MaxStechiometricFactor<1e-16)               259   if (MaxStechiometricFactor<1e-16)
344     {                                             260     {
345       G4ExceptionDescription ed;                  261       G4ExceptionDescription ed;
346       ed << "Inconsistent data of atomic compo << 262       ed << "Inconsistent data of atomic composition for " << 
347   material->GetName() << G4endl;                  263   material->GetName() << G4endl;
348       G4Exception("G4PenelopeRayleighModel::Bu    264       G4Exception("G4PenelopeRayleighModel::BuildFormFactorTable()",
349       "em2042",FatalException,ed);                265       "em2042",FatalException,ed);
350     }                                             266     }
351   //Normalize                                     267   //Normalize
352   for (std::size_t i=0;i<nElements;++i)        << 268   for (G4int i=0;i<nElements;i++)
353     (*StechiometricFactors)[i] /=  MaxStechiom    269     (*StechiometricFactors)[i] /=  MaxStechiometricFactor;
354                                                << 270  
                                                   >> 271   // Equivalent atoms per molecule
                                                   >> 272   G4double atomsPerMolecule = 0;
                                                   >> 273   for (G4int i=0;i<nElements;i++)
                                                   >> 274     atomsPerMolecule += (*StechiometricFactors)[i]; 
                                                   >> 275  
355   /*                                              276   /*
356     CREATE THE FORM FACTOR TABLE                  277     CREATE THE FORM FACTOR TABLE
357   */                                              278   */
358   G4PhysicsFreeVector* theFFVec = new G4Physic << 279   G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(logQSquareGrid.size());
                                                   >> 280   theFFVec->SetSpline(true);
359                                                   281 
360   for (std::size_t k=0;k<fLogQSquareGrid.size( << 282   for (size_t k=0;k<logQSquareGrid.size();k++)
361     {                                             283     {
362       G4double ff2 = 0; //squared form factor     284       G4double ff2 = 0; //squared form factor
363       for (std::size_t i=0;i<nElements;++i)    << 285       for (G4int i=0;i<nElements;i++)
364   {                                               286   {
365     G4int iZ = (*elementVector)[i]->GetZasInt( << 287     G4int iZ = (G4int) (*elementVector)[i]->GetZ();
366     G4PhysicsFreeVector* theAtomVec = fAtomicF << 288     G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second;
367     G4double f = (*theAtomVec)[k]; //the q-gri << 289     G4double f = (*theAtomVec)[k]; //the q-grid is always the same      
368     ff2 += f*f*(*StechiometricFactors)[i];        290     ff2 += f*f*(*StechiometricFactors)[i];
369   }                                               291   }
370       if (ff2)                                    292       if (ff2)
371   theFFVec->PutValue(k,fLogQSquareGrid[k],G4Lo << 293   theFFVec->PutValue(k,logQSquareGrid[k],std::log(ff2)); //NOTICE: THIS IS log(Q^2) vs. log(F^2)
372     }                                             294     }
373   theFFVec->FillSecondDerivatives(); //vector  << 295   logFormFactorTable->insert(std::make_pair(material,theFFVec));
374   fLogFormFactorTable->insert(std::make_pair(m << 
375                                                   296 
376   delete StechiometricFactors;                    297   delete StechiometricFactors;
                                                   >> 298   
377   return;                                         299   return;
378 }                                                 300 }
379                                                   301 
                                                   >> 302 
380 //....oooOO0OOooo........oooOO0OOooo........oo    303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381                                                   304 
382 void G4PenelopeRayleighModel::SampleSecondarie    305 void G4PenelopeRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
383             const G4MaterialCutsCouple* couple    306             const G4MaterialCutsCouple* couple,
384             const G4DynamicParticle* aDynamicG    307             const G4DynamicParticle* aDynamicGamma,
385             G4double,                             308             G4double,
386             G4double)                             309             G4double)
387 {                                                 310 {
388   // Sampling of the Rayleigh final state (nam << 311   // Sampling of the Rayleigh final state (namely, scattering angle of the photon) 
389   // from the Penelope2008 model. The scatteri << 312   // from the Penelope2008 model. The scattering angle is sampled from the atomic 
390   // cross section dOmega/d(cosTheta) from Bor << 313   // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding 
391   // anomalous scattering effects. The Form Fa << 314   // anomalous scattering effects. The Form Factor F(Q) function which appears in the 
392   // analytical cross section is retrieved via << 315   // analytical cross section is retrieved via the method GetFSquared(); atomic data 
393   // are tabulated for F(Q). Form factor for c << 316   // are tabulated for F(Q). Form factor for compounds is calculated according to 
394   // the additivity rule. The sampling from th << 317   // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse 
395   // Transform with Aliasing (RITA) algorithm; << 318   // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once 
396   // for each material and managed by G4Penelo    319   // for each material and managed by G4PenelopeSamplingData objects.
397   // The sampling algorithm (rejection method) << 320   // The sampling algorithm (rejection method) has efficiency 67% at low energy, and 
398   // increases with energy. For E=100 keV the  << 321   // increases with energy. For E=100 keV the efficiency is 100% and 86% for 
399   // hydrogen and uranium, respectively.          322   // hydrogen and uranium, respectively.
400                                                   323 
401   if (fVerboseLevel > 3)                       << 324   if (verboseLevel > 3)
402     G4cout << "Calling SamplingSecondaries() o    325     G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl;
403                                                   326 
404   G4double photonEnergy0 = aDynamicGamma->GetK    327   G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
405                                                << 328  
406   if (photonEnergy0 <= fIntrinsicLowEnergyLimi    329   if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
407     {                                             330     {
408       fParticleChange->ProposeTrackStatus(fSto    331       fParticleChange->ProposeTrackStatus(fStopAndKill);
409       fParticleChange->SetProposedKineticEnerg    332       fParticleChange->SetProposedKineticEnergy(0.);
410       fParticleChange->ProposeLocalEnergyDepos    333       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
411       return ;                                    334       return ;
412     }                                             335     }
413                                                   336 
414   G4ParticleMomentum photonDirection0 = aDynam    337   G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
415                                                << 338   
416   const G4Material* theMat = couple->GetMateri    339   const G4Material* theMat = couple->GetMaterial();
417                                                << 340   
418   //1) Verify if tables are ready                 341   //1) Verify if tables are ready
419   //Either Initialize() was not called, or we  << 342   if (!pMaxTable || !samplingTable)
420   //not invoked                                << 343     {
421   if (!fPMaxTable || !fSamplingTable || !fLogF << 344       G4Exception("G4PenelopeRayleighModel::SampleSecondaries()",
422     {                                          << 345       "em2043",FatalException,"Invalid model initialization");    
423       //create a **thread-local** version of t << 346       return;
424       //Unit Tests                             << 347     }
425       fLocalTable = true;                      << 348   
426       if (!fLogFormFactorTable)                << 349   //2) retrieve or build the sampling table
427   fLogFormFactorTable = new std::map<const G4M << 350   if (!(samplingTable->count(theMat)))
428       if (!fPMaxTable)                         << 351     InitializeSamplingAlgorithm(theMat);
429   fPMaxTable = new std::map<const G4Material*, << 352   G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second;
430       if (!fSamplingTable)                     << 353   
431   fSamplingTable = new std::map<const G4Materi << 354   //3) retrieve or build the pMax data
432     }                                          << 355   if (!pMaxTable->count(theMat))
433                                                << 356     GetPMaxTable(theMat);
434   if (!fSamplingTable->count(theMat))          << 357   G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second;
435     {                                          << 
436       //If we are here, it means that Initiali << 
437       //not filled up. This can happen in a Un << 
438       if (fVerboseLevel > 0)                   << 
439   {                                            << 
440     //Issue a G4Exception (warning) only in ve << 
441     G4ExceptionDescription ed;                 << 
442     ed << "Unable to find the fSamplingTable d << 
443       theMat->GetName() << G4endl;             << 
444     ed << "This can happen only in Unit Tests" << 
445     G4Exception("G4PenelopeRayleighModel::Samp << 
446           "em2019",JustWarning,ed);            << 
447   }                                            << 
448       const G4ElementVector* theElementVector  << 
449       //protect file reading via autolock      << 
450       G4AutoLock lock(&PenelopeRayleighModelMu << 
451       for (std::size_t j=0;j<theMat->GetNumber << 
452   {                                            << 
453     G4int iZ = theElementVector->at(j)->GetZas << 
454     if (!fLogAtomicCrossSection[iZ])           << 
455       {                                        << 
456         lock.lock();                           << 
457         ReadDataFile(iZ);                      << 
458         lock.unlock();                         << 
459       }                                        << 
460   }                                            << 
461       lock.lock();                             << 
462       //1) If the table has not been built for << 
463       if (!fLogFormFactorTable->count(theMat)) << 
464   BuildFormFactorTable(theMat);                << 
465                                                << 
466       //2) retrieve or build the sampling tabl << 
467       if (!(fSamplingTable->count(theMat)))    << 
468   InitializeSamplingAlgorithm(theMat);         << 
469                                                << 
470       //3) retrieve or build the pMax data     << 
471       if (!fPMaxTable->count(theMat))          << 
472   GetPMaxTable(theMat);                        << 
473       lock.unlock();                           << 
474     }                                          << 
475                                                << 
476   //Ok, restart the job                        << 
477   G4PenelopeSamplingData* theDataTable = fSamp << 
478   G4PhysicsFreeVector* thePMax = fPMaxTable->f << 
479                                                   358 
480   G4double cosTheta = 1.0;                        359   G4double cosTheta = 1.0;
481                                                << 360   
482   //OK, ready to go!                              361   //OK, ready to go!
483   G4double qmax = 2.0*photonEnergy0/electron_m    362   G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now
484                                                   363 
485   if (qmax < 1e-10) //very low momentum transf    364   if (qmax < 1e-10) //very low momentum transfer
486     {                                             365     {
487       G4bool loopAgain=false;                     366       G4bool loopAgain=false;
488       do                                          367       do
489   {                                               368   {
490     loopAgain = false;                            369     loopAgain = false;
491     cosTheta = 1.0-2.0*G4UniformRand();           370     cosTheta = 1.0-2.0*G4UniformRand();
492     G4double G = 0.5*(1+cosTheta*cosTheta);       371     G4double G = 0.5*(1+cosTheta*cosTheta);
493     if (G4UniformRand()>G)                        372     if (G4UniformRand()>G)
494       loopAgain = true;                           373       loopAgain = true;
495   }while(loopAgain);                              374   }while(loopAgain);
496     }                                             375     }
497   else //larger momentum transfer                 376   else //larger momentum transfer
498     {                                             377     {
499       std::size_t nData = theDataTable->GetNum << 378       size_t nData = theDataTable->GetNumberOfStoredPoints();
500       G4double LastQ2inTheTable = theDataTable    379       G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
501       G4double q2max = std::min(qmax*qmax,Last    380       G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
502                                                   381 
503       G4bool loopAgain = false;                   382       G4bool loopAgain = false;
504       G4double MaxPValue = thePMax->Value(phot    383       G4double MaxPValue = thePMax->Value(photonEnergy0);
505       G4double xx=0;                              384       G4double xx=0;
506                                                << 385       
507       //Sampling by rejection method. The reje << 386       //Sampling by rejection method. The rejection function is 
508       //G = 0.5*(1+cos^2(theta))                  387       //G = 0.5*(1+cos^2(theta))
509       //                                          388       //
510       do{                                         389       do{
511   loopAgain = false;                              390   loopAgain = false;
512   G4double RandomMax = G4UniformRand()*MaxPVal    391   G4double RandomMax = G4UniformRand()*MaxPValue;
513   xx = theDataTable->SampleValue(RandomMax);      392   xx = theDataTable->SampleValue(RandomMax);
514   //xx is a random value of q^2 in (0,q2max),s << 393   //xx is a random value of q^2 in (0,q2max),sampled according to 
515   //F(Q^2) via the RITA algorithm                 394   //F(Q^2) via the RITA algorithm
516   if (xx > q2max)                                 395   if (xx > q2max)
517     loopAgain = true;                             396     loopAgain = true;
518   cosTheta = 1.0-2.0*xx/q2max;                    397   cosTheta = 1.0-2.0*xx/q2max;
519   G4double G = 0.5*(1+cosTheta*cosTheta);         398   G4double G = 0.5*(1+cosTheta*cosTheta);
520   if (G4UniformRand()>G)                          399   if (G4UniformRand()>G)
521     loopAgain = true;                             400     loopAgain = true;
522       }while(loopAgain);                          401       }while(loopAgain);
523     }                                             402     }
524                                                << 403   
525   G4double sinTheta = std::sqrt(1-cosTheta*cos    404   G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
526                                                << 405  
527   // Scattered photon angles. ( Z - axis along    406   // Scattered photon angles. ( Z - axis along the parent photon)
528   G4double phi = twopi * G4UniformRand() ;        407   G4double phi = twopi * G4UniformRand() ;
529   G4double dirX = sinTheta*std::cos(phi);         408   G4double dirX = sinTheta*std::cos(phi);
530   G4double dirY = sinTheta*std::sin(phi);         409   G4double dirY = sinTheta*std::sin(phi);
531   G4double dirZ = cosTheta;                       410   G4double dirZ = cosTheta;
532                                                << 411   
533   // Update G4VParticleChange for the scattere << 412   // Update G4VParticleChange for the scattered photon 
534   G4ThreeVector photonDirection1(dirX, dirY, d    413   G4ThreeVector photonDirection1(dirX, dirY, dirZ);
535                                                   414 
536   photonDirection1.rotateUz(photonDirection0);    415   photonDirection1.rotateUz(photonDirection0);
537   fParticleChange->ProposeMomentumDirection(ph    416   fParticleChange->ProposeMomentumDirection(photonDirection1) ;
538   fParticleChange->SetProposedKineticEnergy(ph    417   fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
539                                                << 418  
540   return;                                         419   return;
541 }                                                 420 }
542                                                   421 
543                                                   422 
544 //....oooOO0OOooo........oooOO0OOooo........oo    423 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
545                                                   424 
546 void G4PenelopeRayleighModel::ReadDataFile(con    425 void G4PenelopeRayleighModel::ReadDataFile(const G4int Z)
547 {                                                 426 {
548   if (fVerboseLevel > 2)                       << 427   if (verboseLevel > 2)
549     {                                             428     {
550       G4cout << "G4PenelopeRayleighModel::Read    429       G4cout << "G4PenelopeRayleighModel::ReadDataFile()" << G4endl;
551       G4cout << "Going to read Rayleigh data f    430       G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl;
552     }                                             431     }
553     const char* path = G4FindDataDir("G4LEDATA << 432 
554     if(!path)                                  << 433   char* path = getenv("G4LEDATA");
                                                   >> 434   if (!path)
555     {                                             435     {
556       G4String excep = "G4LEDATA environment v    436       G4String excep = "G4LEDATA environment variable not set!";
557       G4Exception("G4PenelopeRayleighModel::Re    437       G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
558       "em0006",FatalException,excep);             438       "em0006",FatalException,excep);
559       return;                                     439       return;
560     }                                             440     }
561                                                   441 
562   /*                                              442   /*
563     Read first the cross section file             443     Read first the cross section file
564   */                                              444   */
565   std::ostringstream ost;                         445   std::ostringstream ost;
566   if (Z>9)                                        446   if (Z>9)
567     ost << path << "/penelope/rayleigh/pdgra"     447     ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08";
568   else                                            448   else
569     ost << path << "/penelope/rayleigh/pdgra0"    449     ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08";
570   std::ifstream file(ost.str().c_str());          450   std::ifstream file(ost.str().c_str());
571   if (!file.is_open())                            451   if (!file.is_open())
572     {                                             452     {
573       G4String excep = "Data file " + G4String    453       G4String excep = "Data file " + G4String(ost.str()) + " not found!";
574       G4Exception("G4PenelopeRayleighModel::Re    454       G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
575       "em0003",FatalException,excep);             455       "em0003",FatalException,excep);
576     }                                             456     }
577   G4int readZ =0;                                 457   G4int readZ =0;
578   std::size_t nPoints= 0;                      << 458   size_t nPoints= 0;
579   file >> readZ >> nPoints;                       459   file >> readZ >> nPoints;
580   //check the right file is opened.               460   //check the right file is opened.
581   if (readZ != Z || nPoints <= 0 || nPoints >=    461   if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
582     {                                             462     {
583       G4ExceptionDescription ed;                  463       G4ExceptionDescription ed;
584       ed << "Corrupted data file for Z=" << Z     464       ed << "Corrupted data file for Z=" << Z << G4endl;
585       G4Exception("G4PenelopeRayleighModel::Re    465       G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
586       "em0005",FatalException,ed);                466       "em0005",FatalException,ed);
587       return;                                     467       return;
588     }                                          << 468     }  
589                                                << 469   G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints);
590   fLogAtomicCrossSection[Z] = new G4PhysicsFre << 
591   G4double ene=0,f1=0,f2=0,xs=0;                  470   G4double ene=0,f1=0,f2=0,xs=0;
592   for (std::size_t i=0;i<nPoints;++i)          << 471   for (size_t i=0;i<nPoints;i++)
593     {                                             472     {
594       file >> ene >> f1 >> f2 >> xs;              473       file >> ene >> f1 >> f2 >> xs;
595       //dimensional quantities                    474       //dimensional quantities
596       ene *= eV;                                  475       ene *= eV;
597       xs *= cm2;                                  476       xs *= cm2;
598       fLogAtomicCrossSection[Z]->PutValue(i,G4 << 477       theVec->PutValue(i,std::log(ene),std::log(xs));
599       if (file.eof() && i != (nPoints-1)) //fi    478       if (file.eof() && i != (nPoints-1)) //file ended too early
600   {                                               479   {
601     G4ExceptionDescription ed ;                << 480     G4ExceptionDescription ed ;   
602     ed << "Corrupted data file for Z=" << Z <<    481     ed << "Corrupted data file for Z=" << Z << G4endl;
603     ed << "Found less than " << nPoints << "en    482     ed << "Found less than " << nPoints << "entries " <<G4endl;
604     G4Exception("G4PenelopeRayleighModel::Read    483     G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
605           "em0005",FatalException,ed);            484           "em0005",FatalException,ed);
606   }                                               485   }
607     }                                             486     }
                                                   >> 487   if (!logAtomicCrossSection)
                                                   >> 488     {
                                                   >> 489       G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
                                                   >> 490       "em2044",FatalException,"Unable to allocate the atomic cross section table");
                                                   >> 491       delete theVec;
                                                   >> 492       return;
                                                   >> 493     }
                                                   >> 494   logAtomicCrossSection->insert(std::make_pair(Z,theVec));
608   file.close();                                   495   file.close();
609                                                   496 
610   /*                                              497   /*
611     Then read the form factor file                498     Then read the form factor file
612   */                                              499   */
613   std::ostringstream ost2;                        500   std::ostringstream ost2;
614   if (Z>9)                                        501   if (Z>9)
615     ost2 << path << "/penelope/rayleigh/pdaff"    502     ost2 << path << "/penelope/rayleigh/pdaff" << Z << ".p08";
616   else                                            503   else
617     ost2 << path << "/penelope/rayleigh/pdaff0    504     ost2 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08";
618   file.open(ost2.str().c_str());                  505   file.open(ost2.str().c_str());
619   if (!file.is_open())                            506   if (!file.is_open())
620     {                                             507     {
621       G4String excep = "Data file " + G4String    508       G4String excep = "Data file " + G4String(ost2.str()) + " not found!";
622       G4Exception("G4PenelopeRayleighModel::Re    509       G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
623       "em0003",FatalException,excep);             510       "em0003",FatalException,excep);
624     }                                             511     }
625   file >> readZ >> nPoints;                       512   file >> readZ >> nPoints;
626   //check the right file is opened.               513   //check the right file is opened.
627   if (readZ != Z || nPoints <= 0 || nPoints >=    514   if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
628     {                                             515     {
629       G4ExceptionDescription ed;                  516       G4ExceptionDescription ed;
630       ed << "Corrupted data file for Z=" << Z     517       ed << "Corrupted data file for Z=" << Z << G4endl;
631       G4Exception("G4PenelopeRayleighModel::Re    518       G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
632       "em0005",FatalException,ed);                519       "em0005",FatalException,ed);
633       return;                                     520       return;
634     }                                          << 521     }  
635   fAtomicFormFactor[Z] = new G4PhysicsFreeVect << 522   G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints);
636   G4double q=0,ff=0,incoh=0;                      523   G4double q=0,ff=0,incoh=0;
637   G4bool fillQGrid = false;                       524   G4bool fillQGrid = false;
638   //fill this vector only the first time.         525   //fill this vector only the first time.
639   if (!fLogQSquareGrid.size())                 << 526   if (!logQSquareGrid.size())
640     fillQGrid = true;                             527     fillQGrid = true;
641   for (std::size_t i=0;i<nPoints;++i)          << 528   for (size_t i=0;i<nPoints;i++)
642     {                                             529     {
643       file >> q >> ff >> incoh;                   530       file >> q >> ff >> incoh;
644       //q and ff are dimensionless (q is in un    531       //q and ff are dimensionless (q is in units of (m_e*c)
645       fAtomicFormFactor[Z]->PutValue(i,q,ff);  << 532       theFFVec->PutValue(i,q,ff);
646       if (fillQGrid)                              533       if (fillQGrid)
647   {                                               534   {
648     fLogQSquareGrid.push_back(2.0*G4Log(q));   << 535     logQSquareGrid.push_back(2.0*std::log(q));
649   }                                               536   }
650       if (file.eof() && i != (nPoints-1)) //fi    537       if (file.eof() && i != (nPoints-1)) //file ended too early
651   {                                               538   {
652     G4ExceptionDescription ed;                 << 539     G4ExceptionDescription ed;    
653     ed << "Corrupted data file for Z=" << Z <<    540     ed << "Corrupted data file for Z=" << Z << G4endl;
654     ed << "Found less than " << nPoints << "en    541     ed << "Found less than " << nPoints << "entries " <<G4endl;
655     G4Exception("G4PenelopeRayleighModel::Read    542     G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
656           "em0005",FatalException,ed);            543           "em0005",FatalException,ed);
657   }                                               544   }
658     }                                             545     }
                                                   >> 546   if (!atomicFormFactor)
                                                   >> 547     {
                                                   >> 548       G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
                                                   >> 549       "em2045",FatalException,
                                                   >> 550       "Unable to allocate the atomicFormFactor data table");
                                                   >> 551       delete theFFVec;
                                                   >> 552       return;
                                                   >> 553     }
                                                   >> 554   atomicFormFactor->insert(std::make_pair(Z,theFFVec));
659   file.close();                                   555   file.close();
660   return;                                         556   return;
661 }                                                 557 }
662                                                   558 
663 //....oooOO0OOooo........oooOO0OOooo........oo    559 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
664                                                   560 
665 G4double G4PenelopeRayleighModel::GetFSquared(    561 G4double G4PenelopeRayleighModel::GetFSquared(const G4Material* mat, const G4double QSquared)
666 {                                                 562 {
667   G4double f2 = 0;                                563   G4double f2 = 0;
668   //Input value QSquared could be zero: protec << 564   //Input value QSquared could be zero: protect the log() below against 
669   //the FPE exception                             565   //the FPE exception
670   //If Q<1e-10, set Q to 1e-10                    566   //If Q<1e-10, set Q to 1e-10
671   G4double logQSquared = (QSquared>1e-10) ? G4 << 567   G4double logQSquared = (QSquared>1e-10) ? std::log(QSquared) : -23.;
672   //last value of the table                       568   //last value of the table
673   G4double maxlogQ2 = fLogQSquareGrid[fLogQSqu << 569   G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1];
674                                                << 570   //If the table has not been built for the material, do it!
                                                   >> 571   if (!logFormFactorTable->count(mat))
                                                   >> 572     BuildFormFactorTable(mat);
                                                   >> 573   
675   //now it should  be all right                   574   //now it should  be all right
676   G4PhysicsFreeVector* theVec = fLogFormFactor << 575   G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
677                                                   576 
678   if (!theVec)                                    577   if (!theVec)
679     {                                             578     {
680       G4ExceptionDescription ed;                  579       G4ExceptionDescription ed;
681       ed << "Unable to retrieve F squared tabl    580       ed << "Unable to retrieve F squared table for " << mat->GetName() << G4endl;
682       G4Exception("G4PenelopeRayleighModel::Ge    581       G4Exception("G4PenelopeRayleighModel::GetFSquared()",
683       "em2046",FatalException,ed);                582       "em2046",FatalException,ed);
684       return 0;                                   583       return 0;
685     }                                             584     }
686   if (logQSquared < -20) // Q < 1e-9              585   if (logQSquared < -20) // Q < 1e-9
687     {                                             586     {
688       G4double logf2 = (*theVec)[0]; //first v    587       G4double logf2 = (*theVec)[0]; //first value of the table
689       f2 = G4Exp(logf2);                       << 588       f2 = std::exp(logf2);
690     }                                             589     }
691   else if (logQSquared > maxlogQ2)                590   else if (logQSquared > maxlogQ2)
692     f2 =0;                                        591     f2 =0;
693   else                                            592   else
694     {                                             593     {
695       //log(Q^2) vs. log(F^2)                     594       //log(Q^2) vs. log(F^2)
696       G4double logf2 = theVec->Value(logQSquar    595       G4double logf2 = theVec->Value(logQSquared);
697       f2 = G4Exp(logf2);                       << 596       f2 = std::exp(logf2);
698                                                   597 
699     }                                             598     }
700   if (fVerboseLevel > 3)                       << 599   if (verboseLevel > 3)
701     {                                             600     {
702       G4cout << "G4PenelopeRayleighModel::GetF << 601       G4cout << "G4PenelopeRayleighModel::GetFSquared() in " << mat->GetName() << G4endl;  
703       G4cout << "Q^2 = " <<  QSquared << " (un    602       G4cout << "Q^2 = " <<  QSquared << " (units of 1/(m_e*c); F^2 = " << f2 << G4endl;
704     }                                             603     }
705   return f2;                                      604   return f2;
706 }                                                 605 }
707                                                   606 
708 //....oooOO0OOooo........oooOO0OOooo........oo    607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
709                                                   608 
710 void G4PenelopeRayleighModel::InitializeSampli    609 void G4PenelopeRayleighModel::InitializeSamplingAlgorithm(const G4Material* mat)
711 {                                                 610 {
712   G4double q2min = 0;                             611   G4double q2min = 0;
713   G4double q2max = 0;                             612   G4double q2max = 0;
714   const std::size_t np = 150; //hard-coded in  << 613   const size_t np = 150; //hard-coded in Penelope
715   //G4cout << "Init N= " << fLogQSquareGrid.si << 614   for (size_t i=1;i<logQSquareGrid.size();i++)
716   for (std::size_t i=1;i<fLogQSquareGrid.size( << 
717     {                                             615     {
718       G4double Q2 = G4Exp(fLogQSquareGrid[i]); << 616       G4double Q2 = std::exp(logQSquareGrid[i]);
719       if (GetFSquared(mat,Q2) >  1e-35)           617       if (GetFSquared(mat,Q2) >  1e-35)
720   {                                               618   {
721     q2max = G4Exp(fLogQSquareGrid[i-1]);       << 619     q2max = std::exp(logQSquareGrid[i-1]);
722   }                                               620   }
723       //G4cout << "Q2= " << Q2 << " q2max= " < << 
724     }                                             621     }
725                                                << 622   
726   std::size_t nReducedPoints = np/4;           << 623   size_t nReducedPoints = np/4;
727                                                   624 
728   //check for errors                              625   //check for errors
729   if (np < 16)                                    626   if (np < 16)
730     {                                             627     {
731       G4Exception("G4PenelopeRayleighModel::In    628       G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
732       "em2047",FatalException,                    629       "em2047",FatalException,
733       "Too few points to initialize the sampli    630       "Too few points to initialize the sampling algorithm");
734     }                                             631     }
735   if (q2min > (q2max-1e-10))                      632   if (q2min > (q2max-1e-10))
736     {                                             633     {
737       G4cout << "q2min= " << q2min << " q2max= << 
738       G4Exception("G4PenelopeRayleighModel::In    634       G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
739       "em2048",FatalException,                    635       "em2048",FatalException,
740       "Too narrow grid to initialize the sampl    636       "Too narrow grid to initialize the sampling algorithm");
741     }                                             637     }
742                                                   638 
743   //This is subroutine RITAI0 of Penelope         639   //This is subroutine RITAI0 of Penelope
744   //Create an object of type G4PenelopeRayleig    640   //Create an object of type G4PenelopeRayleighSamplingData --> store in a map::Material*
745                                                   641 
746   //temporary vectors --> Then everything is s    642   //temporary vectors --> Then everything is stored in G4PenelopeSamplingData
747   G4DataVector* x = new G4DataVector();           643   G4DataVector* x = new G4DataVector();
748                                                << 644   
749   /*******************************************    645   /*******************************************************************************
750     Start with a grid of NUNIF points uniforml    646     Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max
751   ********************************************    647   ********************************************************************************/
752   std::size_t NUNIF = std::min(std::max(((std: << 648   size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2);
753   const G4int nip = 51; //hard-coded in Penelo    649   const G4int nip = 51; //hard-coded in Penelope
754                                                   650 
755   G4double dx = (q2max-q2min)/((G4double) NUNI << 651   G4double dx = (q2max-q2min)/((G4double) NUNIF-1);  
756   x->push_back(q2min);                         << 652   x->push_back(q2min); 
757   for (std::size_t i=1;i<NUNIF-1;++i)          << 653   for (size_t i=1;i<NUNIF-1;i++)
758     {                                             654     {
759       G4double app = q2min + i*dx;                655       G4double app = q2min + i*dx;
760       x->push_back(app); //increase            << 656       x->push_back(app); //increase 
761     }                                             657     }
762   x->push_back(q2max);                            658   x->push_back(q2max);
763                                                << 659   
764   if (fVerboseLevel> 3)                        << 660   if (verboseLevel> 3)
765     G4cout << "Vector x has " << x->size() <<     661     G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl;
766                                                   662 
767   //Allocate temporary storage vectors            663   //Allocate temporary storage vectors
768   G4DataVector* area = new G4DataVector();        664   G4DataVector* area = new G4DataVector();
769   G4DataVector* a = new G4DataVector();           665   G4DataVector* a = new G4DataVector();
770   G4DataVector* b = new G4DataVector();           666   G4DataVector* b = new G4DataVector();
771   G4DataVector* c = new G4DataVector();           667   G4DataVector* c = new G4DataVector();
772   G4DataVector* err = new G4DataVector();         668   G4DataVector* err = new G4DataVector();
773                                                   669 
774   for (std::size_t i=0;i<NUNIF-1;++i) //build  << 670   for (size_t i=0;i<NUNIF-1;i++) //build all points but the last
775     {                                          << 671     {      
776       //Temporary vectors for this loop           672       //Temporary vectors for this loop
777       G4DataVector* pdfi = new G4DataVector();    673       G4DataVector* pdfi = new G4DataVector();
778       G4DataVector* pdfih = new G4DataVector()    674       G4DataVector* pdfih = new G4DataVector();
779       G4DataVector* sumi = new G4DataVector();    675       G4DataVector* sumi = new G4DataVector();
780                                                   676 
781       G4double dxi = ((*x)[i+1]-(*x)[i])/(G4do    677       G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
782       G4double pdfmax = 0;                        678       G4double pdfmax = 0;
783       for (G4int k=0;k<nip;k++)                   679       for (G4int k=0;k<nip;k++)
784   {                                               680   {
785     G4double xik = (*x)[i]+k*dxi;              << 681     G4double xik = (*x)[i]+k*dxi; 
786     G4double pdfk = std::max(GetFSquared(mat,x    682     G4double pdfk = std::max(GetFSquared(mat,xik),0.);
787     pdfi->push_back(pdfk);                        683     pdfi->push_back(pdfk);
788     pdfmax = std::max(pdfmax,pdfk);            << 684     pdfmax = std::max(pdfmax,pdfk); 
789     if (k < (nip-1))                              685     if (k < (nip-1))
790       {                                           686       {
791         G4double xih = xik + 0.5*dxi;             687         G4double xih = xik + 0.5*dxi;
792         G4double pdfIK = std::max(GetFSquared(    688         G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
793         pdfih->push_back(pdfIK);                  689         pdfih->push_back(pdfIK);
794         pdfmax = std::max(pdfmax,pdfIK);          690         pdfmax = std::max(pdfmax,pdfIK);
795       }                                           691       }
796   }                                               692   }
797                                                << 693     
798       //Simpson's integration                     694       //Simpson's integration
799       G4double cons = dxi*0.5*(1./3.);            695       G4double cons = dxi*0.5*(1./3.);
800       sumi->push_back(0.);                        696       sumi->push_back(0.);
801       for (G4int k=1;k<nip;k++)                   697       for (G4int k=1;k<nip;k++)
802   {                                               698   {
803     G4double previous = (*sumi)[k-1];             699     G4double previous = (*sumi)[k-1];
804     G4double next = previous + cons*((*pdfi)[k    700     G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
805     sumi->push_back(next);                        701     sumi->push_back(next);
806   }                                               702   }
807                                                << 703     
808       G4double lastIntegral = (*sumi)[sumi->si    704       G4double lastIntegral = (*sumi)[sumi->size()-1];
809       area->push_back(lastIntegral);              705       area->push_back(lastIntegral);
810       //Normalize cumulative function             706       //Normalize cumulative function
811       G4double factor = 1.0/lastIntegral;         707       G4double factor = 1.0/lastIntegral;
812       for (std::size_t k=0;k<sumi->size();++k) << 708       for (size_t k=0;k<sumi->size();k++)
813   (*sumi)[k] *= factor;                           709   (*sumi)[k] *= factor;
814                                                << 710       
815       //When the PDF vanishes at one of the in    711       //When the PDF vanishes at one of the interval end points, its value is modified
816       if ((*pdfi)[0] < 1e-35)                  << 712       if ((*pdfi)[0] < 1e-35) 
817   (*pdfi)[0] = 1e-5*pdfmax;                       713   (*pdfi)[0] = 1e-5*pdfmax;
818       if ((*pdfi)[pdfi->size()-1] < 1e-35)        714       if ((*pdfi)[pdfi->size()-1] < 1e-35)
819   (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;          715   (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
820                                                   716 
821       G4double pli = (*pdfi)[0]*factor;           717       G4double pli = (*pdfi)[0]*factor;
822       G4double pui = (*pdfi)[pdfi->size()-1]*f    718       G4double pui = (*pdfi)[pdfi->size()-1]*factor;
823       G4double B_temp = 1.0-1.0/(pli*pui*dx*dx    719       G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
824       G4double A_temp = (1.0/(pli*dx))-1.0-B_t    720       G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
825       G4double C_temp = 1.0+A_temp+B_temp;        721       G4double C_temp = 1.0+A_temp+B_temp;
826       if (C_temp < 1e-35)                         722       if (C_temp < 1e-35)
827   {                                               723   {
828     a->push_back(0.);                             724     a->push_back(0.);
829     b->push_back(0.);                             725     b->push_back(0.);
830     c->push_back(1.);                          << 726     c->push_back(1.);   
831   }                                               727   }
832       else                                        728       else
833   {                                               729   {
834     a->push_back(A_temp);                         730     a->push_back(A_temp);
835     b->push_back(B_temp);                         731     b->push_back(B_temp);
836     c->push_back(C_temp);                         732     c->push_back(C_temp);
837   }                                               733   }
838                                                   734 
839       //OK, now get ERR(I), the integral of th << 735       //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation 
840       //and the true pdf, extended over the in    736       //and the true pdf, extended over the interval (X(I),X(I+1))
841       G4int icase = 1; //loop code                737       G4int icase = 1; //loop code
842       G4bool reLoop = false;                      738       G4bool reLoop = false;
843       err->push_back(0.);                         739       err->push_back(0.);
844       do                                          740       do
845   {                                               741   {
846     reLoop = false;                               742     reLoop = false;
847     (*err)[i] = 0.; //zero variable               743     (*err)[i] = 0.; //zero variable
848     for (G4int k=0;k<nip;k++)                     744     for (G4int k=0;k<nip;k++)
849       {                                           745       {
850         G4double rr = (*sumi)[k];                 746         G4double rr = (*sumi)[k];
851         G4double pap = (*area)[i]*(1.0+((*a)[i    747         G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
852     ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*    748     ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
853         if (k == 0 || k == nip-1)                 749         if (k == 0 || k == nip-1)
854     (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k])    750     (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
855         else                                      751         else
856     (*err)[i] += std::fabs(pap-(*pdfi)[k]);       752     (*err)[i] += std::fabs(pap-(*pdfi)[k]);
857       }                                           753       }
858     (*err)[i] *= dxi;                             754     (*err)[i] *= dxi;
859                                                << 755       
860     //If err(I) is too large, the pdf is appro    756     //If err(I) is too large, the pdf is approximated by a uniform distribution
861     if ((*err)[i] > 0.1*(*area)[i] && icase == << 757     if ((*err)[i] > 0.1*(*area)[i] && icase == 1) 
862       {                                           758       {
863         (*b)[i] = 0;                              759         (*b)[i] = 0;
864         (*a)[i] = 0;                              760         (*a)[i] = 0;
865         (*c)[i] = 1.;                             761         (*c)[i] = 1.;
866         icase = 2;                                762         icase = 2;
867         reLoop = true;                            763         reLoop = true;
868       }                                           764       }
869   }while(reLoop);                                 765   }while(reLoop);
                                                   >> 766 
870       delete pdfi;                                767       delete pdfi;
871       delete pdfih;                               768       delete pdfih;
872       delete sumi;                                769       delete sumi;
873     } //end of first loop over i                  770     } //end of first loop over i
874                                                   771 
875   //Now assign last point                         772   //Now assign last point
876   (*x)[x->size()-1] = q2max;                      773   (*x)[x->size()-1] = q2max;
877   a->push_back(0.);                               774   a->push_back(0.);
878   b->push_back(0.);                               775   b->push_back(0.);
879   c->push_back(0.);                               776   c->push_back(0.);
880   err->push_back(0.);                             777   err->push_back(0.);
881   area->push_back(0.);                            778   area->push_back(0.);
882                                                   779 
883   if (x->size() != NUNIF || a->size() != NUNIF << 780   if (x->size() != NUNIF || a->size() != NUNIF || 
884       err->size() != NUNIF || area->size() !=     781       err->size() != NUNIF || area->size() != NUNIF)
885     {                                             782     {
886       G4ExceptionDescription ed;                  783       G4ExceptionDescription ed;
887       ed << "Problem in building the Table for    784       ed << "Problem in building the Table for Sampling: array dimensions do not match" << G4endl;
888       G4Exception("G4PenelopeRayleighModel::In    785       G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
889       "em2049",FatalException,ed);                786       "em2049",FatalException,ed);
890     }                                             787     }
891                                                << 788   
892   /*******************************************    789   /*******************************************************************************
893    New grid points are added by halving the su    790    New grid points are added by halving the sub-intervals with the largest absolute error
894   This is done up to np=150 points in the grid    791   This is done up to np=150 points in the grid
895   ********************************************    792   ********************************************************************************/
896   do                                              793   do
897     {                                             794     {
898       G4double maxError = 0.0;                    795       G4double maxError = 0.0;
899       std::size_t iErrMax = 0;                 << 796       size_t iErrMax = 0;
900       for (std::size_t i=0;i<err->size()-2;++i << 797       for (size_t i=0;i<err->size()-2;i++) 
901   {                                               798   {
902     //maxError is the lagest of the interval e    799     //maxError is the lagest of the interval errors err[i]
903     if ((*err)[i] > maxError)                     800     if ((*err)[i] > maxError)
904       {                                           801       {
905         maxError = (*err)[i];                     802         maxError = (*err)[i];
906         iErrMax = i;                              803         iErrMax = i;
907       }                                           804       }
908   }                                               805   }
909                                                << 806       
910       //OK, now I have to insert one new point    807       //OK, now I have to insert one new point in the position iErrMax
911       G4double newx = 0.5*((*x)[iErrMax]+(*x)[    808       G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
912                                                << 809       
913       x->insert(x->begin()+iErrMax+1,newx);       810       x->insert(x->begin()+iErrMax+1,newx);
914       //Add place-holders in the other vectors    811       //Add place-holders in the other vectors
915       area->insert(area->begin()+iErrMax+1,0.)    812       area->insert(area->begin()+iErrMax+1,0.);
916       a->insert(a->begin()+iErrMax+1,0.);         813       a->insert(a->begin()+iErrMax+1,0.);
917       b->insert(b->begin()+iErrMax+1,0.);         814       b->insert(b->begin()+iErrMax+1,0.);
918       c->insert(c->begin()+iErrMax+1,0.);         815       c->insert(c->begin()+iErrMax+1,0.);
919       err->insert(err->begin()+iErrMax+1,0.);     816       err->insert(err->begin()+iErrMax+1,0.);
920                                                << 817         
921       //Now calculate the other parameters        818       //Now calculate the other parameters
922       for (std::size_t i=iErrMax;i<=iErrMax+1; << 819       for (size_t i=iErrMax;i<=iErrMax+1;i++)
923   {                                               820   {
924     //Temporary vectors for this loop             821     //Temporary vectors for this loop
925     G4DataVector* pdfi = new G4DataVector();      822     G4DataVector* pdfi = new G4DataVector();
926     G4DataVector* pdfih = new G4DataVector();     823     G4DataVector* pdfih = new G4DataVector();
927     G4DataVector* sumi = new G4DataVector();      824     G4DataVector* sumi = new G4DataVector();
928                                                << 825     
929     G4double dxLocal = (*x)[i+1]-(*x)[i];      << 826     G4double dx = (*x)[i+1]-(*x)[i];
930     G4double dxi = ((*x)[i+1]-(*x)[i])/(G4doub    827     G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
931     G4double pdfmax = 0;                          828     G4double pdfmax = 0;
932     for (G4int k=0;k<nip;k++)                     829     for (G4int k=0;k<nip;k++)
933       {                                           830       {
934         G4double xik = (*x)[i]+k*dxi;             831         G4double xik = (*x)[i]+k*dxi;
935         G4double pdfk = std::max(GetFSquared(m    832         G4double pdfk = std::max(GetFSquared(mat,xik),0.);
936         pdfi->push_back(pdfk);                    833         pdfi->push_back(pdfk);
937         pdfmax = std::max(pdfmax,pdfk);        << 834         pdfmax = std::max(pdfmax,pdfk); 
938         if (k < (nip-1))                          835         if (k < (nip-1))
939     {                                             836     {
940       G4double xih = xik + 0.5*dxi;               837       G4double xih = xik + 0.5*dxi;
941       G4double pdfIK = std::max(GetFSquared(ma    838       G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
942       pdfih->push_back(pdfIK);                    839       pdfih->push_back(pdfIK);
943       pdfmax = std::max(pdfmax,pdfIK);            840       pdfmax = std::max(pdfmax,pdfIK);
944     }                                             841     }
945       }                                           842       }
946                                                << 843     
947     //Simpson's integration                       844     //Simpson's integration
948     G4double cons = dxi*0.5*(1./3.);              845     G4double cons = dxi*0.5*(1./3.);
949     sumi->push_back(0.);                          846     sumi->push_back(0.);
950     for (G4int k=1;k<nip;k++)                     847     for (G4int k=1;k<nip;k++)
951       {                                           848       {
952         G4double previous = (*sumi)[k-1];         849         G4double previous = (*sumi)[k-1];
953         G4double next = previous + cons*((*pdf    850         G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
954         sumi->push_back(next);                    851         sumi->push_back(next);
955       }                                           852       }
956     G4double lastIntegral = (*sumi)[sumi->size    853     G4double lastIntegral = (*sumi)[sumi->size()-1];
957     (*area)[i] = lastIntegral;                    854     (*area)[i] = lastIntegral;
958                                                << 855     
959     //Normalize cumulative function               856     //Normalize cumulative function
960     G4double factor = 1.0/lastIntegral;           857     G4double factor = 1.0/lastIntegral;
961     for (std::size_t k=0;k<sumi->size();++k)   << 858     for (size_t k=0;k<sumi->size();k++)
962       (*sumi)[k] *= factor;                       859       (*sumi)[k] *= factor;
963                                                << 860     
964     //When the PDF vanishes at one of the inte    861     //When the PDF vanishes at one of the interval end points, its value is modified
965     if ((*pdfi)[0] < 1e-35)                    << 862     if ((*pdfi)[0] < 1e-35) 
966       (*pdfi)[0] = 1e-5*pdfmax;                   863       (*pdfi)[0] = 1e-5*pdfmax;
967     if ((*pdfi)[pdfi->size()-1] < 1e-35)          864     if ((*pdfi)[pdfi->size()-1] < 1e-35)
968       (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;      865       (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
969                                                << 866     
970     G4double pli = (*pdfi)[0]*factor;             867     G4double pli = (*pdfi)[0]*factor;
971     G4double pui = (*pdfi)[pdfi->size()-1]*fac    868     G4double pui = (*pdfi)[pdfi->size()-1]*factor;
972     G4double B_temp = 1.0-1.0/(pli*pui*dxLocal << 869     G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
973     G4double A_temp = (1.0/(pli*dxLocal))-1.0- << 870     G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
974     G4double C_temp = 1.0+A_temp+B_temp;          871     G4double C_temp = 1.0+A_temp+B_temp;
975     if (C_temp < 1e-35)                           872     if (C_temp < 1e-35)
976       {                                           873       {
977         (*a)[i]= 0.;                              874         (*a)[i]= 0.;
978         (*b)[i] = 0.;                             875         (*b)[i] = 0.;
979         (*c)[i] = 1;                              876         (*c)[i] = 1;
980       }                                           877       }
981     else                                          878     else
982       {                                           879       {
983         (*a)[i]= A_temp;                          880         (*a)[i]= A_temp;
984         (*b)[i] = B_temp;                         881         (*b)[i] = B_temp;
985         (*c)[i] = C_temp;                         882         (*c)[i] = C_temp;
986       }                                           883       }
987     //OK, now get ERR(I), the integral of the  << 884     //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation 
988     //and the true pdf, extended over the inte    885     //and the true pdf, extended over the interval (X(I),X(I+1))
989     G4int icase = 1; //loop code                  886     G4int icase = 1; //loop code
990     G4bool reLoop = false;                        887     G4bool reLoop = false;
991     do                                            888     do
992       {                                           889       {
993         reLoop = false;                           890         reLoop = false;
994         (*err)[i] = 0.; //zero variable           891         (*err)[i] = 0.; //zero variable
995         for (G4int k=0;k<nip;k++)                 892         for (G4int k=0;k<nip;k++)
996     {                                             893     {
997       G4double rr = (*sumi)[k];                << 894       G4double rr = (*sumi)[k];       
998       G4double pap = (*area)[i]*(1.0+((*a)[i]+    895       G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
999         ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1    896         ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1000       if (k == 0 || k == nip-1)                  897       if (k == 0 || k == nip-1)
1001         (*err)[i] += 0.5*std::fabs(pap-(*pdfi    898         (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1002       else                                       899       else
1003         (*err)[i] += std::fabs(pap-(*pdfi)[k]    900         (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1004     }                                            901     }
1005         (*err)[i] *= dxi;                        902         (*err)[i] *= dxi;
1006                                               << 903         
1007         //If err(I) is too large, the pdf is     904         //If err(I) is too large, the pdf is approximated by a uniform distribution
1008         if ((*err)[i] > 0.1*(*area)[i] && ica << 905         if ((*err)[i] > 0.1*(*area)[i] && icase == 1) 
1009     {                                            906     {
1010       (*b)[i] = 0;                               907       (*b)[i] = 0;
1011       (*a)[i] = 0;                               908       (*a)[i] = 0;
1012       (*c)[i] = 1.;                              909       (*c)[i] = 1.;
1013       icase = 2;                                 910       icase = 2;
1014       reLoop = true;                             911       reLoop = true;
1015     }                                            912     }
1016       }while(reLoop);                            913       }while(reLoop);
1017     delete pdfi;                                 914     delete pdfi;
1018     delete pdfih;                                915     delete pdfih;
1019     delete sumi;                                 916     delete sumi;
1020   }                                              917   }
1021     }while(x->size() < np);                      918     }while(x->size() < np);
1022                                                  919 
1023   if (x->size() != np || a->size() != np ||   << 920   if (x->size() != np || a->size() != np || 
1024       err->size() != np || area->size() != np    921       err->size() != np || area->size() != np)
1025     {                                            922     {
1026       G4Exception("G4PenelopeRayleighModel::I    923       G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1027       "em2050",FatalException,                   924       "em2050",FatalException,
1028       "Problem in building the extended Table    925       "Problem in building the extended Table for Sampling: array dimensions do not match ");
1029     }                                            926     }
1030                                                  927 
1031   /******************************************    928   /*******************************************************************************
1032    Renormalization                               929    Renormalization
1033   *******************************************    930   ********************************************************************************/
1034   G4double ws = 0;                               931   G4double ws = 0;
1035   for (std::size_t i=0;i<np-1;++i)            << 932   for (size_t i=0;i<np-1;i++)
1036     ws += (*area)[i];                            933     ws += (*area)[i];
1037   ws = 1.0/ws;                                   934   ws = 1.0/ws;
1038   G4double errMax = 0;                           935   G4double errMax = 0;
1039   for (std::size_t i=0;i<np-1;++i)            << 936   for (size_t i=0;i<np-1;i++)
1040     {                                            937     {
1041       (*area)[i] *= ws;                          938       (*area)[i] *= ws;
1042       (*err)[i] *= ws;                           939       (*err)[i] *= ws;
1043       errMax = std::max(errMax,(*err)[i]);       940       errMax = std::max(errMax,(*err)[i]);
1044     }                                            941     }
1045                                                  942 
1046   //Vector with the normalized cumulative dis    943   //Vector with the normalized cumulative distribution
1047   G4DataVector* PAC = new G4DataVector();        944   G4DataVector* PAC = new G4DataVector();
1048   PAC->push_back(0.);                            945   PAC->push_back(0.);
1049   for (std::size_t i=0;i<np-1;++i)            << 946   for (size_t i=0;i<np-1;i++)
1050     {                                            947     {
1051       G4double previous = (*PAC)[i];             948       G4double previous = (*PAC)[i];
1052       PAC->push_back(previous+(*area)[i]);       949       PAC->push_back(previous+(*area)[i]);
1053     }                                            950     }
1054   (*PAC)[PAC->size()-1] = 1.;                    951   (*PAC)[PAC->size()-1] = 1.;
1055                                               << 952          
1056   /******************************************    953   /*******************************************************************************
1057   Pre-calculated limits for the initial binar    954   Pre-calculated limits for the initial binary search for subsequent sampling
1058   *******************************************    955   ********************************************************************************/
1059   std::vector<std::size_t> *ITTL = new std::v << 956 
1060   std::vector<std::size_t> *ITTU = new std::v << 957   //G4DataVector* ITTL = new G4DataVector();
                                                   >> 958   std::vector<size_t> *ITTL = new std::vector<size_t>;
                                                   >> 959   std::vector<size_t> *ITTU = new std::vector<size_t>;
1061                                                  960 
1062   //Just create place-holders                    961   //Just create place-holders
1063   for (std::size_t i=0;i<np;++i)              << 962   for (size_t i=0;i<np;i++)
1064     {                                            963     {
1065       ITTL->push_back(0);                        964       ITTL->push_back(0);
1066       ITTU->push_back(0);                        965       ITTU->push_back(0);
1067     }                                            966     }
1068                                                  967 
1069   G4double bin = 1.0/(np-1);                     968   G4double bin = 1.0/(np-1);
1070   (*ITTL)[0]=0;                                  969   (*ITTL)[0]=0;
1071   for (std::size_t i=1;i<(np-1);++i)          << 970   for (size_t i=1;i<(np-1);i++)
1072     {                                            971     {
1073       G4double ptst = i*bin;                  << 972       G4double ptst = i*bin; 
1074       G4bool found = false;                      973       G4bool found = false;
1075       for (std::size_t j=(*ITTL)[i-1];j<np && << 974       for (size_t j=(*ITTL)[i-1];j<np && !found;j++)
1076   {                                              975   {
1077     if ((*PAC)[j] > ptst)                        976     if ((*PAC)[j] > ptst)
1078       {                                          977       {
1079         (*ITTL)[i] = j-1;                        978         (*ITTL)[i] = j-1;
1080         (*ITTU)[i-1] = j;                        979         (*ITTU)[i-1] = j;
1081         found = true;                            980         found = true;
1082       }                                          981       }
1083   }                                              982   }
1084     }                                            983     }
1085   (*ITTU)[ITTU->size()-2] = ITTU->size()-1;      984   (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1086   (*ITTU)[ITTU->size()-1] = ITTU->size()-1;      985   (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1087   (*ITTL)[ITTL->size()-1] = ITTU->size()-2;      986   (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1088                                                  987 
1089   if (ITTU->size() != np || ITTU->size() != n    988   if (ITTU->size() != np || ITTU->size() != np)
1090     {                                            989     {
1091       G4Exception("G4PenelopeRayleighModel::I    990       G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1092       "em2051",FatalException,                   991       "em2051",FatalException,
1093       "Problem in building the Limit Tables f    992       "Problem in building the Limit Tables for Sampling: array dimensions do not match");
1094     }                                            993     }
1095                                                  994 
                                                   >> 995 
1096   /******************************************    996   /********************************************************************************
1097     Copy tables                                  997     Copy tables
1098   *******************************************    998   ********************************************************************************/
1099   G4PenelopeSamplingData* theTable = new G4Pe    999   G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np);
1100   for (std::size_t i=0;i<np;++i)              << 1000   for (size_t i=0;i<np;i++)
1101     {                                            1001     {
1102       theTable->AddPoint((*x)[i],(*PAC)[i],(*    1002       theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
1103     }                                            1003     }
1104                                                  1004 
1105   if (fVerboseLevel > 2)                      << 1005   if (verboseLevel > 2)
1106     {                                            1006     {
1107       G4cout << "**************************** << 1007       G4cout << "*************************************************************************" << 
1108   G4endl;                                        1008   G4endl;
1109       G4cout << "Sampling table for Penelope     1009       G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl;
1110       theTable->DumpTable();                     1010       theTable->DumpTable();
1111     }                                         << 1011     } 
1112   fSamplingTable->insert(std::make_pair(mat,t << 1012   samplingTable->insert(std::make_pair(mat,theTable));
1113                                                  1013 
                                                   >> 1014  
1114   //Clean up temporary vectors                   1015   //Clean up temporary vectors
1115   delete x;                                      1016   delete x;
1116   delete a;                                      1017   delete a;
1117   delete b;                                      1018   delete b;
1118   delete c;                                      1019   delete c;
1119   delete err;                                    1020   delete err;
1120   delete area;                                   1021   delete area;
1121   delete PAC;                                    1022   delete PAC;
1122   delete ITTL;                                   1023   delete ITTL;
1123   delete ITTU;                                   1024   delete ITTU;
1124                                                  1025 
1125   //DONE!                                        1026   //DONE!
1126   return;                                        1027   return;
                                                   >> 1028   
1127 }                                                1029 }
1128                                                  1030 
1129 //....oooOO0OOooo........oooOO0OOooo........o    1031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1130                                                  1032 
1131 void G4PenelopeRayleighModel::GetPMaxTable(co    1033 void G4PenelopeRayleighModel::GetPMaxTable(const G4Material* mat)
1132 {                                                1034 {
1133   if (!fPMaxTable)                            << 1035   if (!pMaxTable)
1134     {                                            1036     {
1135       G4cout << "G4PenelopeRayleighModel::Bui    1037       G4cout << "G4PenelopeRayleighModel::BuildPMaxTable" << G4endl;
1136       G4cout << "Going to instanziate the fPM << 1038       G4cout << "Going to instanziate the pMaxTable !" << G4endl;
1137       G4cout << "That should _not_ be here! "    1039       G4cout << "That should _not_ be here! " << G4endl;
1138       fPMaxTable = new std::map<const G4Mater << 1040       pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
1139     }                                            1041     }
1140   //check if the table is already there          1042   //check if the table is already there
1141   if (fPMaxTable->count(mat))                 << 1043   if (pMaxTable->count(mat))
1142     return;                                      1044     return;
1143                                                  1045 
1144   //otherwise build it                           1046   //otherwise build it
1145   if (!fSamplingTable)                        << 1047   if (!samplingTable)
1146     {                                            1048     {
1147       G4Exception("G4PenelopeRayleighModel::G    1049       G4Exception("G4PenelopeRayleighModel::GetPMaxTable()",
1148       "em2052",FatalException,                   1050       "em2052",FatalException,
1149       "SamplingTable is not properly instanti    1051       "SamplingTable is not properly instantiated");
1150       return;                                    1052       return;
1151     }                                            1053     }
1152                                                  1054 
1153   //This should not be: the sampling table is << 1055   if (!samplingTable->count(mat))
1154   if (!fSamplingTable->count(mat))            << 1056     InitializeSamplingAlgorithm(mat);      
1155     {                                         << 1057   
1156        G4ExceptionDescription ed;             << 1058   G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second;
1157        ed << "Sampling table for material " < << 1059   size_t tablePoints = theTable->GetNumberOfStoredPoints();
1158        G4Exception("G4PenelopeRayleighModel:: << 
1159                   "em2052",FatalException,    << 
1160                   ed);                        << 
1161        return;                                << 
1162     }                                         << 
1163                                                  1060 
1164   G4PenelopeSamplingData *theTable = fSamplin << 1061   size_t nOfEnergyPoints = logEnergyGridPMax.size();
1165   std::size_t tablePoints = theTable->GetNumb << 
1166                                               << 
1167   std::size_t nOfEnergyPoints = fLogEnergyGri << 
1168   G4PhysicsFreeVector* theVec = new G4Physics    1062   G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints);
1169                                                  1063 
1170   const std::size_t nip = 51; //hard-coded in << 1064   const size_t nip = 51; //hard-coded in Penelope
1171                                                  1065 
1172   for (std::size_t ie=0;ie<fLogEnergyGridPMax << 1066   for (size_t ie=0;ie<logEnergyGridPMax.size();ie++)
1173     {                                            1067     {
1174       G4double energy = G4Exp(fLogEnergyGridP << 1068       G4double energy = std::exp(logEnergyGridPMax[ie]);
1175       G4double Qm = 2.0*energy/electron_mass_    1069       G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now
1176       G4double Qm2 = Qm*Qm;                      1070       G4double Qm2 = Qm*Qm;
1177       G4double firstQ2 = theTable->GetX(0);      1071       G4double firstQ2 = theTable->GetX(0);
1178       G4double lastQ2 = theTable->GetX(tableP    1072       G4double lastQ2 = theTable->GetX(tablePoints-1);
1179       G4double thePMax = 0;                      1073       G4double thePMax = 0;
1180                                               << 1074       
1181       if (Qm2 > firstQ2)                         1075       if (Qm2 > firstQ2)
1182   {                                              1076   {
1183     if (Qm2 < lastQ2)                            1077     if (Qm2 < lastQ2)
1184       {                                          1078       {
1185         //bisection to look for the index of     1079         //bisection to look for the index of Qm
1186         std::size_t lowerBound = 0;           << 1080         size_t lowerBound = 0;
1187         std::size_t upperBound = tablePoints- << 1081         size_t upperBound = tablePoints-1;
1188         while (lowerBound <= upperBound)         1082         while (lowerBound <= upperBound)
1189     {                                            1083     {
1190       std::size_t midBin = (lowerBound + uppe << 1084       size_t midBin = (lowerBound + upperBound)/2;
1191       if( Qm2 < theTable->GetX(midBin))          1085       if( Qm2 < theTable->GetX(midBin))
1192         { upperBound = midBin-1; }               1086         { upperBound = midBin-1; }
1193       else                                       1087       else
1194         { lowerBound = midBin+1; }               1088         { lowerBound = midBin+1; }
1195     }                                            1089     }
1196         //upperBound is the output (but also     1090         //upperBound is the output (but also lowerBounf --> should be the same!)
1197         G4double Q1 = theTable->GetX(upperBou    1091         G4double Q1 = theTable->GetX(upperBound);
1198         G4double Q2 = Qm2;                       1092         G4double Q2 = Qm2;
1199         G4double DQ = (Q2-Q1)/((G4double)(nip    1093         G4double DQ = (Q2-Q1)/((G4double)(nip-1));
1200         G4double theA = theTable->GetA(upperB    1094         G4double theA = theTable->GetA(upperBound);
1201         G4double theB = theTable->GetB(upperB    1095         G4double theB = theTable->GetB(upperBound);
1202         G4double thePAC = theTable->GetPAC(up    1096         G4double thePAC = theTable->GetPAC(upperBound);
1203         G4DataVector* fun = new G4DataVector(    1097         G4DataVector* fun = new G4DataVector();
1204         for (std::size_t k=0;k<nip;++k)       << 1098         for (size_t k=0;k<nip;k++)
1205     {                                            1099     {
1206       G4double qi = Q1 + k*DQ;                   1100       G4double qi = Q1 + k*DQ;
1207       G4double tau = (qi-Q1)/                    1101       G4double tau = (qi-Q1)/
1208         (theTable->GetX(upperBound+1)-Q1);       1102         (theTable->GetX(upperBound+1)-Q1);
1209       G4double con1 = 2.0*theB*tau;           << 1103       G4double con1 = 2.0*theB*tau; 
1210       G4double ci = 1.0+theA+theB;               1104       G4double ci = 1.0+theA+theB;
1211       G4double con2 = ci-theA*tau;               1105       G4double con2 = ci-theA*tau;
1212       G4double etap = 0;                         1106       G4double etap = 0;
1213       if (std::fabs(con1) > 1.0e-16*std::fabs    1107       if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
1214         etap = con2*(1.0-std::sqrt(1.0-2.0*ta    1108         etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
1215       else                                       1109       else
1216         etap = tau/con2;                         1110         etap = tau/con2;
1217       G4double theFun = (theTable->GetPAC(upp    1111       G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)*
1218         (1.0+(theA+theB*etap)*etap)*(1.0+(the    1112         (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
1219         ((1.0-theB*etap*etap)*ci*(theTable->G    1113         ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1));
1220       fun->push_back(theFun);                    1114       fun->push_back(theFun);
1221     }                                            1115     }
1222         //Now intergrate numerically the fun     1116         //Now intergrate numerically the fun Cavalieri-Simpson's method
1223         G4DataVector* sum = new G4DataVector;    1117         G4DataVector* sum = new G4DataVector;
1224         G4double CONS = DQ*(1./12.);             1118         G4double CONS = DQ*(1./12.);
1225         G4double HCONS = 0.5*CONS;               1119         G4double HCONS = 0.5*CONS;
1226         sum->push_back(0.);                      1120         sum->push_back(0.);
1227         G4double secondPoint = (*sum)[0] +    << 1121         G4double secondPoint = (*sum)[0] + 
1228     (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*C    1122     (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
1229         sum->push_back(secondPoint);             1123         sum->push_back(secondPoint);
1230         for (std::size_t hh=2;hh<nip-1;++hh)  << 1124         for (size_t hh=2;hh<nip-1;hh++)
1231     {                                            1125     {
1232       G4double previous = (*sum)[hh-1];          1126       G4double previous = (*sum)[hh-1];
1233       G4double next = previous+(13.0*((*fun)[    1127       G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
1234               (*fun)[hh+1]-(*fun)[hh-2])*HCON    1128               (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
1235       sum->push_back(next);                      1129       sum->push_back(next);
1236     }                                            1130     }
1237         G4double last = (*sum)[nip-2]+(5.0*(*    1131         G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
1238                (*fun)[nip-3])*CONS;              1132                (*fun)[nip-3])*CONS;
1239         sum->push_back(last);                 << 1133         sum->push_back(last);  
1240         thePMax = thePAC + (*sum)[sum->size()    1134         thePMax = thePAC + (*sum)[sum->size()-1]; //last point
1241         delete fun;                              1135         delete fun;
1242         delete sum;                              1136         delete sum;
1243       }                                          1137       }
1244     else                                         1138     else
1245       {                                          1139       {
1246         thePMax = 1.0;                           1140         thePMax = 1.0;
1247       }                                       << 1141       }  
1248   }                                              1142   }
1249       else                                       1143       else
1250   {                                              1144   {
1251     thePMax = theTable->GetPAC(0);               1145     thePMax = theTable->GetPAC(0);
1252   }                                              1146   }
1253                                                  1147 
1254       //Write number in the table                1148       //Write number in the table
1255       theVec->PutValue(ie,energy,thePMax);       1149       theVec->PutValue(ie,energy,thePMax);
1256   }                                              1150   }
1257                                               << 1151   
1258   fPMaxTable->insert(std::make_pair(mat,theVe << 1152   pMaxTable->insert(std::make_pair(mat,theVec));
1259   return;                                        1153   return;
                                                   >> 1154 
1260 }                                                1155 }
1261                                                  1156 
1262 //....oooOO0OOooo........oooOO0OOooo........o    1157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1263                                                  1158 
1264 void G4PenelopeRayleighModel::DumpFormFactorT    1159 void G4PenelopeRayleighModel::DumpFormFactorTable(const G4Material* mat)
1265 {                                                1160 {
1266   G4cout << "********************************    1161   G4cout << "*****************************************************************" << G4endl;
1267   G4cout << "G4PenelopeRayleighModel: Form Fa    1162   G4cout << "G4PenelopeRayleighModel: Form Factor Table for " << mat->GetName() << G4endl;
1268   //try to use the same format as Penelope-Fo    1163   //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
1269   G4cout <<  "Q/(m_e*c)                 F(Q)     1164   G4cout <<  "Q/(m_e*c)                 F(Q)     " << G4endl;
1270   G4cout << "********************************    1165   G4cout << "*****************************************************************" << G4endl;
1271   if (!fLogFormFactorTable->count(mat))       << 1166   if (!logFormFactorTable->count(mat))
1272     BuildFormFactorTable(mat);                   1167     BuildFormFactorTable(mat);
1273                                               << 1168   
1274   G4PhysicsFreeVector* theVec = fLogFormFacto << 1169   G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
1275   for (std::size_t i=0;i<theVec->GetVectorLen << 1170   for (size_t i=0;i<theVec->GetVectorLength();i++)
1276     {                                            1171     {
1277       G4double logQ2 = theVec->GetLowEdgeEner    1172       G4double logQ2 = theVec->GetLowEdgeEnergy(i);
1278       G4double Q = G4Exp(0.5*logQ2);          << 1173       G4double Q = std::exp(0.5*logQ2);
1279       G4double logF2 = (*theVec)[i];             1174       G4double logF2 = (*theVec)[i];
1280       G4double F = G4Exp(0.5*logF2);          << 1175       G4double F = std::exp(0.5*logF2);
1281       G4cout << Q << "              " << F <<    1176       G4cout << Q << "              " << F << G4endl;
1282     }                                            1177     }
1283   //DONE                                         1178   //DONE
1284   return;                                        1179   return;
1285 }                                             << 
1286                                               << 
1287 //....oooOO0OOooo........oooOO0OOooo........o << 
1288                                               << 
1289 void G4PenelopeRayleighModel::SetParticle(con << 
1290 {                                             << 
1291   if(!fParticle) {                            << 
1292     fParticle = p;                            << 
1293   }                                           << 
1294 }                                                1180 }
1295                                                  1181