genie::Spline Class Reference

A numeric analysis tool class for interpolating 1-D functions. More...

#include <Spline.h>

Inheritance diagram for genie::Spline:

Collaboration diagram for genie::Spline:

Public Member Functions
	Spline ()
	Spline (string filename, string xtag="", string ytag="", bool is_xml=false)
	Spline (TNtupleD *ntuple, string xy, string cut="")
	Spline (TTree *tree, string xy, string cut="")
	Spline (TSQLServer *db, string query)
	Spline (int nentries, double x[], double y[])
	Spline (int nentries, float x[], float y[])
	Spline (const Spline &spline)
	Spline (const TSpline3 &spline, int nknots)
virtual	~Spline ()
bool	LoadFromXmlFile (string filename, string xtag, string ytag)
bool	LoadFromAsciiFile (string filename)
bool	LoadFromNtuple (TNtupleD *nt, string xy, string cut="")
bool	LoadFromTree (TTree *tr, string xy, string cut="")
bool	LoadFromDBase (TSQLServer *db, string query)
bool	LoadFromTSpline3 (const TSpline3 &spline, int nknots)
int	NKnots (void) const
void	GetKnot (int iknot, double &x, double &y) const
double	GetKnotX (int iknot) const
double	GetKnotY (int iknot) const
double	XMin (void) const
double	XMax (void) const
double	YMax (void) const
double	Evaluate (double x) const
bool	IsWithinValidRange (double x) const
void	SetName (string name)
string	Name (void) const
void	YCanBeNegative (bool tf)
void	SaveAsXml (string filename, string xtag, string ytag, string name="") const
void	SaveAsXml (ofstream &str, string xtag, string ytag, string name="") const
void	SaveAsText (string filename, string format="%10.6f\t%10.6f") const
void	SaveAsROOT (string filename, string name="", bool recreate=false) const
TGraph *	GetAsTGraph (int np=500, bool xscaling=false, bool inlog=false, double fx=1., double fy=1.) const
TSpline3 *	GetAsTSpline (void) const
void	FindClosestKnot (double x, double &xknot, double &yknot, Option_t *opt="-+") const
bool	ClosestKnotValueIsZero (double x, Option_t *opt="-+") const
void	Add (const Spline &spl, double c=1)
void	Multiply (const Spline &spl, double c=1)
void	Divide (const Spline &spl, double c=1)
void	Add (double a)
void	Multiply (double a)
void	Divide (double a)
void	SetType (string type)
string	GetType (void)
void	Print (ostream &stream) const

Private Member Functions
void	InitSpline (void)
void	ResetSpline (void)
void	BuildSpline (int nentries, double x[], double y[])

Private Attributes
string	fName
int	fNKnots
double	fXMin
double	fXMax
double	fYMax
TSpline3 *	fInterpolator
bool	fYCanBeNegative
TSpline5 *	fInterpolator5
ROOT::Math::Interpolator *	fGSLInterpolator
string	fInterpolatorType

Friends
ostream &	operator<< (ostream &stream, const Spline &spl)

Detailed Description

A numeric analysis tool class for interpolating 1-D functions.

      Uses ROOT's TSpline3 for the actual interpolation and can retrieve
      function (x,y(x)) pairs from an XML file, a flat ascii file, a
      TNtuple, a TTree or an SQL database.\n

      Update May 15, 2022 IK: 
      Adding as extra interpolators TSpline5 and 
      ROOT::Math::GSLInterpolator (LINEAR, POLYNOMIAL, CSPLINE, CSPLINE_PERIODIC,
      AKIMA, AKIMA_PERIODIC)

References:

[1] GENIE docdb 297

Author

Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool
Igor Kakorin kakor.nosp@m.in@j.nosp@m.inr.r.nosp@m.u Joint Institute for Nuclear Research

Created:

May 04, 2004

License:

Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 58 of file Spline.h.

Constructor & Destructor Documentation

Spline::Spline

(

)

Definition at line 52 of file Spline.cxx.

References InitSpline().

{
  this->InitSpline();
}

Spline::Spline	(	string	filename,
		string	xtag = "",
		string	ytag = "",
		bool	is_xml = false
	)

Definition at line 57 of file Spline.cxx.

References InitSpline(), LoadFromAsciiFile(), LoadFromXmlFile(), LOG, and pDEBUG.

                                                                     :
TObject()
{
  string fmt = (is_xml) ? "XML" : "ASCII";

  LOG("Spline", pDEBUG)
      << "Constructing spline from data in " << fmt << " file: " << filename;

  this->InitSpline();

  if(is_xml)
     this->LoadFromXmlFile(filename, xtag, ytag);
  else
     this->LoadFromAsciiFile(filename);
}

Spline::Spline	(	TNtupleD *	ntuple,
		string	xy,
		string	cut = ""
	)

Definition at line 73 of file Spline.cxx.

References InitSpline(), LoadFromNtuple(), LOG, and pDEBUG.

                                                        :
TObject()
{
  LOG("Spline", pDEBUG) << "Constructing spline from data in a TNtuple";

  this->InitSpline();
  this->LoadFromNtuple(ntuple, var, cut);
}

Spline::Spline	(	TTree *	tree,
		string	xy,
		string	cut = ""
	)

