HistogramUtils.h File Reference

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

Include dependency graph for HistogramUtils.h:

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

Go to the source code of 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, const std::string &TempName, const std::vector< double > &xbins, const bool computeErrors=false)
	WP: Poly Projectors. More...
TH1D *	PolyProjectionY (TObject *poly, const 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)
	Make Poisson fluctuation of TH1D hist using default fast method. More...
void	MakeFluctuatedHistogramAlternative (TH1D *FluctHist, TH1D *PolyHist)
	Make Poisson fluctuation of TH1D hist using slow method which is only for cross-check. More...
void	MakeFluctuatedHistogramStandard (TH2D *FluctHist, TH2D *PolyHist)
	Make Poisson fluctuation of TH2D hist using default fast method. More...
void	MakeFluctuatedHistogramAlternative (TH2D *FluctHist, TH2D *PolyHist)
	Make Poisson fluctuation of TH2D hist. More...
void	MakeFluctuatedHistogramStandard (TH2Poly *FluctHist, TH2Poly *PolyHist)
	Make Poisson fluctuation of TH2Poly hist using default fast method. More...
void	MakeFluctuatedHistogramAlternative (TH2Poly *FluctHist, TH2Poly *PolyHist)
	Make Poisson fluctuation of TH2Poly hist using slow method which is only for cross-check. More...
int	GetRandomPoly2 (const TH2Poly *PolyHist)
	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...
std::vector< double >	BuildBinEdgesFromNode (YAML::Node const &bin_edges_node, bool &found_range_specifier)
	Builds a single dimension's bin edges from YAML::Node. More...
std::vector< double >	BinRangeToBinEdges (YAML::Node const &bin_range)
	Converts a range (linspace/logspace) to a std::vector<double> 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

BinRangeToBinEdges()

std::vector<double> BinRangeToBinEdges

(

YAML::Node const &

bin_range

)

Converts a range (linspace/logspace) to a std::vector<double>

Parameters

bin_range

Yaml node containing binning

