Plots.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Core/FwdDeclare.h"
#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Utilities.h"

#include "CAFAna/Prediction/IPrediction.h"

#include <map>
#include <vector>

#include "Rtypes.h"

class TCanvas;
class TGraph;
class TGraphAsymmErrors;
class TH1;
class TH2;
class THStack;
class TLegend;
class TMarker;
class TPad;
class TPaveText;

namespace ana
{
  class IPrediction;
  class SystShifts;

  /// Overlay MC spectrum with data spectrum
  ///
  /// \return The first histogram drawn so you can alter axis labels etc
  TH1* DataMCComparison(const Spectrum& data, const Spectrum& mc, EBinType bintype = kBinContent);


  /// Overlay MC spectrum with data spectrum
  ///
  /// \return The first histogram drawn so you can alter axis labels etc
  TH1* DataMCComparison(const Spectrum& data,
                        const IPrediction* mc,
                        osc::IOscCalc* calc,
                        const SystShifts & shifts = kNoShift,
                        EBinType bintype = kBinContent);

  /// Overlay MC spectrum with data spectrum, area normalized
  ///
  /// \return The first histogram drawn so you can alter axis labels etc
  TH1* DataMCComparisonAreaNormalized(const Spectrum& data, const Spectrum& mc);

  /// Overlay MC spectrum with data spectrum, area normalized
  ///
  /// \return The first histogram drawn so you can alter axis labels etc
  TH1* DataMCComparisonAreaNormalized(const Spectrum& data,
                   const IPrediction* mc,
                   osc::IOscCalc* calc);


  /// \brief Plot MC broken down into flavour components, overlayed with data
  ///
  /// \return The first histogram drawn so you can alter axis labels etc
  TH1* DataMCComparisonComponents(const Spectrum& data,
                                  const IPrediction* mc,
                                  osc::IOscCalc* calc);

  /// Plot data/MC ratio for the given spectrum. Normalize MC to Data by area
  void DataMCAreaNormalizedRatio(const Spectrum& data, const Spectrum& mc,
                                 double miny = 0, double maxy = 3);

  /// Plot data/MC ratio for the given spectrum. Normalize MC to Data by area
  void DataMCAreaNormalizedRatio(const Spectrum& data,
                                 const IPrediction* mc,
                                 osc::IOscCalc* calc,
                                 double miny = 0, double maxy = 3);

  /// Plot data/MC ratio for the given spectrum. Normalize MC to Data by POT
  void DataMCRatio(const Spectrum& data, const Spectrum& mc,
                   double miny = 0, double maxy = 3);

  /// Plot data/MC ratio for the given spectrum. Normalize MC to Data by POT
  void DataMCRatio(const Spectrum& data,
                   const IPrediction* mc,
                   osc::IOscCalc* calc,
                   double miny = 0, double maxy = 3);

  /// Plot data/expected, compared with fit/expected
  void RatioPlot(const Spectrum& data,
                 const Spectrum& expected,
                 const Spectrum& fit,
                 double miny = 0, double maxy = 1.2);

  /// Return TH1 vectors of the +/-1sigma syst shifted prediction
  ///
  /// \param calc  Oscillation parameters
  /// \param pred  The prediction. Must be able to generate shifted spectra
  /// \param systs List of systematics to include in error band
  /// \param hUps  Empty TH1* vector to return the +1sigma shifted hists
  /// \param hDns  Empty TH1* vector to return the -1sigma shifted hists
  /// \param pot   POT to scale the prediction to
  void GetSystBands(osc::IOscCalc *calc,
                    const IPrediction* pred,
                    std::vector<const ISyst*> systs,
                    std::vector<TH1*> &hUps,
                    std::vector<TH1*> &hDns,
                    double pot,
                    Flavors::Flavors_t flav = Flavors::kAll,
                    Current::Current_t curr = Current::kBoth,
                    Sign::Sign_t sign = Sign::kBoth);