Definition at line 82 of file Spline.cxx.

References InitSpline(), LoadFromTree(), LOG, and pDEBUG.

                                                   :
TObject()
{
  LOG("Spline", pDEBUG) << "Constructing spline from data in a TTree";

  this->InitSpline();
  this->LoadFromTree(tree, var, cut);
}

Spline::Spline	(	TSQLServer *	db,
		string	query
	)

Definition at line 91 of file Spline.cxx.

References InitSpline(), LoadFromDBase(), LOG, and pDEBUG.

                                            :
TObject()
{
  LOG("Spline", pDEBUG) << "Constructing spline from data in a MySQL server";

  this->InitSpline();
  this->LoadFromDBase(db, query);
}

Spline::Spline	(	int	nentries,
		double	x[],
		double	y[]
	)

Definition at line 100 of file Spline.cxx.

References BuildSpline(), InitSpline(), LOG, and pDEBUG.

                                                   :
TObject()
{
  LOG("Spline", pDEBUG)
                 << "Constructing spline from the arrays passed to the ctor";

  this->InitSpline();
  this->BuildSpline(nentries, x, y);
}

Spline::Spline	(	int	nentries,
		float	x[],
		float	y[]
	)

Definition at line 110 of file Spline.cxx.

References BuildSpline(), InitSpline(), LOG, and pDEBUG.

                                                 :
TObject()
{
  LOG("Spline", pDEBUG)
                 << "Constructing spline from the arrays passed to the ctor";

  double * dblx = new double[nentries];
  double * dbly = new double[nentries];

  for(int i = 0; i < nentries; i++) {
     dblx[i] = double( x[i] );
     dbly[i] = double( y[i] );
  }

  this->InitSpline();
  this->BuildSpline(nentries, dblx, dbly);

  delete [] dblx;
  delete [] dbly;
}

Spline::Spline

(

const Spline &

spline

)

Definition at line 131 of file Spline.cxx.

References GetAsTSpline(), InitSpline(), LoadFromTSpline3(), LOG, NKnots(), and pDEBUG.

                                    :
  TObject(), fInterpolator(0)
{
  LOG("Spline", pDEBUG) << "Spline copy constructor";
  this->InitSpline();
  this->LoadFromTSpline3( *spline.GetAsTSpline(), spline.NKnots() );
}

Spline::Spline	(	const TSpline3 &	spline,
		int	nknots
	)

Definition at line 139 of file Spline.cxx.

References InitSpline(), LoadFromTSpline3(), LOG, and pDEBUG.

                                                  :
  TObject(), fInterpolator(0)
{
  LOG("Spline", pDEBUG)
                    << "Constructing spline from the input TSpline3 object";
  this->InitSpline();
  this->LoadFromTSpline3( spline, nknots );
}

Spline::~Spline ( )

virtual

Definition at line 148 of file Spline.cxx.

References fGSLInterpolator, fInterpolator, and fInterpolator5.

{
  if(fInterpolator) delete fInterpolator;
  if(fInterpolator5) delete fInterpolator5;
  if(fGSLInterpolator) delete fGSLInterpolator;
}

Member Function Documentation

void Spline::Add	(	const Spline &	spl,
		double	c = 1
	)

Definition at line 587 of file Spline.cxx.

References BuildSpline(), Evaluate(), GetKnot(), IsWithinValidRange(), LOG, NKnots(), pERROR, ResetSpline(), XMax(), and XMin().

{
  // continue only if the input spline is defined at a wider x-range
  double xmin = this->XMin();
  double xmax = this->XMax();
  bool ok = spl.IsWithinValidRange(xmin) && spl.IsWithinValidRange(xmax);
  if(!ok) {
    LOG("Spline", pERROR) << "** Can't add splines. X-range mismatch!";
    return;
  }

  int nknots = this->NKnots();
  double * x = new double[nknots];
  double * y = new double[nknots];

  for(int i=0; i<nknots; i++) {
    this->GetKnot(i,x[i],y[i]);
    y[i] += (c * spl.Evaluate(x[i]));
  }
  this->ResetSpline();
  this->BuildSpline(nknots,x,y);
  delete [] x;
  delete [] y;
}

void Spline::Add

(

double

a

)

Definition at line 670 of file Spline.cxx.

References a, BuildSpline(), GetKnot(), NKnots(), and ResetSpline().

{
  int nknots = this->NKnots();
  double * x = new double[nknots];
  double * y = new double[nknots];

  for(int i=0; i<nknots; i++) {
    this->GetKnot(i,x[i],y[i]);
    y[i]+=a;
  }
  this->ResetSpline();
  this->BuildSpline(nknots,x,y);
  delete [] x;
  delete [] y;
}

void Spline::BuildSpline ( int  nentries, double  x[], double  y[]  )

private

Definition at line 750 of file Spline.cxx.

References fInterpolator, fNKnots, fXMax, fXMin, fYMax, LOG, and pDEBUG.

Referenced by Add(), Divide(), LoadFromTree(), LoadFromTSpline3(), LoadFromXmlFile(), Multiply(), and Spline().

