Classes | Public Member Functions | Private Member Functions | Private Attributes | List of all members
INukeOsetTable Class Reference

Table-based implementation of Oset model. More...

#include <INukeOsetTable.h>

Inheritance diagram for INukeOsetTable:
INukeOset

Classes

struct  PointHandler
 handle table's index and weights for given density and energy More...
 

Public Member Functions

 INukeOsetTable (const char *filename)
 constructor More...
 
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) More...
 
- Public Member Functions inherited from INukeOset
 INukeOset ()
 contructor More...
 
double getTotalCrossSection () const
 return total = (qel+cex+abs) cross section More...
 
double getCexCrossSection () const
 return cex cross section More...
 
double getAbsorptionCrossSection () const
 return absorption cross section More...
 
double getCexFraction () const
 return fraction of cex events More...
 
double getAbsorptionFraction () const
 return fraction of absorption events More...
 

Private Member Functions

double interpolate (const std::vector< double > &data) const
 interpolate cross section (method fixed for Oset tables) More...
 
int processLine (const std::string &line)
 process single line from table file, push values to proper vector (method fixed for Oset tables) More...
 
int checkIntegrity (const double &densityValue, const double &energyValue)
 check if data in file is consistent (method fixed for Oset tables) More...
 
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) More...
 
void setCrossSections ()
 calculalte cross sections for each channel More...
 

Private Attributes

std::vector< double > fQelCrossSectionTable [fNChannels]
 quasi-elastic piN cross section More...
 
std::vector< double > fCexCrossSectionTable [fNChannels]
 charge-exchange piN cross section More...
 
std::vector< double > fAbsorptionCrossSectionTable
 pi absorption cross section More...
 
unsigned int fNDensityBins
 number of denisty bins More...
 
unsigned int fNEnergyBins
 number of energy bins More...
 
double fDensityBinWidth
 density step (must be fixed) More...
 
double fEnergyBinWidth
 energy step (must be fixed) More...
 
PointHandler fDensityHandler
 nuclear density handler More...
 
PointHandler fEnergyHandler
 pion kinetic energy handler More...
 

Additional Inherited Members

- Protected Member Functions inherited from INukeOset
void setCrossSections (const int &pionPDG, const double &protonFraction)
 calculate avg cross sections according to proton / neutron fraction More...
 
- Protected Attributes inherited from INukeOset
double fNuclearDensity
 nuclear density in fm-3 More...
 
double fPionKineticEnergy
 pion kinetic energy in MeV More...
 
double fTotalCrossSection
 el+cex+abs cross section (averaged over proton / neutron fraction) More...
 
double fCexCrossSection
 cex cross section (averaged over proton / neutron fraction) More...
 
double fAbsorptionCrossSection
 absorption cross section (averaged over proton / neutron fraction) More...
 
double fQelCrossSections [fNChannels]
 total qel (el+cex) cross section for each channel More...
 
double fCexCrossSections [fNChannels]
 cex cross section for each channel More...
 
- Static Protected Attributes inherited from INukeOset
static const unsigned int fNChannels = 3
 number of possible channels: pi+n, pi+p, pi0 More...
 

Detailed Description

Table-based implementation of Oset model.

Author
Tomasz Golan
Date
2015
Warning
Applicable for pion with Tk < 350 MeV
Remarks
Tables taken from NuWro's implementation of: E. Oset et al., Nucl. Phys. A484 (1988) 557-592

This implementation is kept from historical reasons. Due to different approach to cascade in GENIE and NuWro there are some normalization issues. Default Oset model can be found in INukeOsetFormula

Definition at line 25 of file INukeOsetTable.h.

Constructor & Destructor Documentation

INukeOsetTable::INukeOsetTable ( const char *  filename)

constructor

load tables from file and check its integrity

Definition at line 10 of file INukeOsetTable.cxx.

10  : fNDensityBins (0), fNEnergyBins (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 }
double fEnergyBinWidth
energy step (must be fixed)
std::string string
Definition: nybbler.cc:12
PointHandler fDensityHandler
nuclear density handler
error
Definition: include.cc:26
string filename
Definition: train.py:213
PointHandler fEnergyHandler
pion kinetic energy handler
void setHandler(const double &width, const int &bins)
constructor
unsigned int fNEnergyBins
number of energy bins
double fDensityBinWidth
density step (must be fixed)
int processLine(const std::string &line)
process single line from table file, push values to proper vector (method fixed for Oset tables) ...
int checkIntegrity(const double &densityValue, const double &energyValue)
check if data in file is consistent (method fixed for Oset tables)
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) ...
unsigned int fNDensityBins
number of denisty bins

Member Function Documentation

void INukeOsetTable::badFile ( const char *  file,
const int &  errorCode,
const int &  line = 0 
) const
private

stop program and through an error if input file is corrupted (method fixed for Oset tables)

possible errors:

  • the file with tables could not be loaded
  • density bin width is not constant
  • there is different number of energy point per density value
  • energy bin width is not constant

Definition at line 141 of file INukeOsetTable.cxx.

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 }
void line(double t, double *p, double &x, double &y, double &z)
int INukeOsetTable::checkIntegrity ( const double &  densityValue,
const double &  energyValue 
)
private

check if data in file is consistent (method fixed for Oset tables)

return 0 if all bins widths are constistent or return error code