  /// Return TH1 vectors of the +/-1sigma syst shifted prediction
  /// with a given set of best fit systematic pulls 
  ///
  /// \param calc     Oscillation parameters
  /// \param pred     The prediction. Must be able to generate shifted spectra
  /// \param systs    List of systematics to include in error band
  /// \param bfshifts Systematic shifts determined by the best fit
  /// \param hUps     Empty TH1* vector to return the +1sigma shifted hists
  /// \param hDns     Empty TH1* vector to return the -1sigma shifted hists
  /// \param pot   POT to scale the prediction to
  void GetBFSystBands(osc::IOscCalc *calc,
                      const IPrediction* pred,
                      std::vector<const ISyst*> systs,
                      const SystShifts &bfshifts,
                      std::vector<TH1*> &hUps,
                      std::vector<TH1*> &hDns,
                      double pot,
                      Flavors::Flavors_t flav = Flavors::kAll,
                      Current::Current_t curr = Current::kBoth,
                      Sign::Sign_t sign = Sign::kBoth);

    
  /// \brief Plot prediction with +/-1sigma error band
  ///
  /// When multiple systematics are used, the errors are the quadrature sum
  ///
  /// \param pred   The prediction. Must be able to generate shifted spectra
  /// \param systs  List of systematics to include in error band
  /// \param calc   Oscillation parameters
  /// \param pot    POT to evaluate prediction at
  /// \param col    Color of the prediction, default kRed
  /// \param errCol Color of the shading, default col-7 (kRed-7 is light red)
  /// \param headroom   Fraction of maximum bin for headroom, default 30%
  /// \param alpha  Control transparency of error bands, default is no transparency
  TGraphAsymmErrors*  PlotWithSystErrorBand(IPrediction* pred,
                                            const std::vector<const ISyst*>& systs,
                                            osc::IOscCalc* calc,
                                            double pot,
                                            int col = -1, int errCol = -1,
                                            float headroom = 1.3,
                                            bool newaxis = true,
                                            EBinType bintype = kBinContent,
              double alpha=1);

  /// \brief Plot prediction with error band
  ///
  /// When multiple systematics are used, the errors are the quadrature sum
  ///
  /// \param nominal    Nominal spectrum
  /// \param upShifts   Vector of spectra which have + shifts
  /// \param downShifts Vector of spectra which have - shifts, same order as +
  /// \param pot        POT to scale spectra to
  /// \param col        Color of the prediction, default kRed
  /// \param errCol     Color of the shading, default col-7(kRed-7 is light red)
  /// \param headroom   Fraction of maximum bin for headroom, default 30%
  /// \param alpha      Control transparency of error bands, default is no transparency
  TGraphAsymmErrors*  PlotWithSystErrorBand(const Spectrum& nominal,
                                            const std::vector<Spectrum>& upShifts,
                                            const std::vector<Spectrum>& downShifts,
                                            double pot,
                                            int col = -1, int errCol = -1,
                                            float headroom = 1.3, bool newaxis=true,
                                            EBinType bintype = kBinContent,
              double alpha=1);

   TGraphAsymmErrors*  PlotWithSystErrorBand(TH1*& nominal,
                                             std::vector<TH1*>& upShifts,
                                             std::vector<TH1*>& downShifts,
                                             int col = -1, int errCol = -1,
                                             float headroom = 1.3,
                                             bool newaxis = true,
               double alpha=1);

   /// \brief Plot two different predictions with +/-1sigma shape-only error bands
   void PlotWithSystErrorBandTwoPreds(const Spectrum& nominal,
                                      const std::vector<Spectrum>& upShifts,
                                      const std::vector<Spectrum>& downShifts,
                                      const Spectrum& nominal2,
                                      const std::vector<Spectrum>& upShifts2,
                                      const std::vector<Spectrum>& downShifts2,
                                      double pot,
                                      int col = -1, int errCol = -1,
                                      float headroom = 1.3, bool newaxis=true,
                                      EBinType bintype = kBinContent);



   /// Equivalent default PlotWithSystErrorBand but adding the option for numu quartiles
   ///
   /// All quartile pads should idealy share the same maximum for easier complarison
   /// TGraph is usually bad behaved when plotted with histograms, so have to force a
   /// maximum from here: changing float headroom to float maxY
   TGraphAsymmErrors*  PlotWithSystErrorBand_Quant(const int quant,
                                                   IPrediction* pred,
                                                   const std::vector<const ISyst*>& systs,
                                                   osc::IOscCalc* calc,
                                                   double pot,
                                                   int col = -1, int errCol = -1,
                                                   float maxy = 8.0,
                                                   bool newaxis = true);


   TGraphAsymmErrors*  PlotWithSystErrorBand_Quant(const int quant,
                                                   const Spectrum& nominal,
                                                   const std::vector<Spectrum>& upShifts,
                                                   const std::vector<Spectrum>& downShifts,
                                                   double pot,
                                                   int col = -1, int errCol = -1,
                                                   float mayy = 8.0,
                                                   bool newaxis=true);