{
  LOG("Spline", pDEBUG) << "Building spline...";

  double xmin = x[0];          // minimum x in spline
  double xmax = x[nentries-1]; // maximum x in spline

  fNKnots = nentries;
  fXMin   = xmin;
  fXMax   = xmax;
  fYMax   = y[ TMath::LocMax(nentries, y) ]; // maximum y in spline

  if(fInterpolator) delete fInterpolator;

  fInterpolator = new TSpline3("spl3", x, y, nentries, "0");
    
  LOG("Spline", pDEBUG) << "...done building spline";
}

bool Spline::ClosestKnotValueIsZero	(	double	x,
		Option_t *	opt = "-+"
	)		const

Definition at line 561 of file Spline.cxx.

References genie::utils::math::AreEqual(), and FindClosestKnot().

Referenced by Evaluate(), and genie::PhysInteractionSelector::SelectInteraction().

{
  double xknot = 0, yknot = 0;
  this->FindClosestKnot(x, xknot, yknot, opt);
  if(utils::math::AreEqual(yknot,0)) return true;
  return false;
}

void Spline::Divide	(	const Spline &	spl,
		double	c = 1
	)

Definition at line 637 of file Spline.cxx.

References BuildSpline(), Evaluate(), GetKnot(), IsWithinValidRange(), LOG, NKnots(), pERROR, ResetSpline(), XMax(), and XMin().

{
  // continue only if the input spline is defined at a wider x-range
  double xmin = this->XMin();
  double xmax = this->XMax();
  bool ok = spl.IsWithinValidRange(xmin) && spl.IsWithinValidRange(xmax);
  if(!ok) {
    LOG("Spline", pERROR) << "** Can't divide splines. X-range mismatch!";
    return;
  }

  int nknots = this->NKnots();
  double * x = new double[nknots];
  double * y = new double[nknots];

  for(int i=0; i<nknots; i++) {
    this->GetKnot(i,x[i],y[i]);
    double denom = c * spl.Evaluate(x[i]);
    bool denom_is_zero = TMath::Abs(denom) < DBL_EPSILON;
    if(denom_is_zero) {
        LOG("Spline", pERROR) << "** Refusing to divide spline knot by 0";
        delete [] x;
        delete [] y;
        return;
    }
    y[i] /= denom;
  }
  this->ResetSpline();
  this->BuildSpline(nknots,x,y);
  delete [] x;
  delete [] y;
}

void Spline::Divide

(

double

a

)

Definition at line 702 of file Spline.cxx.

References a, BuildSpline(), GetKnot(), LOG, NKnots(), pERROR, and ResetSpline().

{
  bool a_is_zero = TMath::Abs(a) < DBL_EPSILON;
  if(a_is_zero==0) {
    LOG("Spline", pERROR) << "** Refusing to divide spline by 0";
    return;
  }
  int nknots = this->NKnots();
  double * x = new double[nknots];
  double * y = new double[nknots];

  for(int i=0; i<nknots; i++) {
    this->GetKnot(i,x[i],y[i]);
    y[i]/=a;
  }
  this->ResetSpline();
  this->BuildSpline(nknots,x,y);
  delete [] x;
  delete [] y;
}

double Spline::Evaluate

(

double

x

)

const

Definition at line 363 of file Spline.cxx.

References ClosestKnotValueIsZero(), fGSLInterpolator, FindClosestKnot(), fInterpolator, fInterpolator5, fInterpolatorType, fXMax, fXMin, fYCanBeNegative, IsWithinValidRange(), LOG, pDEBUG, and pINFO.

Referenced by Add(), BuildSpectrum(), CalcTotalXSec(), genie::GMCJDriver::ComputeInteractionProbabilities(), genie::GMCJDriver::ComputeProbScales(), genie::NuclearData::DeuteriumSuppressionFactor(), Divide(), genie::GiBUURESFormFactor::FormFactors::FFRes(), genie::KineGeneratorWithCache::FindMaxXSec(), genie::AGCharm2019::GenerateCharmHadron(), GetAsTGraph(), genie::DISXSec::Integrate(), genie::HEDISXSec::Integrate(), genie::DMDISXSec::Integrate(), genie::HELeptonXSec::Integrate(), genie::AlamSimoAtharVacasSKXSec::Integrate(), genie::SPPXSec::Integrate(), genie::ReinSehgalRESXSec::Integrate(), genie::ReinSehgalRESXSecFast::Integrate(), genie::SpectralFunc1d::LoadConfig(), Multiply(), genie::CacheBranchFx::operator()(), SaveGraphsToRootFile(), genie::PhysInteractionSelector::SelectInteraction(), genie::ReinSehgalRESPXSec::XSec(), and genie::GEVGDriver::XSecSum().

