#include "Samples/HistogramUtils.h"
#include "TList.h"
#include "TObjArray.h"

Include dependency graph for HistogramUtils.cpp:

Namespaces
namespace	M3

Functions
_MaCh3_Safe_Include_Start_ _MaCh3_Safe_Include_End_ void	CheckTH2PolyFileVersion (TFile *file)
	KS: ROOT changes something with binning when moving from ROOT 5 to ROOT 6. If you open ROOT5 produced file with ROOT6 you will be missing 9 last bins.

double	OverflowIntegral (TH2Poly *poly)
	WP: Helper function for calculating unbinned Integral of TH2Poly i.e including overflow.

double	NoOverflowIntegral (TH2Poly *poly)
	WP: Helper function for calculating binned Integral of TH2Poly i.e not including overflow.

TH1D *	PolyProjectionX (TObject *poly, std::string TempName, const std::vector< double > &xbins, const bool computeErrors)
	WP: Poly Projectors.

TH1D *	PolyProjectionY (TObject *poly, std::string TempName, const std::vector< double > &ybins, const bool computeErrors)
	WP: Poly Projectors.

TH2Poly *	NormalisePoly (TH2Poly *Histogram)
	WP: Helper to Normalise histograms.

void	NormaliseTH2Poly (TH2Poly *Histogram)
	Helper to Normalise histograms.

template<class HistType >
HistType *	RatioHists (HistType NumHist, HistType DenomHist)
	Helper to make ratio histograms.

TH2Poly *	RatioPolys (TH2Poly NumHist, TH2Poly DenomHist)
	Helper to make ratio of TH2Polys.

TH2D *	ConvertTH2PolyToTH2D (TH2Poly poly, TH2D h2dhist)
	KS: Convert TH2D to TH2Poly.

TH2Poly *	ConvertTH2DtoTH2Poly (TH2D *hist)
	KS: Convert TH2Poly to TH2D.

TH2Poly *	PolyScaleWidth (TH2Poly *Histogram, double scale)
	WP: Helper to scale th2poly analogous to th2d scale with option "width".

double	PolyIntegralWidth (TH2Poly *Histogram)
	WP: Helper to calc integral of th2poly analogous to th2d integra; with option "width".

TH2Poly *	MakePolyHist (const std::string &name, const std::vector< double > &BinArray_x, const std::vector< double > &BinArray_y)
	WP: Helper function to create TH2Poly histogram with uniform binning.

void	RemoveFitter (TH1D *hist, const std::string &name)
	KS: Remove fitted TF1 from hist to make comparison easier.

void	MakeFluctuatedHistogramStandard (TH1D FluctHist, TH1D PolyHist, TRandom3 *rand)
	Make Poisson fluctuation of TH1D hist using default fast method.

void	MakeFluctuatedHistogramStandard (TH2Poly FluctHist, TH2Poly PolyHist, TRandom3 *rand)
	Make Poisson fluctuation of TH2Poly hist using default fast method.

void	MakeFluctuatedHistogramAlternative (TH1D FluctHist, TH1D PolyHist, TRandom3 *rand)
	Make Poisson fluctuation of TH1D hist using slow method which is only for cross-check.

int	GetRandomPoly2 (const TH2Poly PolyHist, TRandom3 rand)
	KS: ROOT developers were too lazy do develop getRanom2 for TH2Poly, this implementation is based on link

void	MakeFluctuatedHistogramAlternative (TH2Poly FluctHist, TH2Poly PolyHist, TRandom3 *rand)
	Make Poisson fluctuation of TH2Poly hist using slow method which is only for cross-check.

std::unique_ptr< TGraphAsymmErrors >	MakeAsymGraph (TH1 sigmaArrayLeft, TH1 sigmaArrayCentr, TH1 *sigmaArrayRight, const std::string &title)
	Used by sigma variation, check how 1 sigma changes spectra.

void	FastViolinFill (TH2D violin, TH1D hist_1d)
	KS: Fill Violin histogram with entry from a toy.

double	returnCherenkovThresholdMomentum (int PDG)
	DB Get the Cherenkov momentum threshold in MeV.

double	CalculateQ2 (double PLep, double PUpd, double EnuTrue, double InitialQ2)
	Recalculate Q^2 after Eb shift. Takes in shifted lepton momentum, lepton angle, and true neutrino energy.