   TGraphAsymmErrors*  PlotWithSystErrorBand_Quant(const int quant,
                                                   TH1*& nominal,
                                                   std::vector<TH1*>& upShifts,
                                                   std::vector<TH1*>& downShifts,
                                                   int col = -1, int errCol = -1,
                                                   float maxy = 8.0,
                                                   bool newaxis = true);


  /// \brief Plot prediction with +/-1sigma shape-only error band
  ///
  /// When multiple systematics are used, the errors are the quadrature sum
  ///
  /// \param pred   The prediction. Must be able to generate shifted spectra
  /// \param systs  List of systematics to include in error band
  /// \param calc   Oscillation parameters
  /// \param pot    POT to evaluate prediction at
  /// \param col    Color of the prediction, default kRed
  /// \param errCol Color of the shading, default col-7 (kRed-7 is light red)
  /// \param headroom   Fraction of maximum bin for headroom, default 30%
  void PlotWithAreaSystErrorBand(IPrediction* pred,
                                 const std::vector<const ISyst*>& systs,
                                 osc::IOscCalc* calc,
                                 double pot,
                                 int col = -1, int errCol = -1,
                                 float headroom = 1.3,
                                 bool newaxis = true,
                                                   EBinType bintype = kBinContent);

  /// \brief Plot prediction with error band where syst(s) is (are) shape-only
  ///
  /// When multiple systematics are used, the errors are the quadrature sum.
  /// Systs are area normalized to nominal.
  ///
  /// \param nominal    Nominal spectrum
  /// \param upShifts   Vector of spectra which have + shifts
  /// \param downShifts Vector of spectra which have - shifts, same order as +
  /// \param pot        POT to scale spectra to
  /// \param col        Color of the prediction, default kRed
  /// \param errCol     Color of the shading, default col-7(kRed-7 is light red)
  /// \param headroom   Fraction of maximum bin for headroom, default 30%
  void PlotWithAreaSystErrorBand(const Spectrum& nominal,
                                 std::vector<Spectrum> upShifts,   // yup, copy them.
                                 std::vector<Spectrum> downShifts, // we need to change their normalizations
                                 double pot,
                                 int col = -1, int errCol = -1,
                                 float headroom = 1.3, bool newaxis=true,
                                                   EBinType bintype = kBinContent);


  /// Can call like ToTHStack({{h1, kRed}, {h2, kBlue}}, pot)
  THStack* ToTHStack(const std::vector<std::pair<const Spectrum&, Color_t>>& s,
                     double pot);


  /// \brief Create a legend, maximizing distance from all histograms
  ///
  /// \param dx Width of the legend, fraction of pad width
  /// \param dy Height of the legend, fraction of pad height
  /// \param yPin, height in NDC (frac) to pin center of legend, i.e only move x
  TLegend* AutoPlaceLegend(double dx, double dy, double yPin = -1);

  /// \brief Make a simple plot of relative size of different errors
  ///
  /// \param systs     Map from error label to -, + size
  /// \param statErr   Percentage statistical error. Optional.
  /// \param label     Label you want on the axis.
  void ErrorBarChart(const std::map<std::string, std::pair<double, double>>& systErrors,
                     const std::pair<double, double> & statErr = std::make_pair(0., 0.),
                     const std::string & label = "Error" );

  /// \brief Make a simple plot of relative size of different errors
  ///
  /// \param systs     Map from error label to size in percent
  /// \param statErr   Percentage statistical error. Optional.
  /// \param bkgdOrSig Label axis for background or signal. Default (false) = bkgd
  void CountingExperimentErrorBarChart(const std::map<std::string, double>& systs, double statErr = 0, bool bkgdOrSig = false, bool shortchart=false);

  /// Calculate statistical errors appropriate for small Poisson numbers
  TGraphAsymmErrors* GraphWithPoissonErrors(const TH1* h, bool noErrorsXaxis = false, bool drawEmptyBins = true);

  /// Same as above but use a reference histogram to determine which empty bins to draw
  TGraphAsymmErrors* GraphWithPoissonErrors2(const TH1* h, const TH1* h2, bool noErrorsXaxis = false, bool drawEmptyBins = true);

  /// Returns graph with the right size of error for a scaled histogram
  /// Mosty for use in numu reconstructed energy spectra when scaled by 0.1 GeV and binwidth
  /// But scaleFactor can be of course modified locally
  /// Takes a the histogram already saled and the original
  TGraph* graphAsymmErrorScaled(TH1* histScaled, TH1* hist, double overallScale = 0.1);

