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

#include <INukeOsetTable.h>

Inheritance diagram for INukeOsetTable:

Collaboration diagram for INukeOsetTable:

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.

References badFile(), checkIntegrity(), fDensityBinWidth, fDensityHandler, fEnergyBinWidth, fEnergyHandler, fNDensityBins, fNEnergyBins, processLine(), and INukeOsetTable::PointHandler::setHandler().

                                                       : fNDensityBins (0), fNEnergyBins (0), 
                                                           fDensityBinWidth (0.0), fEnergyBinWidth (0.0)
 {
   // open file with Oset table
   std::ifstream fileWithTables (filename);
   // check if file is open
   if (not fileWithTables.is_open()) badFile (filename, 1);
   // temporary string for getline
   std::string tempLine;
   // line counter (used in the case of error)
   int lineCounter = 0;
 
   // skip comments lines at the beginning of the file
   while (fileWithTables.get() == '#' and ++lineCounter)
     getline (fileWithTables, tempLine);
   // process other lines
   while (getline (fileWithTables, tempLine) and ++lineCounter)
     if (int error = processLine (tempLine))   // get exit code
       badFile (filename, error, lineCounter); // stop if exit code != 0
 
   // reset counter (in the case other file will be loaded)
   checkIntegrity (-1.0, -1.0);
 
   fDensityHandler.setHandler (fDensityBinWidth, fNDensityBins);
   fEnergyHandler.setHandler  (fEnergyBinWidth,  fNEnergyBins);
 }

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.

Referenced by INukeOsetTable().

 {
   std::cerr << "\nERROR: " << file << " is corrupted (";
   
   switch (errorCode)
   {
     case 1: std::cerr << "could not be opened).\n"; exit(errorCode);
     case 2: std::cerr << "wrong density"; break;
     case 3: std::cerr << "different number of energy bins per density"; break;
     case 4: std::cerr << "wrong energy"; break;
     default: std::cerr << "undefined error";
   }
 
   std::cerr << " in line " << line << ").\n";
   exit(errorCode);
 }

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.

References fDensityBinWidth, fEnergyBinWidth, fNDensityBins, fNEnergyBins, INukeOset::fNuclearDensity, INukeOset::fPionKineticEnergy, and osetUtils::isEqual().

Referenced by INukeOsetTable(), and processLine().

 {
   static unsigned int energyBinCounter = 0; // #energy bins for current density
   
   if (densityValue < 0) // checkIntegrity(-1.0, -1.0) is called to reset counter
     energyBinCounter = -1;
   else if (not isEqual (fNuclearDensity, densityValue)) // on new density value
   {
     fNDensityBins++; // update nof density bins
 
     if (fNuclearDensity >= 0) // skip first entry (density = -1.0)
     {
       // calculate current bin width
       double currentDensityBinWidth = densityValue - fNuclearDensity;
 
       // assign first density bin width
       // or check if next bins are the same as first
       // stop program with corresponding error message (exit code 2) if not ==
       if (fDensityBinWidth == 0.0) fDensityBinWidth = currentDensityBinWidth;
       else if (not isEqual (fDensityBinWidth, currentDensityBinWidth)) return 2;
       
       // assign energy bins for first density value
       // or check if next bins are the same as first
       if (fNEnergyBins == 0) fNEnergyBins = energyBinCounter;
       else if (fNEnergyBins != energyBinCounter) return 3; // exit code 3
 
       // reset counter and last energy
       energyBinCounter = 0;
       fPionKineticEnergy  = -1;
     }
   }
   else if (fPionKineticEnergy >= 0) // check energy width consistency (skip first)
   {
     // calculate current energy bin width
     double currentEnergyBinWidth = energyValue - fPionKineticEnergy;
 
     // assign first energy bin width
     // or check if next bins are the same as first
     // stop program with corresponding error message (exit code 4) if not ==
     if (fEnergyBinWidth == 0.0) fEnergyBinWidth = currentEnergyBinWidth;
     else if (not isEqual(fEnergyBinWidth, currentEnergyBinWidth)) return 4;
   }
 
   // increase energy bin counter (for current density)
   energyBinCounter++;
   // update last values
   fPionKineticEnergy  = energyValue;
   fNuclearDensity = densityValue;
 
   return 0; // no errors
 }

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.

References fDensityBinWidth, fDensityHandler, fEnergyBinWidth, fEnergyHandler, fNEnergyBins, INukeOsetTable::PointHandler::highWeight, INukeOsetTable::PointHandler::index, INukeOsetTable::PointHandler::isEdge, and INukeOsetTable::PointHandler::lowWeight.

