#include "Samples/SampleStructs.h"
#include "Parameters/ParameterStructs.h"
#include "Manager/YamlHelper.h"

Include dependency graph for HistogramUtils.h:

This graph shows which files directly or indirectly include this file:

Namespaces
	M3
	Main namespace for MaCh3 software.

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

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

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

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

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

TH2Poly *	ConvertTH2DtoTH2Poly (TH2D *TH2Dhist)
	KS: Convert TH2Poly to TH2D. More...

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

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

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

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

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

TH2Poly *	RatioPolys (TH2Poly NumPoly, TH2Poly DenomPoly)
	Helper to make ratio of TH2Polys. More...

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. More...

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. More...

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

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

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

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

void	MakeFluctuatedHistogramAlternative (TH2D FluctHist, TH2D PolyHist, TRandom3 *rand)
	Make Poisson fluctuation of TH2D hist. More...

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

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

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

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

template<typename Derived , typename Base >
std::vector< std::unique_ptr< Base > >	CastVector (std::vector< std::unique_ptr< Derived >> &&inputVec)
	Converts a vector of pointers from a derived type to a base type. More...

std::string	file_exists (std::string filename)
	Helper to check if files exist or not. More...

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

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

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

std::unique_ptr< TH1D >	MakeSummaryFromSpectra (const TH2D *Spectra, const std::string &name)
	Build a 1D posterior-predictive summary from a violin spectrum. More...

template<typename ObjectType >
std::unique_ptr< ObjectType >	M3::Clone (const ObjectType *obj, const std::string &name="")
	KS: Creates a copy of a ROOT-like object and wraps it in a smart pointer. More...

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. More...

void	M3::ScaleHistogram (TH1 *Sample_Hist, const double scale)
	Scale histogram to get divided by bin width. More...

void	M3::CheckBinningMatch (const TH1D Hist1, const TH1D Hist2, const std::string &File, const int Line)
	KS: Helper function check if data and MC binning matches. More...

void	M3::CheckBinningMatch (const TH2D Hist1, const TH2D Hist2, const std::string &File, const int Line)
	KS: Helper function check if data and MC binning matches. More...

void	M3::CheckBinningMatch (TH2Poly Hist1, TH2Poly Hist2, const std::string &File, const int Line)
	KS: Helper function check if data and MC binning matches. More...

YAML::Node	M3::PolyToYaml (TH2Poly *Hist, const std::string &YamlName, const std::string &File, const int Line)
	KS: Convert TH2Poly into yaml config accepted by MaCh3. More...

Detailed Description

Author: Will Parker; Kamil Skwarczynski

Definition in file HistogramUtils.h.

Function Documentation

◆ CalculateEnu()

double CalculateEnu	(	double	PLep,
		double	cosTheta,
		double	EB,
		bool	neutrino
	)

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

Parameters

PLep	Shifted outgoing lepton momentum (MeV/c).
cosTheta	Cosine of the lepton scattering angle.
EB	Binding energy shift applied to the interaction (MeV).
neutrino	True if the interaction is neutrino, false if antineutrino.

Todo:: WARNING this is hardcoded

Definition at line 604 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;
   constexpr double mLep2 = mLep * mLep;
  
   const double eLep = sqrt(PLep * PLep + mLep2);
   double Enu = (2 * mNeff * eLep - mLep2 + mNoth * mNoth - mNeff * mNeff) /(2 * (mNeff - eLep + PLep * costh));
  
   return Enu;
 }

◆ CalculateQ2()

double CalculateQ2	(	double	PLep,
		double	PUpd,
		double	EnuTrue,
		double	InitialQ2 = `0.0`
	)

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

Parameters

PLep	Shifted outgoing lepton momentum (MeV/c).
PUpd	Updated lepton angle or related kinematic quantity used in the recalculation.
EnuTrue	True neutrino energy (MeV).
InitialQ2	Optional initial Q^2 value (used as seed or fallback, default = 0).

Definition at line 584 of file HistogramUtils.cpp.

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

◆ CastVector()

template<typename Derived , typename Base >

std::vector<std::unique_ptr<Base> > CastVector ( std::vector< std::unique_ptr< Derived >> && inputVec )

Converts a vector of pointers from a derived type to a base type.

Template Parameters

Derived	The derived class type.
Base	The base class type.

Parameters

inputVec A std::vector of pointers to Derived objects.

Returns: A std::vector of pointers to Base objects.

Definition at line 95 of file HistogramUtils.h.

                                                          {
   std::vector<std::unique_ptr<Base>> outputVec;
   outputVec.reserve(inputVec.size());
   for (auto& ptr : inputVec) {
     outputVec.push_back(std::unique_ptr<Base>(std::move(ptr)));
   }
   return outputVec;
 }

◆ CheckTH2PolyFileVersion()

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 14 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 * TH2Dhist )

KS: Convert TH2Poly to TH2D.

Definition at line 285 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 *	TH2Dhist
	)

KS: Convert TH2D to TH2Poly.

Definition at line 251 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 562 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);
   }
 }

◆ file_exists()

std::string file_exists ( std::string filename )

inline

Helper to check if files exist or not.

Definition at line 105 of file HistogramUtils.h.

                                                  {
   std::ifstream infile(filename.c_str());
   if (!infile.good()) {
     MACH3LOG_ERROR("File {} does not exist", filename);
     MACH3LOG_ERROR("Please try again");
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   return filename;
 }

◆ 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 511 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;
 }

◆ MakeFluctuatedHistogramAlternative() [1/3]

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 469 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/3]

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

Make Poisson fluctuation of TH2D hist.