double	CalculateEnu (double PLep, double costh, double Eb, bool neutrino)
	Recalculate Enu after Eb shift. Takes in shifted lepton momentum, lepton angle, and binding energy change, and if nu/anu.

TFile *	M3::Open (const std::string &Name, const std::string &Type, const std::string &File, const int Line)
	Opens a ROOT file with the given name and mode.

Function Documentation

◆ CalculateEnu()

double CalculateEnu	(	double	PLep,
		double	costh,
		double	Eb,
		bool	neutrino
	)

Recalculate Enu after Eb shift. Takes in shifted lepton momentum, lepton angle, and binding energy change, and if nu/anu.

Todo:: WARNING this is hardcoded

Definition at line 589 of file HistogramUtils.cpp.

                                                                        {
// ***************************************************************************
  double mNeff = 0.93956536 - Eb / 1000.;
  double mNoth = 0.93827203;
 
  if (!neutrino) {
    mNeff = 0.93827203 - Eb / 1000.;
    mNoth = 0.93956536;
  }
  constexpr double mLep = 0.10565837;
  double eLep = sqrt(PLep * PLep + mLep * mLep);
 
  double Enu = (2 * mNeff * eLep - mLep * mLep + mNoth * mNoth - mNeff * mNeff) /(2 * (mNeff - eLep + PLep * costh));
 
  return Enu;
}

◆ CalculateQ2()

double CalculateQ2	(	double	PLep,
		double	PUpd,
		double	EnuTrue,
		double	InitialQ2
	)

Recalculate Q^2 after Eb shift. Takes in shifted lepton momentum, lepton angle, and true neutrino energy.

Definition at line 570 of file HistogramUtils.cpp.

                                                                              {
// ***************************************************************************
  constexpr double MLep = 0.10565837;
 
  // Caluclate muon energy
  double ELep = sqrt((MLep*MLep)+(PLep*PLep));
 
  double CosTh = (2*EnuTrue*ELep - MLep*MLep - InitialQ2)/(2*EnuTrue*PLep);
 
  ELep = sqrt((MLep*MLep)+(PUpd*PUpd));
 
  // Calculate the new Q2
  double Q2Upd = -(MLep*MLep) + 2.0*EnuTrue*(ELep - PUpd*CosTh);
 
  return Q2Upd - InitialQ2;
}

◆ CheckTH2PolyFileVersion()

_MaCh3_Safe_Include_Start_ _MaCh3_Safe_Include_End_ void CheckTH2PolyFileVersion ( TFile * file )

KS: ROOT changes something with binning when moving from ROOT 5 to ROOT 6. If you open ROOT5 produced file with ROOT6 you will be missing 9 last bins.

However if you use ROOT6 and have ROOT6 file exactly the same code will work. Something have changed with how TH2Poly bins are stored in TFile

Parameters

file	ROOT file that we will make version checks