{
  LOG("Spline", pDEBUG) << "Evaluating spline at point x = " << x;
  assert(!TMath::IsNaN(x));

  double y = 0;
  if( this->IsWithinValidRange(x) ) {

    // we can interpolate within the range of spline knots - be careful with
    // strange cubic spline behaviour when close to knots with y=0
    bool is0p = this->ClosestKnotValueIsZero(x, "+");
    bool is0n = this->ClosestKnotValueIsZero(x, "-");

    if(!is0p && !is0n) {
      // both knots (on the left and right are non-zero) - just interpolate
      LOG("Spline", pDEBUG) << "Point is between non-zero knots";
      if (fInterpolatorType == "TSpline3")
        y = fInterpolator->Eval(x);
      else if (fInterpolatorType == "TSpline5")
        y = fInterpolator5->Eval(x);
      else
        y = fGSLInterpolator->Eval(x);
    } else {
      // at least one of the neighboring knots has y=0
      if(is0p && is0n) {
        // both neighboring knots have y=0
        LOG("Spline", pDEBUG) << "Point is between zero knots";
        y=0;
      } else {
        // just 1 neighboring knot has y=0 - do a linear interpolation
        LOG("Spline", pDEBUG)
          << "Point has zero" << (is0n ? " left " : " right ") << "knot";
        double xpknot=0, ypknot=0, xnknot=0, ynknot=0;
        this->FindClosestKnot(x, xnknot, ynknot, "-");
        this->FindClosestKnot(x, xpknot, ypknot, "+");
        if(is0n) y = ypknot * (x-xnknot)/(xpknot-xnknot);
        else     y = ynknot * (x-xnknot)/(xpknot-xnknot);
      }
    }

  } else {
    LOG("Spline", pDEBUG) << "x = " << x
     << " is not within spline range [" << fXMin << ", " << fXMax << "]";
  }

  if(y<0 && !fYCanBeNegative) {
    LOG("Spline", pINFO) << "Negative y (" << y << ")";
    LOG("Spline", pINFO) << "x = " << x;
    LOG("Spline", pINFO) << "spline range [" << fXMin << ", " << fXMax << "]";
  }

  LOG("Spline", pDEBUG) << "Spline(x = " << x << ") = " << y;

  return y;
}

void Spline::FindClosestKnot	(	double	x,
		double &	xknot,
		double &	yknot,
		Option_t *	opt = "-+"
	)		const

Definition at line 534 of file Spline.cxx.

References fInterpolator.

Referenced by ClosestKnotValueIsZero(), and Evaluate().

{
  string option(opt);

  bool pos = (option.find("+") != string::npos);
  bool neg = (option.find("-") != string::npos);

  if(!pos && !neg) return;

  int iknot = fInterpolator->FindX(x);

  double xp=0, yp=0, xn=0, yn=0;
  fInterpolator->GetKnot(iknot,  xn,yn);
  fInterpolator->GetKnot(iknot+1,xp,yp);

  bool p = (TMath::Abs(x-xp) < TMath::Abs(x-xn));

  if(pos&&neg) {
    if(p) { xknot = xp; yknot = yp; }
    else  { xknot = xn; yknot = yn; }
  } else {
    if(pos) { xknot = xp; yknot = yp; }
    if(neg) { xknot = xn; yknot = yn; }
  }
}

TGraph * Spline::GetAsTGraph	(	int	np = 500,
		bool	xscaling = false,
		bool	inlog = false,
		double	fx = 1.,
		double	fy = 1.
	)		const

Definition at line 496 of file Spline.cxx.

References Evaluate(), fXMax, and fXMin.

Referenced by SaveToPsFile().

{
  double xmin = fXMin;
  double xmax = fXMax;

  np = TMath::Max(np,2);

  bool use_log = in_log && (fXMin>0) && (fXMax>0);

  if(use_log) {
    xmin = TMath::Log10(xmin);
    xmax = TMath::Log10(xmax);
  }

  double dx = (xmax - xmin) / (np-1);

  double * x = new double[np];
  double * y = new double[np];

  for(int i=0; i<np; i++) {
      x[i] = ( (use_log) ? TMath::Power(10, xmin+i*dx) : xmin + i*dx );
      y[i] = this->Evaluate( x[i] );

      // scale with x if needed
      if (scale_with_x) y[i] /= x[i];

      // apply x,y scaling
      y[i] *= fy;
      x[i] *= fx;
  }

  TGraph * graph = new TGraph(np, x, y);
  delete[] x;
  delete[] y;
  return graph;
}

TSpline3* genie::Spline::GetAsTSpline ( void  ) const

inline

Definition at line 108 of file Spline.h.

References fInterpolator.

Referenced by Spline().

{ return fInterpolator; }

void Spline::GetKnot	(	int	iknot,
		double &	x,
		double &	y
	)		const

Definition at line 326 of file Spline.cxx.

References fInterpolator, LOG, and pWARN.

Referenced by Add(), Divide(), Multiply(), and Print().

{
  if(!fInterpolator) {
     LOG("Spline", pWARN) << "Spline has not been built yet!";
     return;
  }
  fInterpolator->GetKnot(iknot,x,y);
}

double Spline::GetKnotX

(

int

iknot

)

const

Definition at line 335 of file Spline.cxx.

References fInterpolator, LOG, and pWARN.

{
  if(!fInterpolator) {
     LOG("Spline", pWARN) << "Spline has not been built yet!";
     return 0;
  }
  double x,y;
  fInterpolator->GetKnot(iknot,x,y);
  return x;
}

double Spline::GetKnotY

(

int

iknot

)

const

Definition at line 346 of file Spline.cxx.

References fInterpolator, LOG, and pWARN.