Definition at line 80 of file INukeOsetTable.cxx.

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 }
bool isEqual(float x1, float x2)
double fEnergyBinWidth
energy step (must be fixed)
double fPionKineticEnergy
pion kinetic energy in MeV
Definition: INukeOset.h:67
unsigned int fNEnergyBins
number of energy bins
double fDensityBinWidth
density step (must be fixed)
double fNuclearDensity
nuclear density in fm-3
Definition: INukeOset.h:66
unsigned int fNDensityBins
number of denisty bins
double INukeOsetTable::interpolate ( const std::vector< double > &  data) const
private

interpolate cross section (method fixed for Oset tables)

make bilinear interpolation between four points around (density, energy)

Definition at line 159 of file INukeOsetTable.cxx.

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 }
double lowWeight
distance from high boundary
double fEnergyBinWidth
energy step (must be fixed)
PointHandler fDensityHandler
nuclear density handler
double highWeight
distance from low boundary
int index
point index = index of low boundary
bool isEdge
true if value is on edge of table (should never happen)
PointHandler fEnergyHandler
pion kinetic energy handler
unsigned int fNEnergyBins
number of energy bins
double fDensityBinWidth
density step (must be fixed)
int INukeOsetTable::processLine ( const std::string line)
private

process single line from table file, push values to proper vector (method fixed for Oset tables)

  • split line read from the file; the following structure is used:

    d; Tk; pi+n: tot fabs fcex; pi+p: tot fabs fcex; pi0: tot fabs fcex

  • check its integrity
  • save values into proper variables

Definition at line 47 of file INukeOsetTable.cxx.

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 }
static const unsigned int fNChannels
number of possible channels: pi+n, pi+p, pi0
Definition: INukeOset.h:81
std::vector< double > fCexCrossSectionTable[fNChannels]
charge-exchange piN cross section
void line(double t, double *p, double &x, double &y, double &z)
int checkIntegrity(const double &densityValue, const double &energyValue)
check if data in file is consistent (method fixed for Oset tables)
std::vector< double > fQelCrossSectionTable[fNChannels]
quasi-elastic piN cross section
std::vector< double > fAbsorptionCrossSectionTable
pi absorption cross section
void INukeOsetTable::setCrossSections ( )
privatevirtual

calculalte cross sections for each channel

assign cross sections values to proper variables using bilinear interpolation

Implements INukeOset.

Definition at line 243 of file INukeOsetTable.cxx.

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 fCexCrossSections[fNChannels]
cex cross section for each channel
Definition: INukeOset.h:85
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
std::vector< double > fCexCrossSectionTable[fNChannels]
charge-exchange piN cross section
std::vector< double > fQelCrossSectionTable[fNChannels]
quasi-elastic piN cross section
std::vector< double > fAbsorptionCrossSectionTable
pi absorption cross section
void INukeOsetTable::setupOset ( const double &  density,
const double &  pionTk,
const int &  pionPDG,
const double &  protonFraction 
)
virtual

use to set up Oset class (assign pion Tk, nuclear density etc)

  • set up density, pion Tk and their handlers (for interpolation)
  • get interpolated cross sections
  • set up proper cross section variables

Implements INukeOset.

Definition at line 229 of file INukeOsetTable.cxx.

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 }
void update(const double &newValue)
update point if changed
PointHandler fDensityHandler
nuclear density handler
double fPionKineticEnergy
pion kinetic energy in MeV
Definition: INukeOset.h:67
PointHandler fEnergyHandler
pion kinetic energy handler
virtual void setCrossSections()=0
calculalte cross sections for each channel
void setCrossSections()
calculalte cross sections for each channel
double fNuclearDensity
nuclear density in fm-3
Definition: INukeOset.h:66

Member Data Documentation

std::vector<double> INukeOsetTable::fAbsorptionCrossSectionTable
private

pi absorption cross section

vector contains values in the following order: d0 e0, d0 e1, ... , d0 en, d1 e0 ...

Definition at line 55 of file INukeOsetTable.h.

std::vector<double> INukeOsetTable::fCexCrossSectionTable[fNChannels]
private

charge-exchange piN cross section

vector contains values in the following order: d0 e0, d0 e1, ... , d0 en, d1 e0 ...
channel = 0 -> pi+n or pi-p, 1 -> pi+p or pi-n, 2 -> pi0

Definition at line 49 of file INukeOsetTable.h.

double INukeOsetTable::fDensityBinWidth
private

density step (must be fixed)

Definition at line 59 of file INukeOsetTable.h.

PointHandler INukeOsetTable::fDensityHandler
private

nuclear density handler

Definition at line 100 of file INukeOsetTable.h.

double INukeOsetTable::fEnergyBinWidth
private

energy step (must be fixed)

Definition at line 60 of file INukeOsetTable.h.

PointHandler INukeOsetTable::fEnergyHandler
private

pion kinetic energy handler

Definition at line 101 of file INukeOsetTable.h.

unsigned int INukeOsetTable::fNDensityBins
private

number of denisty bins

Definition at line 57 of file INukeOsetTable.h.

unsigned int INukeOsetTable::fNEnergyBins
private

number of energy bins

Definition at line 58 of file INukeOsetTable.h.

std::vector<double> INukeOsetTable::fQelCrossSectionTable[fNChannels]
private

quasi-elastic piN cross section

vector contains values in the following order: d0 e0, d0 e1, ... , d0 en, d1 e0 ...
channel = 0 -> pi+n or pi-p, 1 -> pi+p or pi-n, 2 -> pi0

Definition at line 42 of file INukeOsetTable.h.


The documentation for this class was generated from the following files: