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


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