plot_quantbound_twosamples.C

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////////////
// plot_quantbound_twosamples.C                                         //
//                                                                      //
// print the quantile boundaries from two samples in the same canvas    //
// fhc = true prints fhc and rhc boundaries with fhc colz background    //
// uses rhc colz if fhc = false                                         //
//                                                                      //
// colzbkg = false prints the two boundaries but no colz background     //
//////////////////////////////////////////////////////////////////////////

#include "includes.h"

using namespace ana;


void plot_quantbound_twosamples(bool fhc = true, bool colzbkg = true){

  std::string bkg;
  std::string horn_one, horn_two;
  std::string name_one, name_two;
  Color_t color_one, color_two;
  Style_t line_one = kSolid, line_two = kSolid;
  int color_pal = kFuchsia;
  //TColor::InvertPalette();

  if(fhc){
    horn_one = "fhc"; name_one = "FHC";
    horn_two = "rhc"; name_two = "RHC";
    color_one = kBlue-7; color_two = kMagenta+1;
  }
  if(!fhc){
    horn_one = "rhc"; name_one = "RHC";
    horn_two = "fhc"; name_two = "FHC";
    color_one = kMagenta+1; color_two = kBlue-7;
  }


  const int nBins= 19;
  const int nQuant= 4;  
  std::string filename1 = "/nova/ana/nu_mu_ana/Ana2018/Quantiles/quantiles__" + horn_one + "_full__numu2018.root";
  std::string filename2 = "/nova/ana/nu_mu_ana/Ana2018/Quantiles/quantiles__" + horn_two + "_full__numu2018.root";

  TFile* infile1= TFile::Open(pnfs2xrootd(filename1).c_str());
  TFile* infile2= TFile::Open(pnfs2xrootd(filename2).c_str());
  TH2* quanthist1 = (TH2*)infile1->Get("FDSpec2D");
  TH2* quanthist2 = (TH2*)infile2->Get("FDSpec2D");
  std::vector< std::vector<double> > QuantBounds1 = GetQuantileBins(quanthist1, kNumuCCOptimisedAxis, kHadEFracAxis, nQuant, 1);
  std::vector< std::vector<double> > QuantBounds2 = GetQuantileBins(quanthist2, kNumuCCOptimisedAxis, kHadEFracAxis, nQuant, 1);

  std::vector< TH1* > hist_bound1;
  std::vector< TH1* > hist_bound2;

  // fill histograms in different loops, otherwise script messes up
  for(int thisQuant = 0; thisQuant < nQuant; thisQuant++){
    TH1* h_temp1 = quanthist1->ProjectionX();
    h_temp1-> Clear();
    for(int bin = 1; bin <= nBins; bin++){
      double bound1 = QuantBounds1[bin][thisQuant];
      h_temp1-> SetBinContent(bin, bound1);
    }
    hist_bound1.push_back(h_temp1);
  }
  for(int thisQuant = 0; thisQuant < nQuant; thisQuant++){
    TH1* h_temp2 = quanthist2->ProjectionX();
    h_temp2-> Clear();
    for(int bin = 1; bin <= nBins; bin++){
      double bound2 = QuantBounds2[bin][thisQuant];
      h_temp2-> SetBinContent(bin, bound2);
    }
    hist_bound2.push_back(h_temp2);
  }


  gStyle->SetPalette(color_pal);
  gStyle->SetPadRightMargin(0.11);
  new TCanvas("canvas");
  gPad->SetFillStyle(0);
  TH2* axes = new TH2F("axes", "", 20, 0.0, 5.0, 10, 0., 1.);
  std::string name_colz;
  if(colzbkg){
    name_colz = "__colzbkg";
    quanthist1-> Draw("colz");
    quanthist1->GetYaxis()->SetTitle("E_{had} / E_{#nu}");
    quanthist1->GetXaxis()->CenterTitle();
    quanthist1->GetYaxis()->CenterTitle();
  }
  else{
    axes->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
    axes->GetYaxis()->SetTitle("E_{had} / E_{#nu}");
    axes->GetXaxis()->CenterTitle();
    axes->GetYaxis()->CenterTitle();
    axes->Draw();
  }   
  for(int thisQuant = 1; thisQuant < nQuant; thisQuant++){
    hist_bound1[thisQuant] -> SetLineColor(color_one);
    hist_bound2[thisQuant] -> SetLineColor(color_two);
    hist_bound1[thisQuant] -> SetLineStyle(line_one);
    hist_bound2[thisQuant] -> SetLineStyle(line_two);
    hist_bound1[thisQuant] -> Draw("hist same][");
    hist_bound2[thisQuant] -> Draw("hist same][");
  }

  TLegend *legend = new TLegend(0.15,0.70,0.30,0.85);
  legend->AddEntry(hist_bound1[1], Form("%s", name_one.c_str()),"l");
  legend->AddEntry(hist_bound2[1], Form("%s", name_two.c_str()),"l");
  legend->SetTextSize(0.05); //no border for legend                                                                                                                                
  legend->SetBorderSize(0); //no border for legend                                                                                                                                 
  legend->SetFillStyle(0);  //fill colour is transparent                                                                                                                           
  legend->Draw();
  gPad->Update();


  TString outname = "ndfd_quantiles__" + horn_one + "_" + horn_two + "_full__numu2018" + name_colz;
  gPad-> Print(outname + ".pdf");
  gPad-> Print(outname + ".png");

}