  /// Returns graph with the right size of error for a scaled histogram
  /// Mosty for use in numu reconstructed energy spectra when scaled by 0.1 GeV and binwidth
  /// But scaleFactor can be of course modified locally
  /// This implementation is slightly different than graphAsymmErrorScaled as here the 
  /// scaling is manually reversed so there is no need for unscaled hists to be included.
  /// Takes scaled hists of data and estimated background.
  TGraph* graphAsymmErrorWithBkgScaled(TH1* data, TH1* bkgd, double overallScale = 0.1);
      
  /// Returns graph with the right size of error for a ratio histogram
  /// Hists need to be unscaled
  TGraph* RatioAsymmError(TH1* data, TH1* pred);

  /// Returns graph with the right size of error for a ratio histogram
  /// For histograms that have been scaled by bin width. This implementation
  /// is slightly different than graphAsymmErrorScaled as here the scaling is
  /// manually reversed so there is no need for unscaled hists to be included
  TGraph* RatioAsymmErrorScaled(TH1* data, TH1* pred, double overallScale = 0.1);

  /// Returns graph with the right size of error for a ratio histogram
  /// where you have a measure of the background. Hists need to be unscaled
  TGraph* RatioAsymmErrorWithBkg(TH1* data, TH1* pred, TH1* bkg);

  /// Returns graph with the right size of error for a ratio histogram
  /// where you have a measure of the background. Hists need to be unscaled
  TGraph* RatioAsymmErrorWithBkgScaled(TH1* data, TH1* pred, TH1* bkg, double overallScale = 0.1);

  /// Gives a TGraph with the area between two histograms. Do Draw("f") to draw
  /// this area. By default it has a lighter version of the colour of hmin
  TGraph* ShadeBetweenHistograms(TH1* hmin, TH1* hmax);

  /// Calculate profile with error bars corresponding to specified quantiles of a 2D distribution (by default, 68% coverage)
  TGraphAsymmErrors * ProfileQuantile(const TH2 * hist, const std::string & axis_name, const std::string & graph_name="", const std::pair<double, double> & quantile_divisions={0.159, 0.841});


  /// Draw a single BF point
  void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size);

  /// Draw best fit at both octants. Truth for neutrinos, not sure about antineutrinos
  void drawBFMirror(double bfSin, double bfDm, Color_t color, Style_t marker, double size);

  /// Make histogram be ready for use in SplitCanvas_Quant
  /// Depending the quantile, this sets the axes size to zero so they dont overlay with each other
  void MakeHistCanvasReady_Quant(const int quant, TH1* hist, double ymax);


  /// Pimp histogram once and for all
  void PimpHist(TH1* hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize);


  /// \brief Split the current pad into two vertically stacked pieces, suitable
  /// for ratio plots and so on.
  ///
  /// This is better than TCanvas::Divide() in that it keeps axis label font
  /// sizes as expected, and doesn't truncate the zero on the top y-axis. It's
  /// worse in that you won't be able to interact graphically with the bottom
  /// pad.
  ///
  /// \param  ysplit The y position (between zero and one) where the current
  /// pad will be split.
  /// \param[out] p1 Returns the top canvas created by the function
  /// \param[out] p2 Returns the bottom canvas created by the function
  void SplitCanvas(double ysplit, TPad*& p1, TPad*& p2);

  /// Split canvas in four pads for easy drawing of 4 quantiles
  /// TO DO: Educate ROOT so it doesnt overlay the x-axes
  void SplitCanvasQuant(TCanvas *& canvas, TPad *& pad1, TPad *& pad2, TPad *& pad3, TPad *& pad4);

  void CenterTitles(TH1* histo);

  /// Create a simple summary of systematic shifts.
  /// Systematic short names in the X axis, and n sigma in the y axis.
  /// Paints a TCanvas, and returns a pointer to the TH1 for further editing
  TH1* PlotSystShifts(const SystShifts & shifts, bool sortName = true);

  /// Join graphs and set a fill color. Useful for contours
  TGraph* JoinGraphs(TGraph* a, TGraph* b, int fillcol = 0);

  /// Add two points to a TGraph so it reaches the top, and set a fill color.
  /// Default is for a nue dCP x ssth23 contour
  TGraph* ExtendGraphToTop(TGraph* g, int col = 0, double xmin = 0,
                           double xmax = 2, double y = 1);

  /// Put the standardized beam label in the left corner of the active canvas
  TPaveText* DrawBeamLabel(bool isFHC = true);

  /// Put hadEQuant labels on the active canvas or pad
  /// quant = 0 is All Quartiles
  TPaveText* DrawQuantLabel(int quant);

  
}