{
  if(!fInterpolator) {
     LOG("Spline", pWARN) << "Spline has not been built yet!";
     return 0;
  }
  double x,y;
  fInterpolator->GetKnot(iknot,x,y);
  return y;
}

string genie::Spline::GetType ( void  )

inline

Definition at line 122 of file Spline.h.

References fInterpolatorType.

{ return fInterpolatorType; }

void Spline::InitSpline ( void  )

private

Definition at line 723 of file Spline.cxx.

References fGSLInterpolator, fInterpolator, fInterpolator5, fInterpolatorType, fName, fXMax, fXMin, fYCanBeNegative, fYMax, LOG, and pDEBUG.

Referenced by ResetSpline(), and Spline().

{
  LOG("Spline", pDEBUG) << "Initializing spline...";

  fName = "genie-spline";
  fXMin = 0.0;
  fXMax = 0.0;
  fYMax = 0.0;
  
  fInterpolator = 0;
  fInterpolator5 = 0;
  fGSLInterpolator = 0;
  fInterpolatorType = "TSpline3";

  fYCanBeNegative = false;

  LOG("Spline", pDEBUG) << "...done initializing spline";
}

bool Spline::IsWithinValidRange

(

double

x

)

const

Definition at line 357 of file Spline.cxx.

References fXMax, and fXMin.

Referenced by Add(), Divide(), Evaluate(), and Multiply().

{
  bool is_in_range = (fXMin <= x && x <= fXMax);
  return is_in_range;
}

bool Spline::LoadFromAsciiFile

(

string

filename

)

Definition at line 242 of file Spline.cxx.

References LoadFromNtuple(), LOG, and pDEBUG.

Referenced by Spline().

{
  LOG("Spline", pDEBUG) << "Retrieving data from file: " << filename;

  TNtupleD nt("ntuple","","x:y");
  nt.ReadFile(filename.c_str());

  this->LoadFromNtuple(&nt, "x:y");
  return true;
}

bool Spline::LoadFromDBase	(	TSQLServer *	db,
		string	query
	)

Definition at line 299 of file Spline.cxx.

References LOG, and pDEBUG.

Referenced by Spline().

{
  LOG("Spline", pDEBUG) << "Retrieving data from data-base: ";
  return false;
}

bool Spline::LoadFromNtuple	(	TNtupleD *	nt,
		string	xy,
		string	cut = ""
	)

Definition at line 253 of file Spline.cxx.

References LoadFromTree().

Referenced by LoadFromAsciiFile(), and Spline().

{
  if(!nt) return false;

  TTree * tree = dynamic_cast<TTree *> (nt);
  return this->LoadFromTree(tree,var,cut);
}

bool Spline::LoadFromTree	(	TTree *	tr,
		string	xy,
		string	cut = ""
	)

Definition at line 261 of file Spline.cxx.

References BuildSpline(), LOG, and pDEBUG.

Referenced by LoadFromNtuple(), and Spline().

{
  LOG("Spline", pDEBUG) << "Retrieving data from tree: " << tree->GetName();

  if(!cut.size()) tree->Draw(var.c_str(), "",          "GOFF");
  else            tree->Draw(var.c_str(), cut.c_str(), "GOFF");

  TH2F * hst = (TH2F*)gROOT->FindObject("htemp");
  if(hst) { hst->SetDirectory(0); delete hst; }

  // Now, take into account that the data retrieved from the ntuple would
  // not be sorted in x and the resulting spline will be bogus...
  // Sort the x array, use the x-sorting to re-arrange the y array and only
  // then build the spline..

  int n = tree->GetSelectedRows();

  int *    idx = new int[n];
  double * x   = new double[n];
  double * y   = new double[n];

  TMath::Sort(n,tree->GetV1(),idx,false);

  for(int i=0; i<n; i++) {
     int ii = idx[i];
     x[i]   = (tree->GetV1())[ii];
     y[i]   = (tree->GetV2())[ii];
  }

  this->BuildSpline(n,x,y);

  delete [] idx;
  delete [] x;
  delete [] y;

  return true;
}

bool Spline::LoadFromTSpline3	(	const TSpline3 &	spline,
		int	nknots
	)

Definition at line 305 of file Spline.cxx.

References BuildSpline().

Referenced by Spline().

{
  int nentries = nknots;

  double * x = new double[nentries];
  double * y = new double[nentries];
  double cx = 0, cy = 0;

  for(int i = 0; i < nentries; i++) {
     spline.GetKnot(i, cx, cy);
     x[i] = cx;
     y[i] = cy;
  }
  this->BuildSpline(nentries, x, y);

  delete [] x;
  delete [] y;

  return true;
}

bool Spline::LoadFromXmlFile	(	string	filename,
		string	xtag,
		string	ytag
	)

Definition at line 155 of file Spline.cxx.

References BuildSpline(), genie::utils::xml::GetAttribute(), LOG, pDEBUG, pERROR, pINFO, genie::utils::str::TrimSpaces(), and genie::utils::xml::TrimSpaces().

Referenced by Spline().