Definition at line 491 of file HistogramUtils.cpp.

                                                                                          {
 // ****************
   FluctHist->Reset();
   FluctHist->Fill(0.0, 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
  
   double x, y;
   for (int count = 0; count < num; ++count) {
     PolyHist->GetRandom2(x, y);
     FluctHist->Fill(x, y);
   }
 }

◆ MakeFluctuatedHistogramAlternative() [3/3]

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 540 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/3]

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

Make Poisson fluctuation of TH1D hist using default fast method.

Definition at line 408 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/3]

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

Make Poisson fluctuation of TH2D hist using default fast method.

Definition at line 446 of file HistogramUtils.cpp.

                                                                                   {
 // ****************
   // Reset the histogram
   FluctHist->Reset("");
   FluctHist->Fill(0.0, 0.0, 0.0);
  
   // Loop over all bins
   const int nBinsX = Hist->GetXaxis()->GetNbins();
   const int nBinsY = Hist->GetYaxis()->GetNbins();
   for (int ix = 1; ix <= nBinsX; ++ix) {
     for (int iy = 1; iy <= nBinsY; ++iy) {
       // Get the original bin content
       const double MeanContent = Hist->GetBinContent(ix, iy);
       // Generate Poisson fluctuation
       const double Random = rand->PoissonD(MeanContent);
       // Set the Random content
       FluctHist->SetBinContent(ix, iy, Random);
     }
   }
 }

◆ MakeFluctuatedHistogramStandard() [3/3]

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

Make Poisson fluctuation of TH2Poly hist using default fast method.

Definition at line 427 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 title of output histogram
BinArray_x	Bin edges for X axis
BinArray_y	Bin edges for Y axis

Definition at line 373 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;
 }

◆ MakeSummaryFromSpectra()

std::unique_ptr<TH1D> MakeSummaryFromSpectra	(	const TH2D *	Spectra,
		const std::string &	name
	)

Build a 1D posterior-predictive summary from a violin spectrum.

Parameters

Spectra	Input TH2D violin/density histogram.
name	Name of the output TH1D.

Definition at line 624 of file HistogramUtils.cpp.

                                                                     {
 // ***********************************************
   const int nBinsX = Spectra->GetNbinsX();
   // Reuse the exact x binning
   std::vector<double> edges(nBinsX + 1);
   for (int i = 0; i < nBinsX; ++i) {
     edges[i] = Spectra->GetXaxis()->GetBinLowEdge(i+1);
   }
   edges[nBinsX] = Spectra->GetXaxis()->GetBinUpEdge(nBinsX);
  
   auto h1 = std::make_unique<TH1D>(name.c_str(), name.c_str(), nBinsX, edges.data());
   h1->SetDirectory(nullptr);
   h1->Sumw2(true);
  
   // ---- Loop X bins and extract Y distribution
   for (int ix = 1; ix <= nBinsX; ++ix) {
     // Project THIS x-bin onto Y
     std::unique_ptr<TH1D> slice(Spectra->ProjectionY(Form("_py_%d", ix), ix, ix));
     slice->SetDirectory(nullptr);
     if (slice->GetEntries() == 0) {
       h1->SetBinContent(ix, 0.0);
       h1->SetBinError(ix, 0.0);
       continue;
     }
  
     const double mean = slice->GetMean();
     const double rms  = slice->GetRMS();
  
     h1->SetBinContent(ix, mean);
     h1->SetBinError(ix, rms);
   }
  
   h1->GetXaxis()->SetTitle(Spectra->GetXaxis()->GetTitle());
   h1->GetYaxis()->SetTitle("Events");
  
   return h1;
 }

◆ NoOverflowIntegral()

double NoOverflowIntegral ( TH2Poly * poly )

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

Definition at line 53 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 201 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 214 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 38 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 354 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 = `false`
	)

WP: Poly Projectors.

Definition at line 66 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 = `false`
	)

WP: Poly Projectors.

Definition at line 134 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 334 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 226 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 *	NumPoly,
		TH2Poly *	DenomPoly
	)

Helper to make ratio of TH2Polys.

Definition at line 238 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 397 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 ( const int PDG )

DB Get the Cherenkov momentum threshold in MeV.

Parameters

PDG	PDG code of the particle for which the Cherenkov threshold is requested.

Definition at line 574 of file HistogramUtils.cpp.

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

Namespaces

Functions

Detailed Description

Function Documentation

◆ CalculateEnu()

◆ CalculateQ2()

◆ CastVector()

◆ CheckTH2PolyFileVersion()

◆ ConvertTH2DtoTH2Poly()

◆ ConvertTH2PolyToTH2D()

◆ FastViolinFill()

◆ file_exists()

◆ GetRandomPoly2()

◆ MakeFluctuatedHistogramAlternative() [1/3]

◆ MakeFluctuatedHistogramAlternative() [2/3]

◆ MakeFluctuatedHistogramAlternative() [3/3]

◆ MakeFluctuatedHistogramStandard() [1/3]

◆ MakeFluctuatedHistogramStandard() [2/3]

◆ MakeFluctuatedHistogramStandard() [3/3]

◆ MakePolyHist()

◆ MakeSummaryFromSpectra()

◆ NoOverflowIntegral()

◆ NormalisePoly()

◆ NormaliseTH2Poly()

◆ OverflowIntegral()

◆ PolyIntegralWidth()

◆ PolyProjectionX()

◆ PolyProjectionY()

◆ PolyScaleWidth()

◆ RatioHists()

◆ RatioPolys()

◆ RemoveFitter()

◆ returnCherenkovThresholdMomentum()