Definition at line 653 of file HistogramUtils.cpp.

                                                                {
// ************************************************
  bool is_lin = true;
  YAML::Node bin_range_specifier;
  if (bin_range["linspace"]) {
    bin_range_specifier = bin_range["linspace"];
  } else if (bin_range["logspace"]) {
    is_lin = false;
    bin_range_specifier = bin_range["logspace"];
  } else {
    std::stringstream ss;
    ss << bin_range;
    MACH3LOG_ERROR("When parsing binning, expected bin range specifier with "
                   "key linspace or logspace, but found,\n{}",
                   ss.str());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  auto nb = Get<int>(bin_range_specifier["nb"], __FILE__, __LINE__);
  auto low = Get<double>(bin_range_specifier["low"], __FILE__, __LINE__);
  auto up = Get<double>(bin_range_specifier["up"], __FILE__, __LINE__);
 
  std::vector<double> edges(nb + 1, low);
  // force the last bin to be exactly as parsed to avoid numerical instabilities
  // not quite lining back up with the end of the range, which could
  // cause spurious errors or infinitesimally small bins for specifications
  // like: [ { logspace: { nb: 10, 1E-1, 10}, 10, 11} ]
  edges.back() = up;
 
  if (is_lin) {
    double bw = (up - low) / nb;
    for (int i = 0; i < (nb - 1); ++i) {
      edges[i + 1] = edges[i] + bw;
    }
  } else {
    double llow = std::log10(low);
    double lup = std::log10(up);
    double lbw = (lup - llow) / nb;
    for (int i = 0; i < (nb - 1); ++i) {
      edges[i + 1] = std::pow(10, llow + (i + 1) * lbw);
    }
  }
 
  return edges;
}

BuildBinEdgesFromNode()

std::vector<double> BuildBinEdgesFromNode	(	YAML::Node const &	bin_edges_node,
		bool &	found_range_specifier
	)

Builds a single dimension's bin edges from YAML::Node.

Parameters

bin_edges_node	Yaml node containing binning
found_range_specifier	Whether or not a range specifier (linspace/logspace) was found

BinEdges: [ <dim0bin0lowedge>, <dim0bin1upedge>, <dim0bin2upedge>, ... <dim0binNupedge> ] BinEdges: { linspace: { nb: 100, low: 0, up: 10} } BinEdges: { logspace: { nb: 100, low: 1E-1, up: 10} } BinEdges: [ { linspace: { nb: 100, low: 0, up: 10} }, 10, 15, { logspace: { nb: 5, low: 15, up: 100} } ]

Definition at line 707 of file HistogramUtils.cpp.

                                                        {
// ************************************************
  if (bin_edges_node.IsMap()) {
    found_range_specifier = true;
    return BinRangeToBinEdges(bin_edges_node);
  }
  std::vector<double> edges_builder;
  if (!bin_edges_node.IsSequence()) {
    std::stringstream ss;
    ss << bin_edges_node;
    MACH3LOG_ERROR(
        "When parsing binning, expected to find a YAML map or sequence, "
        "but found:\n{}",
        ss.str());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  for (auto const &it : bin_edges_node) {
    if (it.IsScalar()) {
      edges_builder.push_back(it.as<double>());
    } else if (it.IsMap()) {
      found_range_specifier = true;
      auto range_edges = BinRangeToBinEdges(it);
      std::copy(range_edges.begin(), range_edges.end(),
                std::back_inserter(edges_builder));
    } else {
      std::stringstream ss;
      ss << bin_edges_node;
      MACH3LOG_ERROR(
          "When parsing binning, expected elements in outer sequence to all be "
          "either scalars or maps, but found:\n{}",
          ss.str());
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
 
  // Check for duplicates or out-of-order bins
  std::vector<double> edges;
  for (size_t eb_it = 0; eb_it < edges_builder.size(); ++eb_it) {
    if (edges.size()) {
      if (edges_builder[eb_it] == edges.back()) { // remove duplicate edges
        continue;
      } else if (edges_builder[eb_it] < edges.back()) {
        std::stringstream ss;
        ss << "[ ";
        for(auto const & e : edges_builder){
          ss << fmt::format("{:.3g} ", e);
        }
        ss << "]";
        MACH3LOG_ERROR(
            "When parsing binning, found edges that were not monotonically "
            "increasing, problem bin at index: {}:\n{}",
            eb_it, ss.str());
        throw MaCh3Exception(__FILE__, __LINE__);
      }
    }
    edges.push_back(edges_builder[eb_it]);
  }
 
  return edges;
}

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 593 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 573 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 94 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 13 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 284 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 250 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 551 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 104 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

)

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

Definition at line 504 of file HistogramUtils.cpp.

                                            {
// ****************
  const int nbins = PolyHist->GetNumberOfBins();
  const double r1 = M3::rand::Uniform();;
 
  std::vector<double> fIntegral(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.data(), r1));
  //KS: Have to increment because TH2Poly has stupid offset arghh
  iBin += 1;
 
  return iBin;
}

MakeFluctuatedHistogramAlternative() [1/3]

void MakeFluctuatedHistogramAlternative	(	TH1D *	FluctHist,
		TH1D *	PolyHist
	)

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

Definition at line 468 of file HistogramUtils.cpp.

                                                                        {
// ****************
  // Make the Poisson fluctuated hist
  FluctHist->Reset("");
  FluctHist->Fill(0.0, 0.0);
 
  const double evrate = PolyHist->Integral();
  const int num = M3::rand::Poisson(evrate);
  int count = 0;
  while(count < num)
  {
    const double candidate = PolyHist->GetRandom();
    FluctHist->Fill(candidate);
    count++;
  }
}

MakeFluctuatedHistogramAlternative() [2/3]

void MakeFluctuatedHistogramAlternative	(	TH2D *	FluctHist,
		TH2D *	PolyHist
	)

Make Poisson fluctuation of TH2D hist.

Definition at line 487 of file HistogramUtils.cpp.

                                                                         {
// ****************
  FluctHist->Reset();
  FluctHist->Fill(0.0, 0.0, 0.0);
 
  const double evrate = PolyHist->Integral();
  const int num = M3::rand::Poisson(evrate);
 
  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
	)

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

Definition at line 532 of file HistogramUtils.cpp.

                                                                              {
// ****************
  // Make the Poisson fluctuated hist
  FluctHist->Reset("");
  FluctHist->Fill(0.0, 0.0, 0.0);
 
  const double evrate = NoOverflowIntegral(PolyHist);
  const int num = M3::rand::Poisson(evrate);
  int count = 0;
  while(count < num)
  {
    const int iBin = GetRandomPoly2(PolyHist);
    FluctHist->SetBinContent(iBin, FluctHist->GetBinContent(iBin) + 1);
    count++;
  }
}

MakeFluctuatedHistogramStandard() [1/3]

void MakeFluctuatedHistogramStandard	(	TH1D *	FluctHist,
		TH1D *	PolyHist
	)

Make Poisson fluctuation of TH1D hist using default fast method.

Definition at line 407 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 = M3::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
	)

Make Poisson fluctuation of TH2D hist using default fast method.

Definition at line 445 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 = M3::rand::PoissonD(MeanContent);
      // Set the Random content
      FluctHist->SetBinContent(ix, iy, Random);
    }
  }
}

MakeFluctuatedHistogramStandard() [3/3]

void MakeFluctuatedHistogramStandard	(	TH2Poly *	FluctHist,
		TH2Poly *	PolyHist
	)

Make Poisson fluctuation of TH2Poly hist using default fast method.

Definition at line 426 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 = M3::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 372 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 613 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 52 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 200 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 213 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 37 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 353 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,
		const std::string &	TempName,
		const std::vector< double > &	xbins,
		const bool	computeErrors = false
	)

WP: Poly Projectors.

Definition at line 65 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,
		const std::string &	TempName,
		const std::vector< double > &	ybins,
		const bool	computeErrors = false
	)

WP: Poly Projectors.

Definition at line 133 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 333 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 225 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 237 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 396 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 563 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;
}