plot_quantile_boundaries_twosamples_2020.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 "3FlavorAna/Cuts/QuantileCuts.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "Utilities/rootlogon.C"

#include "TCanvas.h"
#include "TColor.h"
#include "TFile.h"
#include "TH2.h"
#include "TLegend.h"
#include "TStyle.h"

#include <iostream>

using namespace ana;

void plot_quantile_boundaries_twosamples_2020(bool isFHC, 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(isFHC){
    horn_one = "fhc"; name_one = "FHC";
    horn_two = "rhc"; name_two = "RHC";
    color_one = kBlue-7; color_two = kMagenta+1;
  } else {
    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 = (std::string)std::getenv("NUMUDATA_DIR")+"lib/ana2020/Quantiles/quantiles_"+horn_one+"_full_numu2020.root";
  std::string filename2 = (std::string)std::getenv("NUMUDATA_DIR")+"lib/ana2020/Quantiles/quantiles_"+horn_two+"_full_numu2020.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->SetTitle("");
    quanthist1->GetYaxis()->SetTitle("E_{had} / E_{#nu}");
    quanthist1->GetXaxis()->CenterTitle();
    quanthist1->GetYaxis()->CenterTitle();
  }
  else{
    axes->SetTitle("");
    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][");
  }
  if (colzbkg) {
    if (isFHC ) CornerLabel( "Neutrino beam"     );
    else        CornerLabel( "Antineutrino beam" );
  }
  Simulation();
  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 = "quantiles_" + horn_one + "_" + horn_two + "_full_numu2020" + name_colz;
  gPad-> Print(outname + ".pdf");
  gPad-> Print(outname + ".png");
}