GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
INukeOsetTable.cxx
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <sstream>
4 #include <algorithm>
5 #include "INukeOsetTable.h"
6 
7 using namespace osetUtils;
8 
9 //! load tables from file and check its integrity
10 INukeOsetTable :: INukeOsetTable (const char *filename) : fNDensityBins (0), fNEnergyBins (0),
11  fDensityBinWidth (0.0), fEnergyBinWidth (0.0)
12 {
13  // open file with Oset table
14  std::ifstream fileWithTables (filename);
15  // check if file is open
16  if (not fileWithTables.is_open()) badFile (filename, 1);
17  // temporary string for getline
18  std::string tempLine;
19  // line counter (used in the case of error)
20  int lineCounter = 0;
21 
22  // skip comments lines at the beginning of the file
23  while (fileWithTables.get() == '#' and ++lineCounter)
24  getline (fileWithTables, tempLine);
25  // process other lines
26  while (getline (fileWithTables, tempLine) and ++lineCounter)
27  if (int error = processLine (tempLine)) // get exit code
28  badFile (filename, error, lineCounter); // stop if exit code != 0
29 
30  // reset counter (in the case other file will be loaded)
31  checkIntegrity (-1.0, -1.0);
32 
35 }
36 
37 // process single line from table file, push values to proper vector
38 /*! <ul>
39  * <li> split line read from the file; the following structure is used:
40  * \n \n
41  * d; Tk; pi+n: tot fabs fcex; pi+p: tot fabs fcex; pi0: tot fabs fcex
42  * \n \n
43  * <li> check its integrity
44  * <li> save values into proper variables
45  * </ul>
46  */
47 int INukeOsetTable :: processLine (const std::string &line)
48 {
49  std::istringstream splitLine (line); // use istringstream to split string
50  // line = d; KE; pi+n: tot fabs fcex; pi+p: tot fabs fcex; pi0: tot fabs fcex
51 
52  double currentDensityValue, currentEnergyValue; // will be extracted from line
53 
54  // get current density and energy value
55  splitLine >> currentDensityValue >> currentEnergyValue;
56 
57  if (int exitCode = checkIntegrity (currentDensityValue, currentEnergyValue))
58  return exitCode; // stop in the case of error (exit code != 0)
59 
60  for (unsigned int i = 0; i < fNChannels; i++) // channel loop
61  {
62  // get qel and cex cross sections for i-th channel
63  double xsecQel, xsecCex;
64  splitLine >> xsecQel >> xsecCex;
65  // save them in proper vector
66  fQelCrossSectionTable[i].push_back (xsecQel);
67  fCexCrossSectionTable[i].push_back (xsecCex);
68  } // channel loop
69 
70  // get absorption cross section
71  double absorption;
72  splitLine >> absorption;
73  // save it in proper vector
74  fAbsorptionCrossSectionTable.push_back (absorption);
75 
76  return 0; // no errors
77 }
78 
79 //! return 0 if all bins widths are constistent or return error code
80 int INukeOsetTable :: checkIntegrity (const double &densityValue, const double &energyValue)
81 {
82  static unsigned int energyBinCounter = 0; // #energy bins for current density
83 
84  if (densityValue < 0) // checkIntegrity(-1.0, -1.0) is called to reset counter
85  energyBinCounter = -1;
86  else if (not isEqual (fNuclearDensity, densityValue)) // on new density value
87  {
88  fNDensityBins++; // update nof density bins
89 
90  if (fNuclearDensity >= 0) // skip first entry (density = -1.0)
91  {
92  // calculate current bin width
93  double currentDensityBinWidth = densityValue - fNuclearDensity;
94 
95  // assign first density bin width
96  // or check if next bins are the same as first
97  // stop program with corresponding error message (exit code 2) if not ==
98  if (fDensityBinWidth == 0.0) fDensityBinWidth = currentDensityBinWidth;
99  else if (not isEqual (fDensityBinWidth, currentDensityBinWidth)) return 2;
100 
101  // assign energy bins for first density value
102  // or check if next bins are the same as first
103  if (fNEnergyBins == 0) fNEnergyBins = energyBinCounter;
104  else if (fNEnergyBins != energyBinCounter) return 3; // exit code 3
105 
106  // reset counter and last energy
107  energyBinCounter = 0;
108  fPionKineticEnergy = -1;
109  }
110  }
111  else if (fPionKineticEnergy >= 0) // check energy width consistency (skip first)
112  {
113  // calculate current energy bin width
114  double currentEnergyBinWidth = energyValue - fPionKineticEnergy;
115 
116  // assign first energy bin width
117  // or check if next bins are the same as first
118  // stop program with corresponding error message (exit code 4) if not ==
119  if (fEnergyBinWidth == 0.0) fEnergyBinWidth = currentEnergyBinWidth;
120  else if (not isEqual(fEnergyBinWidth, currentEnergyBinWidth)) return 4;
121  }
122 
123  // increase energy bin counter (for current density)
124  energyBinCounter++;
125  // update last values
126  fPionKineticEnergy = energyValue;
127  fNuclearDensity = densityValue;
128 
129  return 0; // no errors
130 }
131 
132 /*! possible errors:
133  *
134  * <ul>
135  * <li> the file with tables could not be loaded
136  * <li> density bin width is not constant
137  * <li> there is different number of energy point per density value
138  * <li> energy bin width is not constant
139  * </ul>
140  */
141 void INukeOsetTable :: badFile (const char* file, const int &errorCode, const int &line) const
142 {
143  std::cerr << "\nERROR: " << file << " is corrupted (";
144 
145  switch (errorCode)
146  {
147  case 1: std::cerr << "could not be opened).\n"; exit(errorCode);
148  case 2: std::cerr << "wrong density"; break;
149  case 3: std::cerr << "different number of energy bins per density"; break;
150  case 4: std::cerr << "wrong energy"; break;
151  default: std::cerr << "undefined error";
152  }
153 
154  std::cerr << " in line " << line << ").\n";
155  exit(errorCode);
156 }
157 
158 //! make bilinear interpolation between four points around (density, energy)
159 double INukeOsetTable :: interpolate (const std::vector<double> &data) const
160 {
161  // take four points adjacent to (density, energy) = (d,E):
162  // (d0, E0), (d1, E0), (d0, E1), (d1, E1)
163  // where d0 < d < d1, E0 < E < E1
164  // each point goes in with weight = proper distance
165 
166  // total bin width for normalization
167  static const double totalBinWidth = fDensityBinWidth * fEnergyBinWidth;
168 
169  // extract low boundary value from table
170  double lowValue = data[fEnergyHandler.index +
171  fDensityHandler.index * fNEnergyBins]; // (d0, E0)
172  // high boundaries = low boundary if on edge
173  double highDensityValue = lowValue; // (d1, E0)
174  double highEnergyValue = lowValue; // (d0, E1)
175  double highValue = lowValue; // (d1, E1)
176 
177  // extract high boundaries values from table (if not on edge)
178  if (not fDensityHandler.isEdge)
179  {
180  highDensityValue = data[fEnergyHandler.index +
182 
183  if (not fEnergyHandler.isEdge)
184  {
185  highEnergyValue = data[fEnergyHandler.index + 1 +
187  highValue = data[fEnergyHandler.index + 1 +
189  }
190  }
191  else if (not fEnergyHandler.isEdge)
192  highEnergyValue = data[fEnergyHandler.index + 1 +
194 
196  highDensityValue *= fDensityHandler.highWeight * fEnergyHandler.lowWeight;
197  highEnergyValue *= fDensityHandler.lowWeight * fEnergyHandler.highWeight;
199 
200  return (lowValue + highDensityValue + highEnergyValue + highValue) /
201  totalBinWidth;
202 }
203 
204 //! set up table index and weights for given point
205 void INukeOsetTable :: PointHandler :: update (const double &newValue)
206 {
207  value = newValue; // update value
208  index = value / binWidth; // update index
209 
210  // in the case value > max value use max; check if it is on edge
211  if (index >= nBins - 1)
212  {
213  index = nBins - 1;
214  isEdge = true;
215  }
216  else isEdge = false;
217  // calcualte weights for boundary values
218  lowWeight = (index + 1) * binWidth - value;
220 }
221 
222 
223 /*! <ul>
224  * <li> set up density, pion Tk and their handlers (for interpolation)
225  * <li> get interpolated cross sections
226  * <li> set up proper cross section variables
227  * </ul>
228  */
229 void INukeOsetTable :: setupOset (const double &density, const double &pionTk, const int &pionPDG,
230  const double &protonFraction)
231 {
232  fNuclearDensity = density;
233  fPionKineticEnergy = pionTk;
234  fDensityHandler.update (fNuclearDensity); // update density index / weights
235  fEnergyHandler.update (fPionKineticEnergy); // update energy index / weights
237  INukeOset::setCrossSections (pionPDG, protonFraction);
238 }
239 
240 /*! assign cross sections values to proper variables
241  * using bilinear interpolation
242  */
244 {
245  for (unsigned int i = 0; i < fNChannels; i++) // channel loop
246  {
249  }
250 
252 }
double interpolate(const std::vector< double > &data) const
interpolate cross section (method fixed for Oset tables)
double lowWeight
distance from high boundary
void update(const double &newValue)
update point if changed
double fEnergyBinWidth
energy step (must be fixed)
bool isEqual(const double &x, const double &y, const double &epsilon)
Definition: INukeOset.h:100
PointHandler fDensityHandler
nuclear density handler
double fPionKineticEnergy
pion kinetic energy in MeV
Definition: INukeOset.h:67
double highWeight
distance from low boundary
int index
point index = index of low boundary
double fCexCrossSections[fNChannels]
cex cross section for each channel
Definition: INukeOset.h:85
bool isEdge
true if value is on edge of table (should never happen)
double fAbsorptionCrossSection
absorption cross section (averaged over proton / neutron fraction)
Definition: INukeOset.h:74
double fQelCrossSections[fNChannels]
total qel (el+cex) cross section for each channel
Definition: INukeOset.h:84
static const unsigned int fNChannels
number of possible channels: pi+n, pi+p, pi0
Definition: INukeOset.h:81
PointHandler fEnergyHandler
pion kinetic energy handler
std::vector< double > fCexCrossSectionTable[fNChannels]
charge-exchange piN cross section
double binWidth
bin width used to calculate distances
void setHandler(const double &width, const int &bins)
constructor
unsigned int fNEnergyBins
number of energy bins
double fDensityBinWidth
density step (must be fixed)
virtual void setCrossSections()=0
calculalte cross sections for each channel
void setupOset(const double &density, const double &pionTk, const int &pionPDG, const double &protonFraction)
use to set up Oset class (assign pion Tk, nuclear density etc)
int processLine(const std::string &line)
process single line from table file, push values to proper vector (method fixed for Oset tables) ...
INukeOsetTable(const char *filename)
constructor
double value
exact value as read from table
void setCrossSections()
calculalte cross sections for each channel
int checkIntegrity(const double &densityValue, const double &energyValue)
check if data in file is consistent (method fixed for Oset tables)
double fNuclearDensity
nuclear density in fm-3
Definition: INukeOset.h:66
void badFile(const char *file, const int &errorCode, const int &line=0) const
stop program and through an error if input file is corrupted (method fixed for Oset tables) ...
int nBins
nBins to check isEdge
std::vector< double > fQelCrossSectionTable[fNChannels]
quasi-elastic piN cross section
std::vector< double > fAbsorptionCrossSectionTable
pi absorption cross section
unsigned int fNDensityBins
number of denisty bins