{
  LOG("Spline", pDEBUG) << "Retrieving data from file: " << filename;

  xmlDocPtr xml_doc = xmlParseFile(filename.c_str());

  if(xml_doc==NULL) {
    LOG("Spline", pERROR)
           << "XML file could not be parsed! [filename: " << filename << "]";
    return false;
  }

  xmlNodePtr xmlCur = xmlDocGetRootElement(xml_doc);

  if(xmlCur==NULL) {
    LOG("Spline", pERROR)
         << "XML doc. has null root element! [filename: " << filename << "]";
    return false;
  }
  if( xmlStrcmp(xmlCur->name, (const xmlChar *) "spline") ) {
    LOG("Spline", pERROR)
      << "XML doc. has invalid root element! [filename: " << filename << "]";
    return false;
  }

  string name = utils::str::TrimSpaces(
                    utils::xml::GetAttribute(xmlCur, "name"));
  string snkn = utils::str::TrimSpaces(
                    utils::xml::GetAttribute(xmlCur, "nknots"));
  int nknots = atoi( snkn.c_str() );

  LOG("Spline", pINFO)
             << "Parsing XML spline: " << name << ", nknots = " << nknots;

  int      iknot = 0;
  double * vx    = new double[nknots];
  double * vy    = new double[nknots];

  xmlNodePtr xmlSplChild = xmlCur->xmlChildrenNode; // <knots>'s

  // loop over all xml tree nodes that are children of the <spline> node
  while (xmlSplChild != NULL) {
     LOG("Spline", pDEBUG)
                 << "Got <spline> children node: " << xmlSplChild->name;

     // enter everytime you find a <knot> tag
     if( (!xmlStrcmp(xmlSplChild->name, (const xmlChar *) "knot")) ) {

        xmlNodePtr xmlKnotChild = xmlSplChild->xmlChildrenNode;

        // loop over all xml tree nodes that are children of this <knot>
        while (xmlKnotChild != NULL) {
           LOG("Spline", pDEBUG)
                 << "Got <knot> children node: "  << xmlKnotChild->name;

           // enter everytime you find a <E> or a <xsec> tag
           const xmlChar * tag = xmlKnotChild->name;
           bool is_xtag = ! xmlStrcmp(tag,(const xmlChar *) xtag.c_str());
           bool is_ytag = ! xmlStrcmp(tag,(const xmlChar *) ytag.c_str());
           if (is_xtag || is_ytag) {
              xmlNodePtr xmlValTagChild = xmlKnotChild->xmlChildrenNode;
              string val = utils::xml::TrimSpaces(
                         xmlNodeListGetString(xml_doc, xmlValTagChild, 1));

              if (is_xtag) vx[iknot] = atof(val.c_str());
              if (is_ytag) vy[iknot] = atof(val.c_str());

              xmlFree(xmlValTagChild);
              LOG("Spline", pDEBUG) << "tag: " << tag << ", value: " << val;
           }//if current tag is <E>,<xsec>

           xmlKnotChild = xmlKnotChild->next;
        }//[end of] loop over tags within <knot>...</knot> tags
        xmlFree(xmlKnotChild);
        iknot++;
     } // if <knot>
     xmlSplChild = xmlSplChild->next;
  } //[end of] loop over tags within <spline>...</spline> tags
  xmlFree(xmlSplChild);

  this->BuildSpline(nknots, vx, vy);
  delete [] vx;
  delete [] vy;

  return true;
}

void Spline::Multiply	(	const Spline &	spl,
		double	c = 1
	)

Definition at line 612 of file Spline.cxx.

References BuildSpline(), Evaluate(), GetKnot(), IsWithinValidRange(), LOG, NKnots(), pERROR, ResetSpline(), XMax(), and XMin().

{
  // continue only if the input spline is defined at a wider x-range
  double xmin = this->XMin();
  double xmax = this->XMax();
  bool ok = spl.IsWithinValidRange(xmin) && spl.IsWithinValidRange(xmax);
  if(!ok) {
    LOG("Spline", pERROR) << "** Can't multiply splines. X-range mismatch!";
    return;
  }

  int nknots = this->NKnots();
  double * x = new double[nknots];
  double * y = new double[nknots];

  for(int i=0; i<nknots; i++) {
    this->GetKnot(i,x[i],y[i]);
    y[i] *= (c * spl.Evaluate(x[i]));
  }
  this->ResetSpline();
  this->BuildSpline(nknots,x,y);
  delete [] x;
  delete [] y;
}

void Spline::Multiply

(

double

a

)

Definition at line 686 of file Spline.cxx.

References a, BuildSpline(), GetKnot(), NKnots(), and ResetSpline().

{
  int nknots = this->NKnots();
  double * x = new double[nknots];
  double * y = new double[nknots];

  for(int i=0; i<nknots; i++) {
    this->GetKnot(i,x[i],y[i]);
    y[i]*=a;
  }
  this->ResetSpline();
  this->BuildSpline(nknots,x,y);
  delete [] x;
  delete [] y;
}

string genie::Spline::Name ( void  ) const

inline

Definition at line 95 of file Spline.h.

References fName.

Referenced by Print().

{ return fName; }

int genie::Spline::NKnots ( void  ) const

inline

Definition at line 84 of file Spline.h.

References fNKnots.

Referenced by Add(), Divide(), Multiply(), Print(), SaveAsText(), SaveAsXml(), and Spline().