Definition at line 11 of file HistogramUtils.cpp.

                                          {
// **************************************************
  int FileROOTVersion = file->GetVersion();
  int MainFileROOTVersion = FileROOTVersion;
 
  // Remove last digit from number
  // till only one digit is left
  while (MainFileROOTVersion >= 10)
      MainFileROOTVersion /= 10;
 
  std::string SystemROOTVersion = std::string(ROOT_RELEASE);
  int MainSystemROOTVersion = SystemROOTVersion.at(0)  - '0';
 
  if(MainFileROOTVersion != MainSystemROOTVersion)
  {
    MACH3LOG_ERROR("File was produced with: {} ROOT version", FileROOTVersion);
    MACH3LOG_ERROR("Found: {} ROOT version in the system", SystemROOTVersion);
    MACH3LOG_ERROR("For some documentation please visit me");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

◆ ConvertTH2DtoTH2Poly()

TH2Poly * ConvertTH2DtoTH2Poly ( TH2D * hist )

KS: Convert TH2Poly to TH2D.

Definition at line 282 of file HistogramUtils.cpp.

                                          {
// ****************
  // Make the x axis from the momentum of lepton
  TAxis* xaxis = hist->GetXaxis();
  // Make the y axis from the cos theta of lepton
  TAxis* yaxis = hist->GetYaxis();
 
  TString histname = hist->GetName();
  // Convert TH2D binning to TH2Poly
  TH2Poly* poly = new TH2Poly();
  poly->SetName(histname);
  poly->SetTitle(histname);
 
  // Copy axis titles
  poly->GetXaxis()->SetTitle(xaxis->GetTitle());
  poly->GetYaxis()->SetTitle(yaxis->GetTitle());
 
  double xmax, xmin, ymax, ymin;
  for (int iy = 1; iy <= yaxis->GetNbins(); iy++) {
    ymax = yaxis->GetBinUpEdge(iy);
    ymin = yaxis->GetBinLowEdge(iy);
    for (int ix = 1; ix <= xaxis->GetNbins(); ix++) {
      xmax = xaxis->GetBinUpEdge(ix);
      xmin = xaxis->GetBinLowEdge(ix);
      double binofx[] = {xmin, xmax, xmax, xmin};
      double binofy[] = {ymin, ymin, ymax, ymax};
      poly->AddBin(4, binofx, binofy);
    }
  }
 
  for (int iy = 1; iy <= yaxis->GetNbins(); iy++) {
    ymax = yaxis->GetBinUpEdge(iy);
    ymin = yaxis->GetBinLowEdge(iy);
    for (int ix = 1; ix <= xaxis->GetNbins(); ix++) {
      xmax = xaxis->GetBinUpEdge(ix);
      xmin = xaxis->GetBinLowEdge(ix);
 
      // Get the content of the corresponding bin in TH2D and set it in TH2Poly
      int bin = hist->GetBin(ix, iy);
      double content = hist->GetBinContent(bin);
      poly->SetBinContent(poly->FindBin((xmin + xmax) / 2, (ymin + ymax) / 2), content);
    }
  }
 
  return poly;
}

◆ ConvertTH2PolyToTH2D()

TH2D * ConvertTH2PolyToTH2D	(	TH2Poly *	poly,
		TH2D *	h2dhist
	)

KS: Convert TH2D to TH2Poly.

Definition at line 248 of file HistogramUtils.cpp.

                                                         {
// ****************
  double xlow, xup, ylow, yup;
  std::string HistTempName = poly->GetName();
 
  HistTempName += "_";
  //make the th2d
  TH2D *hist = static_cast<TH2D*>(h2dhist->Clone());
  hist->SetNameTitle(HistTempName.c_str(), HistTempName.c_str());
 
  for(int ix = 0; ix < hist->GetNbinsX() + 2; ix++) {
    for(int iy = 0; iy < hist->GetNbinsY() + 2; iy++) {
      hist->SetBinContent(ix, iy, 0);
    }
  }
  //Loop over poly bins, find the corresponding th2d and setbincontent!
  for(int i = 0; i< poly->GetNumberOfBins(); i++){
    TH2PolyBin & polybin = static_cast<TH2PolyBin &>(*poly->GetBins()->At(i));
    xlow = polybin.GetXMin();
    xup = polybin.GetXMax();
    ylow = polybin.GetYMin();
    yup = polybin.GetYMax();
    int xbin, ybin;
 
    xbin = hist->GetXaxis()->FindBin(xlow+(xup-xlow)/2);
    ybin = hist->GetYaxis()->FindBin(ylow+(yup-ylow)/2);
 
    MACH3LOG_TRACE("Poly bin {}, xlow: {}, xup: {}, ylow: {}, yup: {}. Finding bin for ({}, {}). Found Bin ({}, {}) with content {}. But Poly content: {}",
                i, xlow, xup, ylow, yup, (xlow + (xup - xlow) / 2), (ylow + (yup - ylow) / 2), xbin, ybin, polybin.GetContent(), poly->GetBinContent(i));
    hist->SetBinContent(xbin, ybin, polybin.GetContent());
  }
  return hist;
}

◆ FastViolinFill()

void FastViolinFill	(	TH2D *	violin,
		TH1D *	hist_1d
	)

KS: Fill Violin histogram with entry from a toy.

Parameters

violin	hist that will be filled
hist_1d	refence hist from which we take entries to be filled

Definition at line 548 of file HistogramUtils.cpp.

                                                {
// ****************
  for (int x = 0; x < violin->GetXaxis()->GetNbins(); ++x)
  {
    const int y = violin->GetYaxis()->FindBin(hist_1d->GetBinContent(x+1));
    violin->SetBinContent(x+1, y,  violin->GetBinContent(x+1, y)+1);
  }
}

◆ GetRandomPoly2()

int GetRandomPoly2	(	const TH2Poly *	PolyHist,
		TRandom3 *	rand
	)

KS: ROOT developers were too lazy do develop getRanom2 for TH2Poly, this implementation is based on link

Definition at line 466 of file HistogramUtils.cpp.

                                                           {
// ****************
  const int nbins = PolyHist->GetNumberOfBins();
  const double r1 = rand->Rndm();
 
  double* fIntegral = new double[nbins+2];
  fIntegral[0] = 0.0;
 
  //KS: This is custom version of ComputeIntegral, once again ROOT was lazy :(
  for (int i = 1; i < nbins+1; ++i)
  {
    fIntegral[i] = 0.0;
    const double content = PolyHist->GetBinContent(i);
    fIntegral[i] += fIntegral[i - 1] + content;
  }
  for (Int_t bin = 1; bin < nbins+1; ++bin)  fIntegral[bin] /= fIntegral[nbins];
  fIntegral[nbins+1] = PolyHist->GetEntries();
 
  //KS: We just return one rather then X and Y, this way we can use SetBinContent rather than Fill, which is faster
  int iBin = int(TMath::BinarySearch(nbins, fIntegral, r1));
  //KS: Have to increment because TH2Poly has stupid offset arghh
  iBin += 1;
 
  delete[] fIntegral;
  return iBin;
}

◆ MakeAsymGraph()

std::unique_ptr< TGraphAsymmErrors > MakeAsymGraph	(	TH1 *	sigmaArrayLeft,
		TH1 *	sigmaArrayCentr,
		TH1 *	sigmaArrayRight,
		const std::string &	title
	)

Used by sigma variation, check how 1 sigma changes spectra.

Parameters

sigmaArrayLeft	sigma var hist at -1 or -3 sigma shift
sigmaArrayCentr	sigma var hist at prior values
sigmaArrayRight	sigma var hist at +1 or +3 sigma shift
title	A tittle for returned object

Returns: A TGraphAsymmErrors object that visualizes the sigma variation of spectra, showing confidence intervals between different sigma shifts.

Definition at line 516 of file HistogramUtils.cpp.

                                                                                                                                        {
// *************************
  auto var = std::make_unique<TGraphAsymmErrors>(sigmaArrayCentr);
  var->SetNameTitle((title).c_str(), (title).c_str());
 
  // Need to draw TGraphs to set axes labels
  var->Draw("AP");
  var->GetXaxis()->SetTitle(sigmaArrayCentr->GetXaxis()->GetTitle());
  var->GetYaxis()->SetTitle("Number of events/bin");
 
  for (int m = 0; m < var->GetN(); ++m)
  {
    double xlow = sigmaArrayLeft->GetBinContent(m+1);
    double xhigh = sigmaArrayRight->GetBinContent(m+1);
    double xtemp;
 
    // Figure out which variation is larger so we set the error correctly
    if (xlow > xhigh)
    {
      xtemp = xlow;
      xlow = xhigh;
      xhigh = xtemp;
    }
 
    var->SetPointEYhigh(m, xhigh - var->GetY()[m]);
    var->SetPointEYlow(m, var->GetY()[m] - xlow);
  }
  return var;
}

◆ MakeFluctuatedHistogramAlternative() [1/2]

void MakeFluctuatedHistogramAlternative	(	TH1D *	FluctHist,
		TH1D *	PolyHist,
		TRandom3 *	rand
	)

Make Poisson fluctuation of TH1D hist using slow method which is only for cross-check.

Definition at line 443 of file HistogramUtils.cpp.

                                                                                        {
// ****************
  // Make the Poisson fluctuated hist
  FluctHist->Reset("");
  FluctHist->Fill(0.0, 0.0);
 
  const double evrate = PolyHist->Integral();
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wconversion"
  const int num = rand->Poisson(evrate);
  #pragma GCC diagnostic pop
  int count = 0;
  while(count < num)
  {
    const double candidate = PolyHist->GetRandom();
    FluctHist->Fill(candidate);
    count++;
  }
}

◆ MakeFluctuatedHistogramAlternative() [2/2]

void MakeFluctuatedHistogramAlternative	(	TH2Poly *	FluctHist,
		TH2Poly *	PolyHist,
		TRandom3 *	rand
	)

Make Poisson fluctuation of TH2Poly hist using slow method which is only for cross-check.

Definition at line 495 of file HistogramUtils.cpp.

                                                                                              {
// ****************
  // Make the Poisson fluctuated hist
  FluctHist->Reset("");
  FluctHist->Fill(0.0, 0.0, 0.0);
 
  const double evrate = NoOverflowIntegral(PolyHist);
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wconversion"
  const int num = rand->Poisson(evrate);
  #pragma GCC diagnostic pop
  int count = 0;
  while(count < num)
  {
    const int iBin = GetRandomPoly2(PolyHist, rand);
    FluctHist->SetBinContent(iBin, FluctHist->GetBinContent(iBin) + 1);
    count++;
  }
}

◆ MakeFluctuatedHistogramStandard() [1/2]

void MakeFluctuatedHistogramStandard	(	TH1D *	FluctHist,
		TH1D *	PolyHist,
		TRandom3 *	rand
	)

Make Poisson fluctuation of TH1D hist using default fast method.

Definition at line 405 of file HistogramUtils.cpp.

                                                                                     {
// ****************
  // Make the Poisson fluctuated hist
  FluctHist->Reset("");
  FluctHist->Fill(0.0, 0.0);
 
  for (int i = 1; i <= PolyHist->GetXaxis()->GetNbins(); ++i)
  {
    // Get the posterior predictive bin content
    const double MeanContent = PolyHist->GetBinContent(i);
    // Get a Poisson fluctuation of the content
    const double Random = rand->PoissonD(MeanContent);
    // Set the fluctuated histogram content to the Poisson variation of the posterior predictive histogram
    FluctHist->SetBinContent(i,Random);
  }
}

◆ MakeFluctuatedHistogramStandard() [2/2]

void MakeFluctuatedHistogramStandard	(	TH2Poly *	FluctHist,
		TH2Poly *	PolyHist,
		TRandom3 *	rand
	)

Make Poisson fluctuation of TH2Poly hist using default fast method.

Definition at line 424 of file HistogramUtils.cpp.

                                                                                            {
// ****************
  // Make the Poisson fluctuated hist
  FluctHist->Reset("");
  FluctHist->Fill(0.0, 0.0, 0.0);
 
  for (int i = 1; i < FluctHist->GetNumberOfBins()+1; ++i)
  {
    // Get the posterior predictive bin content
    const double MeanContent = PolyHist->GetBinContent(i);
    // Get a Poisson fluctuation of the content
    const double Random = rand->PoissonD(MeanContent);
    // Set the fluctuated histogram content to the Poisson variation of the posterior predictive histogram
    FluctHist->SetBinContent(i,Random);
  }
}

◆ MakePolyHist()

TH2Poly * MakePolyHist	(	const std::string &	name,
		const std::vector< double > &	BinArray_x,
		const std::vector< double > &	BinArray_y
	)

WP: Helper function to create TH2Poly histogram with uniform binning.

Parameters

name	This will be tittle of output histogram
BinArray_x	Bin edges for X axis
BinArray_y	Bin edges for Y axis

Definition at line 370 of file HistogramUtils.cpp.

                                                                                                                       {
// *********************
  TH2Poly* poly = new TH2Poly();
  poly->SetName(name.c_str());
  poly->SetTitle(name.c_str());
  double xmax, xmin, ymax, ymin;
  for(unsigned int iy = 0; iy < BinArray_y.size()-1; iy++)
  {
    ymax = BinArray_y[iy+1];
    ymin = BinArray_y[iy];
    for(unsigned int ix = 0; ix < BinArray_x.size()-1; ix++)
    {
      xmax = BinArray_x[ix+1];
      xmin = BinArray_x[ix];
      double binofx[] = {xmin, xmax, xmax, xmin};
      double binofy[] = {ymin, ymin, ymax, ymax};
      poly->AddBin(4,binofx,binofy);
    }
  }
  return poly;
}

◆ NoOverflowIntegral()

double NoOverflowIntegral ( TH2Poly * poly )

WP: Helper function for calculating binned Integral of TH2Poly i.e not including overflow.

Definition at line 50 of file HistogramUtils.cpp.

                                         {
// **************************************************
  double integral = 0.;
  for(int i = 1; i < poly->GetNumberOfBins()+1; i++)
  {
    integral += poly->GetBinContent(i);
  }
 
  return integral;
} // end function

◆ NormalisePoly()

TH2Poly * NormalisePoly ( TH2Poly * Histogram )

WP: Helper to Normalise histograms.

Definition at line 198 of file HistogramUtils.cpp.

                                           {
// ****************
  TH2Poly* HistCopy = static_cast<TH2Poly*>(Histogram->Clone());
  double IntegralWidth = PolyIntegralWidth(HistCopy);
  HistCopy = PolyScaleWidth(HistCopy, IntegralWidth);
  std::string title = std::string(HistCopy->GetName())+"_norm";
  HistCopy->SetNameTitle(title.c_str(), title.c_str());
 
  return HistCopy;
}

◆ NormaliseTH2Poly()

void NormaliseTH2Poly ( TH2Poly * Histogram )

Helper to Normalise histograms.

Parameters

Histogram hist which we normalise

Definition at line 211 of file HistogramUtils.cpp.

                                         {
// ****************
  const double Integral = NoOverflowIntegral(Histogram);
  for(int j = 1; j < Histogram->GetNumberOfBins()+1; j++)
  {
    Histogram->SetBinContent(j, Histogram->GetBinContent(j)/Integral);
  }
}

◆ OverflowIntegral()

double OverflowIntegral ( TH2Poly * poly )

WP: Helper function for calculating unbinned Integral of TH2Poly i.e including overflow.

Definition at line 35 of file HistogramUtils.cpp.

                                       {
// **************************************************
  double overflow = 0.;
  //TH2Polys have 9 overflow bins
  for(int iOverflow = -1; iOverflow > -10; iOverflow--)
  {
    overflow += poly->GetBinContent(iOverflow);
  }
  double IntegralUn = NoOverflowIntegral(poly) + overflow;
 
  return IntegralUn;
} // end function

◆ PolyIntegralWidth()

double PolyIntegralWidth ( TH2Poly * Histogram )

WP: Helper to calc integral of th2poly analogous to th2d integra; with option "width".

Definition at line 351 of file HistogramUtils.cpp.

                                             {
// ****************
  double integral = 0;
  double xlow, xup, ylow, yup, area;
 
  for(int i = 1; i < Histogram->GetNumberOfBins()+1; i++)
  {
    TH2PolyBin* bin = static_cast<TH2PolyBin*>(Histogram->GetBins()->At(i - 1));
    xlow = bin->GetXMin();
    xup = bin->GetXMax();
    ylow = bin->GetYMin();
    yup = bin->GetYMax();
    area = (xup-xlow)*(yup-ylow);
    integral += Histogram->GetBinContent(i)*area;
  }
  return integral;
}

◆ PolyProjectionX()

TH1D * PolyProjectionX	(	TObject *	poly,
		std::string	TempName,
		const std::vector< double > &	xbins,
		const bool	computeErrors
	)

WP: Poly Projectors.

Definition at line 63 of file HistogramUtils.cpp.

                                                                                                                   {
// **************************************************
  TH1D* hProjX = new TH1D((TempName+"_x").c_str(),(TempName+"_x").c_str(), int(xbins.size()-1), &xbins[0]);
 
  //KS: Temp Histogram to store error, use double as this is thread safe
  std::vector<double> hProjX_Error(hProjX->GetXaxis()->GetNbins() + 1, 0.0);
  double xlow, xup, frac = 0.;
 
  //loop over bins in the poly
  for (int i = 0; i < static_cast<TH2Poly*>(poly)->GetNumberOfBins(); i++)
  {
    //get bin and its edges
    TH2PolyBin* bin = static_cast<TH2PolyBin*>(static_cast<TH2Poly*>(poly)->GetBins()->At(i));
    xlow = bin->GetXMin();
    xup = bin->GetXMax();
 
    //Loop over projected bins, find fraction of poly bin in each
    for(int dx = 0; dx < int(xbins.size()); dx++)
    {
      if(xbins[dx+1] <= xlow || xbins[dx] >= xup)
      {
        frac = 0;
      }
      else if(xbins[dx] <= xlow && xbins[dx+1] >= xup)
      {
        frac = 1;
      }
      else if(xbins[dx] <= xlow && xbins[dx+1] <= xup)
      {
        frac = (xbins[dx+1]-xlow)/(xup-xlow);
      }
      else if(xbins[dx] >= xlow && xbins[dx+1] >= xup)
      {
        frac = (xup-xbins[dx])/(xup-xlow);
      }
      else if(xbins[dx] >= xlow && xbins[dx+1] <= xup)
      {
        frac = (xbins[dx+1]-xbins[dx])/(xup-xlow);
      }
      else
      {
        frac = 0;
      }
      hProjX->SetBinContent(dx+1,hProjX->GetBinContent(dx+1)+frac*bin->GetContent());
      //KS: Follow ROOT implementation and sum up the variance
      if(computeErrors)
      {
        //KS: TH2PolyBin doesn't have GetError so we have to use TH2Poly,
        //but numbering of GetBinError is different than GetBins...
        const double Temp_Err = frac*static_cast<TH2Poly*>(poly)->GetBinError(i+1)*frac*static_cast<TH2Poly*>(poly)->GetBinError(i+1);
        hProjX_Error[dx+1] += Temp_Err;
      }
    }
  }
  //KS: The error is sqrt(summed variance)) https://root.cern.ch/doc/master/TH2_8cxx_source.html#l02266
  if(computeErrors)
  {
    for (int i = 1; i <= hProjX->GetXaxis()->GetNbins(); ++i)
    {
      const double Error = std::sqrt(hProjX_Error[i]);
      hProjX->SetBinError(i, Error);
    }
  }
  return hProjX;
} // end project poly X function

◆ PolyProjectionY()

TH1D * PolyProjectionY	(	TObject *	poly,
		std::string	TempName,
		const std::vector< double > &	ybins,
		const bool	computeErrors
	)

WP: Poly Projectors.

Definition at line 131 of file HistogramUtils.cpp.

                                                                                                                   {
// **************************************************
  TH1D* hProjY = new TH1D((TempName+"_y").c_str(),(TempName+"_y").c_str(),int(ybins.size()-1),&ybins[0]);
  //KS: Temp Histogram to store error, use double as this is thread safe
  std::vector<double> hProjY_Error(hProjY->GetXaxis()->GetNbins() + 1, 0.0);
  double ylow, yup, frac = 0.;
 
  //loop over bins in the poly
  for (int i = 0; i < static_cast<TH2Poly*>(poly)->GetNumberOfBins(); i++)
  {
    //get bin and its edges
    TH2PolyBin* bin = static_cast<TH2PolyBin*>(static_cast<TH2Poly*>(poly)->GetBins()->At(i));
    ylow = bin->GetYMin();
    yup = bin->GetYMax();
 
    //Loop over projected bins, find fraction of poly bin in each
    for(int dy = 0; dy < int(ybins.size()); dy++)
    {
      if(ybins[dy+1]<=ylow || ybins[dy] >= yup)
      {
        frac = 0;
      }
      else if(ybins[dy] <= ylow && ybins[dy+1] >= yup)
      {
        frac = 1;
      }
      else if(ybins[dy] <= ylow && ybins[dy+1] <= yup)
      {
        frac = (ybins[dy+1]-ylow)/(yup-ylow);
      }
      else if(ybins[dy] >= ylow && ybins[dy+1] >= yup)
      {
        frac = (yup-ybins[dy])/(yup-ylow);
      }
      else if(ybins[dy] >= ylow && ybins[dy+1] <= yup)
      {
        frac = (ybins[dy+1]-ybins[dy])/(yup-ylow);
      }
      else
      {
        frac = 0;
      }
      hProjY->SetBinContent(dy+1,hProjY->GetBinContent(dy+1)+frac*bin->GetContent());
      //KS: Follow ROOT implementation and sum up the variance
      if(computeErrors)
      {
        //KS: TH2PolyBin doesn't have GetError so we have to use TH2Poly,
        //but numbering of GetBinError is different than GetBins...
        const double Temp_Err = frac*static_cast<TH2Poly*>(poly)->GetBinError(i+1)*frac*static_cast<TH2Poly*>(poly)->GetBinError(i+1);
        hProjY_Error[dy+1] += Temp_Err;
      }
    }
  }
  //KS: The error is sqrt(summed variance)) https://root.cern.ch/doc/master/TH2_8cxx_source.html#l02266
  if(computeErrors)
  {
    for (int i = 1; i <= hProjY->GetXaxis()->GetNbins(); ++i)
    {
      const double Error = std::sqrt(hProjY_Error[i]);
      hProjY->SetBinError(i, Error);
    }
  }
  return hProjY;
} // end project poly Y function

◆ PolyScaleWidth()

TH2Poly * PolyScaleWidth	(	TH2Poly *	Histogram,
		double	scale
	)

WP: Helper to scale th2poly analogous to th2d scale with option "width".

Definition at line 331 of file HistogramUtils.cpp.

                                                          {
// ****************
  TH2Poly* HistCopy = static_cast<TH2Poly*>(Histogram->Clone());
  double xlow, xup, ylow, yup, area;
 
  for(int i = 1; i < HistCopy->GetNumberOfBins()+1; i++)
  {
    TH2PolyBin* bin = static_cast<TH2PolyBin*>(Histogram->GetBins()->At(i - 1));
    xlow = bin->GetXMin();
    xup = bin->GetXMax();
    ylow = bin->GetYMin();
    yup = bin->GetYMax();
    area = (xup-xlow)*(yup-ylow);
    HistCopy->SetBinContent(i, Histogram->GetBinContent(i)/(area*scale));
  }
  return HistCopy;
}

◆ RatioHists()

template<class HistType >

HistType * RatioHists	(	HistType *	NumHist,
		HistType *	DenomHist
	)

Helper to make ratio histograms.

Definition at line 223 of file HistogramUtils.cpp.

                                                             {
// ****************
  HistType *NumCopy = static_cast<HistType*>(NumHist->Clone());
  std::string title = std::string(DenomHist->GetName()) + "_ratio";
  NumCopy->SetNameTitle(title.c_str(), title.c_str());
  NumCopy->Divide(DenomHist);
 
  return NumCopy;
}

◆ RatioPolys()

TH2Poly * RatioPolys	(	TH2Poly *	NumHist,
		TH2Poly *	DenomHist
	)

Helper to make ratio of TH2Polys.

Definition at line 235 of file HistogramUtils.cpp.

                                                          {
// ****************
  TH2Poly *NumCopy = static_cast<TH2Poly*>(NumHist->Clone());
  std::string title = std::string(DenomHist->GetName()) + "_ratio";
  NumCopy->SetNameTitle(title.c_str(), title.c_str());
 
  for(int i = 1; i < NumCopy->GetNumberOfBins()+1; ++i) {
    NumCopy->SetBinContent(i,NumHist->GetBinContent(i)/DenomHist->GetBinContent(i));
  }
  return NumCopy;
}

◆ RemoveFitter()

void RemoveFitter	(	TH1D *	hist,
		const std::string &	name
	)

KS: Remove fitted TF1 from hist to make comparison easier.

Definition at line 394 of file HistogramUtils.cpp.

                                                     {
// *********************
  TList *listOfFunctions = hist->GetListOfFunctions();
  TF1 *fitter = dynamic_cast<TF1*>(listOfFunctions->FindObject(name.c_str()));
 
  listOfFunctions->Remove(fitter);
  delete fitter;
}

◆ returnCherenkovThresholdMomentum()

double returnCherenkovThresholdMomentum ( int PDG )

DB Get the Cherenkov momentum threshold in MeV.

Definition at line 560 of file HistogramUtils.cpp.

                                                 {
// ****************
  constexpr double refractiveIndex = 1.334; //DB From https://github.com/fiTQun/fiTQun/blob/646cf9c8ba3d4f7400bcbbde029d5ca15513a3bf/fiTQun_shared.cc#L757
  double mass =  MaCh3Utils::GetMassFromPDG(PDG)*1e3;
  double momentumThreshold = mass/sqrt(refractiveIndex*refractiveIndex-1.);
  return momentumThreshold;
}

Namespaces

Functions

Function Documentation

◆ CalculateEnu()

◆ CalculateQ2()

◆ CheckTH2PolyFileVersion()

◆ ConvertTH2DtoTH2Poly()

◆ ConvertTH2PolyToTH2D()

◆ FastViolinFill()

◆ GetRandomPoly2()

◆ MakeAsymGraph()

◆ MakeFluctuatedHistogramAlternative() [1/2]

◆ MakeFluctuatedHistogramAlternative() [2/2]

◆ MakeFluctuatedHistogramStandard() [1/2]

◆ MakeFluctuatedHistogramStandard() [2/2]

◆ MakePolyHist()

◆ NoOverflowIntegral()

◆ NormalisePoly()

◆ NormaliseTH2Poly()

◆ OverflowIntegral()

◆ PolyIntegralWidth()

◆ PolyProjectionX()

◆ PolyProjectionY()

◆ PolyScaleWidth()

◆ RatioHists()

◆ RatioPolys()

◆ RemoveFitter()

◆ returnCherenkovThresholdMomentum()