Geant4 Cross Reference

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


  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: G4PenelopeCrossSection.cc 95950 2016-03-03 10:42:48Z gcosmo $
 26 //                                                 27 //
 27 // Author: Luciano Pandola                         28 // Author: Luciano Pandola
 28 //                                                 29 //
 29 // History:                                        30 // History:
 30 // --------                                        31 // --------
 31 // 18 Mar 2010   L Pandola    First implementa     32 // 18 Mar 2010   L Pandola    First implementation
 32 // 09 Mar 2012   L. Pandola   Add public metho     33 // 09 Mar 2012   L. Pandola   Add public method (and machinery) to return
 33 //                            the absolute and     34 //                            the absolute and the normalized shell cross
 34 //                            sections indepen     35 //                            sections independently.
 35 //                                                 36 //
 36 #include "G4PenelopeCrossSection.hh"               37 #include "G4PenelopeCrossSection.hh"
 37 #include "G4SystemOfUnits.hh"                      38 #include "G4SystemOfUnits.hh"
 38 #include "G4PhysicsTable.hh"                       39 #include "G4PhysicsTable.hh"
 39 #include "G4PhysicsFreeVector.hh"                  40 #include "G4PhysicsFreeVector.hh"
 40 #include "G4DataVector.hh"                         41 #include "G4DataVector.hh"
 41 #include "G4Exp.hh"                                42 #include "G4Exp.hh"
 42 #include "G4Log.hh"                            << 
 43                                                    43 
 44 //....oooOO0OOooo........oooOO0OOooo........oo     44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
 45 G4PenelopeCrossSection::G4PenelopeCrossSection     45 G4PenelopeCrossSection::G4PenelopeCrossSection(size_t nPointsE,size_t nShells) :
 46   fSoftCrossSections(nullptr),                 <<  46   numberOfEnergyPoints(nPointsE),numberOfShells(nShells),softCrossSections(0),
 47   fHardCrossSections(nullptr),fShellCrossSecti <<  47   hardCrossSections(0),shellCrossSections(0),shellNormalizedCrossSections(0)
 48   fShellNormalizedCrossSections(nullptr),      << 
 49   fNumberOfEnergyPoints(nPointsE),fNumberOfShe << 
 50 {                                                  48 {
 51   //check the number of points is not zero         49   //check the number of points is not zero
 52   if (!fNumberOfEnergyPoints)                  <<  50   if (!numberOfEnergyPoints)
 53     {                                              51     {
 54       G4ExceptionDescription ed;                   52       G4ExceptionDescription ed;
 55       ed << "G4PenelopeCrossSection: invalid n     53       ed << "G4PenelopeCrossSection: invalid number of energy points " << G4endl;
 56       G4Exception("G4PenelopeCrossSection::G4P     54       G4Exception("G4PenelopeCrossSection::G4PenelopeCrossSection()",
 57       "em2017",FatalException,ed);                 55       "em2017",FatalException,ed);
 58     }                                              56     }
 59                                                    57 
 60   fIsNormalized = false;                       <<  58   isNormalized = false;
 61                                                    59 
 62   // 1) soft XS table                              60   // 1) soft XS table
 63   fSoftCrossSections = new G4PhysicsTable();   <<  61   softCrossSections = new G4PhysicsTable();
 64   //the table contains 3 G4PhysicsFreeVectors,     62   //the table contains 3 G4PhysicsFreeVectors,
 65   //(fSoftCrossSections)[0] -->  log XS0 vs. l <<  63   //(softCrossSections)[0] -->  log XS0 vs. log E
 66   //(fSoftCrossSections)[1] -->  log XS1 vs. l <<  64   //(softCrossSections)[1] -->  log XS1 vs. log E
 67   //(fSoftCrossSections)[2] -->  log XS2 vs. l <<  65   //(softCrossSections)[2] -->  log XS2 vs. log E
 68                                                    66 
 69   //everything is log-log                          67   //everything is log-log
 70   for (size_t i=0;i<3;i++)                         68   for (size_t i=0;i<3;i++)
 71     fSoftCrossSections->push_back(new G4Physic <<  69     softCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
 72                                                    70 
 73   //2) hard XS table                               71   //2) hard XS table
 74   fHardCrossSections = new G4PhysicsTable();   <<  72   hardCrossSections = new G4PhysicsTable();
 75   //the table contains 3 G4PhysicsFreeVectors,     73   //the table contains 3 G4PhysicsFreeVectors,
 76   //(fHardCrossSections)[0] -->  log XH0 vs. l <<  74   //(hardCrossSections)[0] -->  log XH0 vs. log E
 77   //(fHardCrossSections)[1] -->  log XH1 vs. l <<  75   //(hardCrossSections)[1] -->  log XH1 vs. log E
 78   //(fHardCrossSections)[2] -->  log XH2 vs. l <<  76   //(hardCrossSections)[2] -->  log XH2 vs. log E
 79                                                    77 
 80   //everything is log-log                          78   //everything is log-log
 81   for (size_t i=0;i<3;i++)                         79   for (size_t i=0;i<3;i++)
 82     fHardCrossSections->push_back(new G4Physic <<  80     hardCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
 83                                                    81 
 84   //3) shell XS table, if it is the case           82   //3) shell XS table, if it is the case
 85   if (fNumberOfShells)                         <<  83   if (numberOfShells)
 86     {                                              84     {
 87       fShellCrossSections = new G4PhysicsTable <<  85       shellCrossSections = new G4PhysicsTable();
 88       fShellNormalizedCrossSections = new G4Ph <<  86       shellNormalizedCrossSections = new G4PhysicsTable();
 89       //the table has to contain numberofShell     87       //the table has to contain numberofShells G4PhysicsFreeVectors,
 90       //(theTable)[ishell] --> cross section f     88       //(theTable)[ishell] --> cross section for shell #ishell
 91       for (size_t i=0;i<fNumberOfShells;i++)   <<  89       for (size_t i=0;i<numberOfShells;i++)
 92   {                                                90   {
 93     fShellCrossSections->push_back(new G4Physi <<  91     shellCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
 94     fShellNormalizedCrossSections->push_back(n <<  92     shellNormalizedCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
 95   }                                                93   }
 96     }                                              94     }
 97 }                                                  95 }
 98                                                    96 
 99 //....oooOO0OOooo........oooOO0OOooo........oo     97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
100 G4PenelopeCrossSection::~G4PenelopeCrossSectio     98 G4PenelopeCrossSection::~G4PenelopeCrossSection()
101 {                                                  99 {
102   //clean up tables                               100   //clean up tables
103   if (fShellCrossSections)                     << 101   if (shellCrossSections)
104     {                                             102     {
105       fShellCrossSections->clearAndDestroy();  << 103       //shellCrossSections->clearAndDestroy();
106       delete fShellCrossSections;              << 104       delete shellCrossSections;
107     }                                             105     }
108   if (fShellNormalizedCrossSections)           << 106   if (shellNormalizedCrossSections)
109     {                                             107     {
110       fShellNormalizedCrossSections->clearAndD << 108       //shellNormalizedCrossSections->clearAndDestroy();
111       delete fShellNormalizedCrossSections;    << 109       delete shellNormalizedCrossSections;
112     }                                             110     }
113   if (fSoftCrossSections)                      << 111   if (softCrossSections)
114     {                                             112     {
115       fSoftCrossSections->clearAndDestroy();   << 113       //softCrossSections->clearAndDestroy();
116       delete fSoftCrossSections;               << 114       delete softCrossSections;
117     }                                             115     }
118   if (fHardCrossSections)                      << 116   if (hardCrossSections)
119     {                                             117     {
120       fHardCrossSections->clearAndDestroy();   << 118       //hardCrossSections->clearAndDestroy();
121       delete fHardCrossSections;               << 119       delete hardCrossSections;
122     }                                             120     }
123 }                                                 121 }
124                                                   122 
125 //....oooOO0OOooo........oooOO0OOooo........oo    123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
126 void G4PenelopeCrossSection::AddCrossSectionPo    124 void G4PenelopeCrossSection::AddCrossSectionPoint(size_t binNumber,G4double energy,
127               G4double XH0,                       125               G4double XH0,
128               G4double XH1, G4double XH2,         126               G4double XH1, G4double XH2,
129               G4double XS0, G4double XS1,         127               G4double XS0, G4double XS1,
130               G4double XS2)                       128               G4double XS2)
131 {                                                 129 {
132   if (!fSoftCrossSections || !fHardCrossSectio << 130   if (!softCrossSections || !hardCrossSections)
133     {                                             131     {
134       G4cout << "Something wrong in G4Penelope    132       G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
135   G4endl;                                         133   G4endl;
136       G4cout << "Trying to fill un-initialized    134       G4cout << "Trying to fill un-initialized tables" << G4endl;
137       return;                                     135       return;
138     }                                             136     }
139                                                   137 
140   //fill vectors                                  138   //fill vectors
141   G4PhysicsFreeVector* theVector = (G4PhysicsF << 139   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[0];
142                                                   140 
143   if (binNumber >= fNumberOfEnergyPoints)      << 141   if (binNumber >= numberOfEnergyPoints)
144     {                                             142     {
145       G4cout << "Something wrong in G4Penelope    143       G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
146   G4endl;                                         144   G4endl;
147       G4cout << "Trying to register more point    145       G4cout << "Trying to register more points than originally declared" << G4endl;
148       return;                                     146       return;
149     }                                             147     }
150    G4double logEne = G4Log(energy);            << 148    G4double logEne = std::log(energy);
151                                                   149 
152    //XS0                                          150    //XS0
153    G4double val = G4Log(std::max(XS0,1e-42*cm2 << 151    G4double val = std::log(std::max(XS0,1e-42*cm2)); //avoid log(0)
154    theVector->PutValues(binNumber,logEne,val); << 152    theVector->PutValue(binNumber,logEne,val);
155                                                   153 
156    //XS1                                          154    //XS1
157    theVector = (G4PhysicsFreeVector*) (*fSoftC << 155    theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
158    val =  G4Log(std::max(XS1,1e-42*eV*cm2)); / << 156    val =  std::log(std::max(XS1,1e-42*eV*cm2)); //avoid log(0)
159    theVector->PutValues(binNumber,logEne,val); << 157    theVector->PutValue(binNumber,logEne,val);
160                                                   158 
161    //XS2                                          159    //XS2
162    theVector = (G4PhysicsFreeVector*) (*fSoftC << 160    theVector = (G4PhysicsFreeVector*) (*softCrossSections)[2];
163    val =  G4Log(std::max(XS2,1e-42*eV*eV*cm2)) << 161    val =  std::log(std::max(XS2,1e-42*eV*eV*cm2)); //avoid log(0)
164    theVector->PutValues(binNumber,logEne,val); << 162    theVector->PutValue(binNumber,logEne,val);
165                                                   163 
166    //XH0                                          164    //XH0
167    theVector = (G4PhysicsFreeVector*) (*fHardC << 165    theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
168    val =  G4Log(std::max(XH0,1e-42*cm2)); //av << 166    val =  std::log(std::max(XH0,1e-42*cm2)); //avoid log(0)
169    theVector->PutValues(binNumber,logEne,val); << 167    theVector->PutValue(binNumber,logEne,val);
170                                                   168 
171    //XH1                                          169    //XH1
172    theVector = (G4PhysicsFreeVector*) (*fHardC << 170    theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[1];
173    val =  G4Log(std::max(XH1,1e-42*eV*cm2)); / << 171    val =  std::log(std::max(XH1,1e-42*eV*cm2)); //avoid log(0)
174    theVector->PutValues(binNumber,logEne,val); << 172    theVector->PutValue(binNumber,logEne,val);
175                                                   173 
176     //XH2                                         174     //XH2
177    theVector = (G4PhysicsFreeVector*) (*fHardC << 175    theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[2];
178    val =  G4Log(std::max(XH2,1e-42*eV*eV*cm2)) << 176    val =  std::log(std::max(XH2,1e-42*eV*eV*cm2)); //avoid log(0)
179    theVector->PutValues(binNumber,logEne,val); << 177    theVector->PutValue(binNumber,logEne,val);
180                                                   178 
181    return;                                        179    return;
182 }                                                 180 }
183                                                   181 
184 //....oooOO0OOooo........oooOO0OOooo........oo    182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
185                                                   183 
186 void G4PenelopeCrossSection::AddShellCrossSect    184 void G4PenelopeCrossSection::AddShellCrossSectionPoint(size_t binNumber,
187                    size_t shellID,                185                    size_t shellID,
188                    G4double energy,               186                    G4double energy,
189                    G4double xs)                   187                    G4double xs)
190 {                                                 188 {
191   if (!fShellCrossSections)                    << 189   if (!shellCrossSections)
192     {                                             190     {
193       G4cout << "Something wrong in G4Penelope    191       G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
194   G4endl;                                         192   G4endl;
195       G4cout << "Trying to fill un-initialized    193       G4cout << "Trying to fill un-initialized table" << G4endl;
196       return;                                     194       return;
197     }                                             195     }
198                                                   196 
199   if (shellID >= fNumberOfShells)              << 197   if (shellID >= numberOfShells)
200     {                                             198     {
201       G4cout << "Something wrong in G4Penelope    199       G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
202   G4endl;                                         200   G4endl;
203       G4cout << "Trying to fill shell #" << sh    201       G4cout << "Trying to fill shell #" << shellID << " while the maximum is "
204        <<  fNumberOfShells-1 << G4endl;        << 202        <<  numberOfShells-1 << G4endl;
205       return;                                     203       return;
206     }                                             204     }
207                                                   205 
208   //fill vector                                   206   //fill vector
209   G4PhysicsFreeVector* theVector = (G4PhysicsF << 207   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
210                                                   208 
211   if (binNumber >= fNumberOfEnergyPoints)      << 209   if (binNumber >= numberOfEnergyPoints)
212     {                                             210     {
213       G4cout << "Something wrong in G4Penelope    211       G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
214   G4endl;                                         212   G4endl;
215       G4cout << "Trying to register more point    213       G4cout << "Trying to register more points than originally declared" << G4endl;
216       return;                                     214       return;
217     }                                             215     }
218    G4double logEne = G4Log(energy);            << 216    G4double logEne = std::log(energy);
219    G4double val = G4Log(std::max(xs,1e-42*cm2) << 217    G4double val = std::log(std::max(xs,1e-42*cm2)); //avoid log(0)
220    theVector->PutValues(binNumber,logEne,val); << 218    theVector->PutValue(binNumber,logEne,val);
221                                                   219 
222    return;                                        220    return;
223 }                                                 221 }
224                                                   222 
225 //....oooOO0OOooo........oooOO0OOooo........oo    223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
226                                                   224 
227 G4double G4PenelopeCrossSection::GetTotalCross    225 G4double G4PenelopeCrossSection::GetTotalCrossSection(G4double energy) const
228 {                                                 226 {
229   G4double result = 0;                            227   G4double result = 0;
230   //take here XS0 + XH0                           228   //take here XS0 + XH0
231   if (!fSoftCrossSections || !fHardCrossSectio << 229   if (!softCrossSections || !hardCrossSections)
232     {                                             230     {
233       G4cout << "Something wrong in G4Penelope    231       G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
234   G4endl;                                         232   G4endl;
235       G4cout << "Trying to retrieve from un-in    233       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
236       return result;                              234       return result;
237     }                                             235     }
238                                                   236 
239   // 1) soft part                                 237   // 1) soft part
240   G4PhysicsFreeVector* theVector = (G4PhysicsF << 238   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[0];
241   if (theVector->GetVectorLength() < fNumberOf << 239   if (theVector->GetVectorLength() < numberOfEnergyPoints)
242     {                                             240     {
243       G4cout << "Something wrong in G4Penelope    241       G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
244   G4endl;                                         242   G4endl;
245       G4cout << "Soft cross section table look    243       G4cout << "Soft cross section table looks not filled" << G4endl;
246       return result;                              244       return result;
247     }                                             245     }
248   G4double logene = G4Log(energy);             << 246   G4double logene = std::log(energy);
249   G4double logXS = theVector->Value(logene);      247   G4double logXS = theVector->Value(logene);
250   G4double softXS = G4Exp(logXS);                 248   G4double softXS = G4Exp(logXS);
251                                                   249 
252    // 2) hard part                                250    // 2) hard part
253   theVector = (G4PhysicsFreeVector*) (*fHardCr << 251   theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
254   if (theVector->GetVectorLength() < fNumberOf << 252   if (theVector->GetVectorLength() < numberOfEnergyPoints)
255     {                                             253     {
256       G4cout << "Something wrong in G4Penelope    254       G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
257   G4endl;                                         255   G4endl;
258       G4cout << "Hard cross section table look    256       G4cout << "Hard cross section table looks not filled" << G4endl;
259       return result;                              257       return result;
260     }                                             258     }
261   logXS = theVector->Value(logene);               259   logXS = theVector->Value(logene);
262   G4double hardXS = G4Exp(logXS);                 260   G4double hardXS = G4Exp(logXS);
263                                                   261 
264   result = hardXS + softXS;                       262   result = hardXS + softXS;
265   return result;                                  263   return result;
                                                   >> 264 
266 }                                                 265 }
267                                                   266 
268 //....oooOO0OOooo........oooOO0OOooo........oo    267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
269                                                   268 
270 G4double G4PenelopeCrossSection::GetHardCrossS    269 G4double G4PenelopeCrossSection::GetHardCrossSection(G4double energy) const
271 {                                                 270 {
272   G4double result = 0;                            271   G4double result = 0;
273   //take here XH0                                 272   //take here XH0
274   if (!fHardCrossSections)                     << 273   if (!hardCrossSections)
275     {                                             274     {
276       G4cout << "Something wrong in G4Penelope    275       G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
277   G4endl;                                         276   G4endl;
278       G4cout << "Trying to retrieve from un-in    277       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
279       return result;                              278       return result;
280     }                                             279     }
281                                                   280 
282   G4PhysicsFreeVector* theVector = (G4PhysicsF << 281   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
283   if (theVector->GetVectorLength() < fNumberOf << 282   if (theVector->GetVectorLength() < numberOfEnergyPoints)
284     {                                             283     {
285       G4cout << "Something wrong in G4Penelope    284       G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
286   G4endl;                                         285   G4endl;
287       G4cout << "Hard cross section table look    286       G4cout << "Hard cross section table looks not filled" << G4endl;
288       return result;                              287       return result;
289     }                                             288     }
290   G4double logene = G4Log(energy);             << 289   G4double logene = std::log(energy);
291   G4double logXS = theVector->Value(logene);      290   G4double logXS = theVector->Value(logene);
292   result = G4Exp(logXS);                          291   result = G4Exp(logXS);
293                                                   292 
294   return result;                                  293   return result;
295 }                                                 294 }
296                                                   295 
                                                   >> 296 
297 //....oooOO0OOooo........oooOO0OOooo........oo    297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
298                                                   298 
299 G4double G4PenelopeCrossSection::GetSoftStoppi    299 G4double G4PenelopeCrossSection::GetSoftStoppingPower(G4double energy) const
300 {                                                 300 {
301   G4double result = 0;                            301   G4double result = 0;
302   //take here XH0                                 302   //take here XH0
303   if (!fSoftCrossSections)                     << 303   if (!softCrossSections)
304     {                                             304     {
305       G4cout << "Something wrong in G4Penelope    305       G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
306   G4endl;                                         306   G4endl;
307       G4cout << "Trying to retrieve from un-in    307       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
308       return result;                              308       return result;
309     }                                             309     }
310                                                   310 
311   G4PhysicsFreeVector* theVector = (G4PhysicsF << 311   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
312   if (theVector->GetVectorLength() < fNumberOf << 312   if (theVector->GetVectorLength() < numberOfEnergyPoints)
313     {                                             313     {
314       G4cout << "Something wrong in G4Penelope    314       G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
315   G4endl;                                         315   G4endl;
316       G4cout << "Soft cross section table look    316       G4cout << "Soft cross section table looks not filled" << G4endl;
317       return result;                              317       return result;
318     }                                             318     }
319   G4double logene = G4Log(energy);             << 319   G4double logene = std::log(energy);
320   G4double logXS = theVector->Value(logene);      320   G4double logXS = theVector->Value(logene);
321   result = G4Exp(logXS);                          321   result = G4Exp(logXS);
322                                                   322 
323   return result;                                  323   return result;
324 }                                                 324 }
325                                                   325 
326 //....oooOO0OOooo........oooOO0OOooo........oo    326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
327                                                   327 
328 G4double G4PenelopeCrossSection::GetShellCross    328 G4double G4PenelopeCrossSection::GetShellCrossSection(size_t shellID,G4double energy) const
329 {                                                 329 {
330   G4double result = 0;                            330   G4double result = 0;
331   if (!fShellCrossSections)                    << 331   if (!shellCrossSections)
332     {                                             332     {
333       G4cout << "Something wrong in G4Penelope    333       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
334   G4endl;                                         334   G4endl;
335       G4cout << "Trying to retrieve from un-in    335       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
336       return result;                              336       return result;
337     }                                             337     }
338   if (shellID >= fNumberOfShells)              << 338   if (shellID >= numberOfShells)
339     {                                             339     {
340       G4cout << "Something wrong in G4Penelope    340       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
341   G4endl;                                         341   G4endl;
342       G4cout << "Trying to retrieve shell #" <    342       G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
343        <<  fNumberOfShells-1 << G4endl;        << 343        <<  numberOfShells-1 << G4endl;
344       return result;                              344       return result;
345     }                                             345     }
346                                                   346 
347   G4PhysicsFreeVector* theVector = (G4PhysicsF << 347   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
348                                                   348 
349   if (theVector->GetVectorLength() < fNumberOf << 349   if (theVector->GetVectorLength() < numberOfEnergyPoints)
350     {                                             350     {
351       G4cout << "Something wrong in G4Penelope    351       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
352   G4endl;                                         352   G4endl;
353       G4cout << "Shell cross section table loo    353       G4cout << "Shell cross section table looks not filled" << G4endl;
354       return result;                              354       return result;
355     }                                             355     }
356   G4double logene = G4Log(energy);             << 356   G4double logene = std::log(energy);
357   G4double logXS = theVector->Value(logene);      357   G4double logXS = theVector->Value(logene);
358   result = G4Exp(logXS);                          358   result = G4Exp(logXS);
359                                                   359 
360   return result;                                  360   return result;
361 }                                                 361 }
362 //....oooOO0OOooo........oooOO0OOooo........oo    362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
363                                                   363 
364 G4double G4PenelopeCrossSection::GetNormalized    364 G4double G4PenelopeCrossSection::GetNormalizedShellCrossSection(size_t shellID,G4double energy) const
365 {                                                 365 {
366   G4double result = 0;                            366   G4double result = 0;
367   if (!fShellNormalizedCrossSections)          << 367   if (!shellNormalizedCrossSections)
368     {                                             368     {
369       G4cout << "Something wrong in G4Penelope    369       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
370   G4endl;                                         370   G4endl;
371       G4cout << "Trying to retrieve from un-in    371       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
372       return result;                              372       return result;
373     }                                             373     }
374                                                   374 
375   if (!fIsNormalized)                          << 375   if (!isNormalized)
376     {                                             376     {
377       G4cout << "Something wrong in G4Penelope    377       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" << G4endl;
378       G4cout << "The table of normalized cross    378       G4cout << "The table of normalized cross section is not initialized" << G4endl;
379     }                                             379     }
380                                                   380 
381   if (shellID >= fNumberOfShells)              << 381 
                                                   >> 382   if (shellID >= numberOfShells)
382     {                                             383     {
383       G4cout << "Something wrong in G4Penelope    384       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
384   G4endl;                                         385   G4endl;
385       G4cout << "Trying to retrieve shell #" <    386       G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
386        <<  fNumberOfShells-1 << G4endl;        << 387        <<  numberOfShells-1 << G4endl;
387       return result;                              388       return result;
388     }                                             389     }
389                                                   390 
390   const G4PhysicsFreeVector* theVector =          391   const G4PhysicsFreeVector* theVector =
391     (G4PhysicsFreeVector*) (*fShellNormalizedC << 392     (G4PhysicsFreeVector*) (*shellNormalizedCrossSections)[shellID];
392                                                   393 
393   if (theVector->GetVectorLength() < fNumberOf << 394   if (theVector->GetVectorLength() < numberOfEnergyPoints)
394     {                                             395     {
395       G4cout << "Something wrong in G4Penelope    396       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
396   G4endl;                                         397   G4endl;
397       G4cout << "Shell cross section table loo    398       G4cout << "Shell cross section table looks not filled" << G4endl;
398       return result;                              399       return result;
399     }                                             400     }
400   G4double logene = G4Log(energy);             << 401   G4double logene = std::log(energy);
401   G4double logXS = theVector->Value(logene);      402   G4double logXS = theVector->Value(logene);
402   result = G4Exp(logXS);                          403   result = G4Exp(logXS);
403                                                   404 
404   return result;                                  405   return result;
405 }                                                 406 }
406                                                   407 
                                                   >> 408 
                                                   >> 409 
407 //....oooOO0OOooo........oooOO0OOooo........oo    410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
408                                                   411 
409 void G4PenelopeCrossSection::NormalizeShellCro    412 void G4PenelopeCrossSection::NormalizeShellCrossSections()
410 {                                                 413 {
411   if (fIsNormalized) //already done!           << 414   if (isNormalized) //already done!
412     {                                             415     {
413       G4cout << "G4PenelopeCrossSection::Norma    416       G4cout << "G4PenelopeCrossSection::NormalizeShellCrossSections()" << G4endl;
414       G4cout << "already invoked. Ignore it" <    417       G4cout << "already invoked. Ignore it" << G4endl;
415       return;                                     418       return;
416     }                                             419     }
417                                                   420 
418   if (!fShellNormalizedCrossSections)          << 421   if (!shellNormalizedCrossSections)
419     {                                             422     {
420       G4cout << "Something wrong in G4Penelope    423       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
421   G4endl;                                         424   G4endl;
422       G4cout << "Trying to retrieve from un-in    425       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
423       return;                                     426       return;
424     }                                             427     }
425                                                   428 
426   for (size_t i=0;i<fNumberOfEnergyPoints;i++) << 429   for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
427     {                                             430     {
428       //energy grid is the same for all shells    431       //energy grid is the same for all shells
429                                                   432 
430       //Recalculate manually the XS factor, to    433       //Recalculate manually the XS factor, to avoid problems with
431       //underflows                                434       //underflows
432       G4double normFactor = 0.;                   435       G4double normFactor = 0.;
433       for (size_t shellID=0;shellID<fNumberOfS << 436       for (size_t shellID=0;shellID<numberOfShells;shellID++)
434   {                                               437   {
435     G4PhysicsFreeVector* theVec =                 438     G4PhysicsFreeVector* theVec =
436       (G4PhysicsFreeVector*) (*fShellCrossSect << 439       (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
437                                                   440 
438     normFactor += G4Exp((*theVec)[i]);            441     normFactor += G4Exp((*theVec)[i]);
439   }                                               442   }
440       G4double logNormFactor = G4Log(normFacto << 443       G4double logNormFactor = std::log(normFactor);
441       //Normalize                                 444       //Normalize
442       for (size_t shellID=0;shellID<fNumberOfS << 445       for (size_t shellID=0;shellID<numberOfShells;shellID++)
443   {                                               446   {
444    G4PhysicsFreeVector* theVec =                  447    G4PhysicsFreeVector* theVec =
445       (G4PhysicsFreeVector*) (*fShellNormalize << 448       (G4PhysicsFreeVector*) (*shellNormalizedCrossSections)[shellID];
446    G4PhysicsFreeVector* theFullVec =              449    G4PhysicsFreeVector* theFullVec =
447      (G4PhysicsFreeVector*) (*fShellCrossSecti << 450      (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
448    G4double previousValue = (*theFullVec)[i];     451    G4double previousValue = (*theFullVec)[i]; //log(XS)
449    G4double logEnergy = theFullVec->GetLowEdge    452    G4double logEnergy = theFullVec->GetLowEdgeEnergy(i);
450    //log(XS/normFactor) = log(XS) - log(normFa    453    //log(XS/normFactor) = log(XS) - log(normFactor)
451    theVec->PutValues(i,logEnergy,previousValue << 454    theVec->PutValue(i,logEnergy,previousValue-logNormFactor);
                                                   >> 455   }
                                                   >> 456     }
                                                   >> 457 
                                                   >> 458   isNormalized = true;
                                                   >> 459 
                                                   >> 460 
                                                   >> 461   /*
                                                   >> 462   //TESTING
                                                   >> 463   for (size_t shellID=0;shellID<numberOfShells;shellID++)
                                                   >> 464     {
                                                   >> 465       G4cout << "SHELL " << shellID << G4endl;
                                                   >> 466       G4PhysicsFreeVector* theVec =
                                                   >> 467   (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
                                                   >> 468       for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
                                                   >> 469   {
                                                   >> 470     G4double logene = theVec->GetLowEdgeEnergy(i);
                                                   >> 471     G4cout << G4Exp(logene)/MeV << " " << G4Exp((*theVec)[i]) << G4endl;
452   }                                               472   }
453     }                                             473     }
454   fIsNormalized = true;                        << 474   */
455                                                   475 
456   return;                                         476   return;
457 }                                                 477 }
458                                                   478