{return fNKnots;}

void Spline::Print

(

ostream &

stream

)

const

Definition at line 569 of file Spline.cxx.

References GetKnot(), Name(), NKnots(), XMax(), and XMin().

Referenced by genie::operator<<().

{
  int    nknots = this->NKnots();
  double xmin   = this->XMin();
  double xmax   = this->XMax();

  stream << endl;
  stream << "** Spline: " << this->Name() << endl;
  stream << "Has " << nknots
        << " knots in the [" << xmin << ", " << xmax << "] range" << endl;
  double x,y;
  for(int i=0; i<nknots; i++) {
    this->GetKnot(i,x,y);
    stream << "-- knot : " << i
           << " -> (x = " << x << ", y = " << y << ")" << endl;
  }
}

void Spline::ResetSpline ( void  )

private

Definition at line 742 of file Spline.cxx.

References fGSLInterpolator, fInterpolator, fInterpolator5, and InitSpline().

Referenced by Add(), Divide(), and Multiply().

{
  if(fInterpolator) delete fInterpolator;
  if(fInterpolator5) delete fInterpolator5;
  if(fGSLInterpolator) delete fGSLInterpolator;
  this->InitSpline();
}

void Spline::SaveAsROOT	(	string	filename,
		string	name = "",
		bool	recreate = false
	)		const

Definition at line 485 of file Spline.cxx.

References fInterpolator, and fName.

{
  string spline_name = (name.size()>0 ? name : fName);

  string opt = ( (recreate) ? "RECREATE" : "UPDATE" );

  TFile f(filename.c_str(), opt.c_str());
  if(fInterpolator) fInterpolator->Write(spline_name.c_str());
  f.Close();
}

void Spline::SaveAsText	(	string	filename,
		string	format = "%10.6f\t%10.6f"
	)		const

Definition at line 464 of file Spline.cxx.

References fInterpolator, LOG, NKnots(), and pERROR.

{
  ofstream outtxt(filename.c_str());
  if(!outtxt.is_open()) {
    LOG("Spline", pERROR) << "Couldn't create file = " << filename;
    return;
  }
  int nknots = this->NKnots();
  outtxt << nknots << endl;

  double x=0, y=0;
  for(int iknot = 0; iknot < nknots; iknot++) {
    fInterpolator->GetKnot(iknot, x, y);
    char line[1024];
    sprintf(line,format.c_str(),x,y);
    outtxt << line << endl;
  }
  outtxt << endl;
  outtxt.close();
}

void Spline::SaveAsXml	(	string	filename,
		string	xtag,
		string	ytag,
		string	name = ""
	)		const

Definition at line 419 of file Spline.cxx.

References fName, LOG, and pERROR.

Referenced by genie::XSecSplineList::SaveAsXml().

{
  ofstream outxml(filename.c_str());
  if(!outxml.is_open()) {
    LOG("Spline", pERROR) << "Couldn't create file = " << filename;
    return;
  }
  string spline_name = (name.size()>0 ? name : fName);

  this->SaveAsXml(outxml, xtag, ytag, spline_name);

  outxml << endl;
  outxml.close();
}

void Spline::SaveAsXml	(	ofstream &	str,
		string	xtag,
		string	ytag,
		string	name = ""
	)		const

Definition at line 435 of file Spline.cxx.

References fInterpolator, fName, and NKnots().

{
  string spline_name = (name.size()>0 ? name : fName);

  // create a spline tag with the number of knots as an attribute
  int nknots = this->NKnots();
  string padding = "    ";
  ofs << padding << "<spline name=\"" << spline_name
      << "\" nknots=\"" << nknots << "\">" << endl;

  // start printing the knots
  double x=0, y=0;
  for(int iknot = 0; iknot < nknots; iknot++)
  {
    fInterpolator->GetKnot(iknot, x, y);

    ofs  << std::fixed << setprecision(5);
    ofs  << "\t<knot>"
         << " <" << xtag << "> " << setfill(' ')
                                 << setw(10) << x << " </" << xtag << ">";
    ofs  << std::scientific << setprecision(10);
    ofs  << " <" << ytag << "> " << setfill(' ')
                                 << setw(10) << y << " </" << ytag << ">"
         << " </knot>" << endl;
  }
  ofs << padding << "</spline>" << endl;
}

void genie::Spline::SetName ( string  name )

inline

Definition at line 94 of file Spline.h.

References fName.

{ fName = name; }

void Spline::SetType

(

string

type

)

Definition at line 769 of file Spline.cxx.

References fGSLInterpolator, fInterpolator, fInterpolator5, fInterpolatorType, fNKnots, LOG, pWARN, and genie::utils::str::ToUpper().

Referenced by genie::CacheBranchFx::CreateSpline().