Referenced by setCrossSections().

 {
   // take four points adjacent to (density, energy) = (d,E):
   // (d0, E0), (d1, E0), (d0, E1), (d1, E1)
   // where d0 < d < d1, E0 < E < E1
   // each point goes in with weight = proper distance
 
   // total bin width for normalization
   static const double totalBinWidth = fDensityBinWidth * fEnergyBinWidth;
   
   // extract low boundary value from table
   double lowValue = data[fEnergyHandler.index +
                          fDensityHandler.index * fNEnergyBins]; // (d0, E0)
   // high boundaries = low boundary if on edge
   double highDensityValue = lowValue; // (d1, E0)
   double highEnergyValue  = lowValue; // (d0, E1)
   double highValue        = lowValue; // (d1, E1)
 
   // extract high boundaries values from table (if not on edge)
   if (not fDensityHandler.isEdge)
   {
     highDensityValue = data[fEnergyHandler.index +
                             (fDensityHandler.index + 1) * fNEnergyBins];
 
     if (not fEnergyHandler.isEdge)
     {
       highEnergyValue = data[fEnergyHandler.index + 1 +
                              fDensityHandler.index * fNEnergyBins];
       highValue       = data[fEnergyHandler.index + 1 +
                              (fDensityHandler.index + 1) * fNEnergyBins];
     }
   }
   else if (not fEnergyHandler.isEdge)
     highEnergyValue = data[fEnergyHandler.index + 1 +
                            fDensityHandler.index * fNEnergyBins];
 
   lowValue         *= fDensityHandler.lowWeight  * fEnergyHandler.lowWeight;
   highDensityValue *= fDensityHandler.highWeight * fEnergyHandler.lowWeight;
   highEnergyValue  *= fDensityHandler.lowWeight  * fEnergyHandler.highWeight;
   highValue        *= fDensityHandler.highWeight * fEnergyHandler.highWeight;
   
   return (lowValue + highDensityValue + highEnergyValue + highValue) /
           totalBinWidth;
 }

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.

References checkIntegrity(), fAbsorptionCrossSectionTable, fCexCrossSectionTable, INukeOset::fNChannels, and fQelCrossSectionTable.

Referenced by INukeOsetTable().

 {
   std::istringstream splitLine (line); // use istringstream to split string
   // line = d; KE; pi+n: tot fabs fcex; pi+p: tot fabs fcex; pi0: tot fabs fcex
 
   double currentDensityValue, currentEnergyValue; // will be extracted from line
 
   // get current density and energy value
   splitLine >> currentDensityValue >> currentEnergyValue;
 
   if (int exitCode = checkIntegrity (currentDensityValue, currentEnergyValue))
     return exitCode; // stop in the case of error (exit code != 0)
 
   for (unsigned int i = 0; i < fNChannels; i++) // channel loop
   {
     // get qel and cex cross sections for i-th channel
     double xsecQel, xsecCex;
     splitLine >> xsecQel >> xsecCex;
     // save them in proper vector
     fQelCrossSectionTable[i].push_back (xsecQel);
     fCexCrossSectionTable[i].push_back (xsecCex);
   } // channel loop
 
   // get absorption cross section
   double absorption;
   splitLine >> absorption;
   // save it in proper vector
   fAbsorptionCrossSectionTable.push_back (absorption);
 
   return 0; // no errors
 }

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.

References INukeOset::fAbsorptionCrossSection, fAbsorptionCrossSectionTable, INukeOset::fCexCrossSections, fCexCrossSectionTable, INukeOset::fNChannels, INukeOset::fQelCrossSections, fQelCrossSectionTable, and interpolate().

Referenced by setupOset().

 {
     for (unsigned int i = 0; i < fNChannels; i++) // channel loop
     {
       fQelCrossSections[i] = interpolate (fQelCrossSectionTable[i]);
       fCexCrossSections[i] = interpolate (fCexCrossSectionTable[i]);
     }
 
     fAbsorptionCrossSection = interpolate (fAbsorptionCrossSectionTable); 
 }

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.

References fDensityHandler, fEnergyHandler, INukeOset::fNuclearDensity, INukeOset::fPionKineticEnergy, setCrossSections(), INukeOset::setCrossSections(), and INukeOsetTable::PointHandler::update().

 {
   fNuclearDensity    = density;
   fPionKineticEnergy = pionTk;
   fDensityHandler.update (fNuclearDensity);    // update density index / weights
   fEnergyHandler.update  (fPionKineticEnergy); // update energy index / weights
   setCrossSections();
   INukeOset::setCrossSections (pionPDG, protonFraction);  
 }