{
  if(!fInterpolator) return;
  
  fInterpolatorType = genie::utils::str::ToUpper(type);
  
  ROOT::Math::Interpolation::Type gsltype;
  
  if ( fInterpolatorType ==  "TSPLINE3" || fInterpolatorType ==  "" ) 
  {
    fInterpolatorType = "TSpline3";
    return;
  }
  else if ( fInterpolatorType == "TSPLINE5" )
  {
    fInterpolatorType = "TSpline5";
    if(fInterpolator5) delete fInterpolator5;
    double x[fNKnots], y[fNKnots];
    for (int i=0; i<fNKnots; i++)
    {
       fInterpolator->GetKnot(i, x[i], y[i]);
    }
    fInterpolator5 = new TSpline5("spl5", x, y, fNKnots, "0");
    return;
  }
  else if ( fInterpolatorType ==  "LINEAR" )
  {
    gsltype = ROOT::Math::Interpolation::kLINEAR;
  }
  else if ( fInterpolatorType == "POLYNOMIAL" )
  {
    gsltype = ROOT::Math::Interpolation::kPOLYNOMIAL;
  }
  else if ( fInterpolatorType == "CSPLINE" )
  {
    gsltype = ROOT::Math::Interpolation::kCSPLINE;
  }
  else if ( fInterpolatorType == "CSPLINE_PERIODIC" )
  {
    gsltype = ROOT::Math::Interpolation::kCSPLINE_PERIODIC;
  }
  else if ( fInterpolatorType == "AKIMA" )
  {
    gsltype = ROOT::Math::Interpolation::kAKIMA;
  }
  else if ( fInterpolatorType == "AKIMA_PERIODIC" )
  {
    gsltype = ROOT::Math::Interpolation::kAKIMA_PERIODIC;
  }
  else
  {
    fInterpolatorType = "TSpline3";
    LOG("Spline", pWARN)
       << "Unknown interpolator type. Setting it to default [TSpline3].";
    return;
  }
  
  if(fGSLInterpolator) delete fGSLInterpolator;
  vector<double> x(fNKnots);
  vector<double> y(fNKnots);
  for (int i=0; i<fNKnots; i++)
  {
     fInterpolator->GetKnot(i, x[i], y[i]);
  }
  fGSLInterpolator = new ROOT::Math::Interpolator(x, y, gsltype);

   
}

double genie::Spline::XMax ( void  ) const

inline

Definition at line 89 of file Spline.h.

References fXMax.

Referenced by Add(), genie::KineGeneratorWithCache::CacheMaxXSec(), Divide(), genie::KineGeneratorWithCache::FindMaxXSec(), Multiply(), and Print().

{return fXMax;  }

double genie::Spline::XMin ( void  ) const

inline

Definition at line 88 of file Spline.h.

References fXMin.

Referenced by Add(), Divide(), genie::KineGeneratorWithCache::FindMaxXSec(), Multiply(), and Print().

{return fXMin;  }

void genie::Spline::YCanBeNegative ( bool  tf )

inline

Definition at line 97 of file Spline.h.

References fYCanBeNegative.

{ fYCanBeNegative = tf; }

double genie::Spline::YMax ( void  ) const

inline

Definition at line 90 of file Spline.h.

References fYMax.

{return fYMax;  }

Friends And Related Function Documentation

ostream& operator<< ( ostream &  stream, const Spline &  spl  )

friend

Definition at line 45 of file Spline.cxx.

  {
     spl.Print(stream);
     return stream;
  }

Member Data Documentation

ROOT::Math::Interpolator* genie::Spline::fGSLInterpolator

private

Definition at line 146 of file Spline.h.

Referenced by Evaluate(), InitSpline(), ResetSpline(), SetType(), and ~Spline().

TSpline3* genie::Spline::fInterpolator

private

Definition at line 143 of file Spline.h.

Referenced by BuildSpline(), Evaluate(), FindClosestKnot(), GetAsTSpline(), GetKnot(), GetKnotX(), GetKnotY(), InitSpline(), ResetSpline(), SaveAsROOT(), SaveAsText(), SaveAsXml(), SetType(), and ~Spline().

TSpline5* genie::Spline::fInterpolator5

private

Definition at line 145 of file Spline.h.

Referenced by Evaluate(), InitSpline(), ResetSpline(), SetType(), and ~Spline().

string genie::Spline::fInterpolatorType

private

Definition at line 147 of file Spline.h.

Referenced by Evaluate(), GetType(), InitSpline(), and SetType().

string genie::Spline::fName

private

Definition at line 138 of file Spline.h.

Referenced by InitSpline(), Name(), SaveAsROOT(), SaveAsXml(), and SetName().

int genie::Spline::fNKnots

private

Definition at line 139 of file Spline.h.

Referenced by BuildSpline(), NKnots(), and SetType().

double genie::Spline::fXMax

private

Definition at line 141 of file Spline.h.

Referenced by BuildSpline(), Evaluate(), GetAsTGraph(), InitSpline(), IsWithinValidRange(), and XMax().

double genie::Spline::fXMin

private

Definition at line 140 of file Spline.h.

Referenced by BuildSpline(), Evaluate(), GetAsTGraph(), InitSpline(), IsWithinValidRange(), and XMin().

bool genie::Spline::fYCanBeNegative

private

Definition at line 144 of file Spline.h.

Referenced by Evaluate(), InitSpline(), and YCanBeNegative().

double genie::Spline::fYMax

private

Definition at line 142 of file Spline.h.

Referenced by BuildSpline(), InitSpline(), and YMax().

---

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

• Spline.h
• Spline.cxx