plot_nd_spectra_2018.C

Go to the documentation of this file.

#pragma once

#include "CAFAna/Analysis/Exposures.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "Utilities/rootlogon.C"

#include "TFile.h"
#include "TH1.h"
#include "THStack.h"
#include "TCanvas.h"
#include "TGaxis.h"
#include "TGraph.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TLine.h"
#include "TSystem.h"
#include "TPaveLabel.h"
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <cstddef> 

using namespace ana;

void SetAxis(TH1* h, bool visible, int color, int style = 1, bool scale = true);
TString POTstr(double pot);
double DerivePowScale(TH1* h);
TH1* MakeRatio(TH1* num, TH1* denom, int color=1);
void One();
void Zero();
void Legend(TH1* h_to_clone);
void RatioLegend();
void NeutrinoMode(std::string mode);


void plot_nd_spectra_2018(std::string mode = "fhc") {
  
  std::string fdir = "/nova/ana/nu_e_ana/Ana2018/DataMC/";

  TFile* fdata = new TFile((fdir+"20180331_Ana2018_"+mode+"_datamc_decomp_data_validation.root").c_str(), "READ");
  TFile *fmc = new TFile((fdir+ "20180331_Ana2018_"+mode+"_datamc_decomp_mc_validation.root").c_str(), "READ");
  TFile* fdecomp = new TFile((fdir + "20180331_Ana2018_"+mode+"_datamc_decomp_mc_decomp_validation.root").c_str(), "READ");

  if(fmc->IsZombie()) {std::cout << "MC File not found "; return; }
  if(fdata->IsZombie()) {std::cout << "Data File not found "; return; }
  if(fdecomp->IsZombie()) {std::cout << "MC Decomp File not found"; return; }

  // Pulling Data        
  auto spec = Spectrum::LoadFrom(fdata, (TString)"Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}");
  double datapot = spec->POT();
  TH1* hdata = spec->ToTH1(datapot);
 
  std::vector <TH1*> mc;
  std::vector <TH1*> decomp;
  std::vector <std::string> comps = {"0_nue", "1_numu", "2_nc", "3_antinue", "4_antinumu"};
  TH1* tempmc;
  TH1* tempdecomp;
  unique_ptr<Spectrum> tempspec;
  // Load spectra from uncorrected mc and decomp corrected mc files
  for(std::string comp : comps ) {
    tempspec = Spectrum::LoadFrom(fmc, (TString::Format("Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat=%s}",comp.c_str())));
    tempmc = tempspec->ToTH1(datapot);
    
    std::cout << "Loaded MC Spectrum Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat="+comp+"}" << endl;; 
    
    tempspec = Spectrum::LoadFrom(fdecomp, (TString::Format("Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat=%s}",comp.c_str())));
    tempdecomp = tempspec->ToTH1(datapot);
    
    std::cout << "Loaded Decomp Spectrum Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat="+comp+"}" << endl;

    mc.push_back(tempmc);
    decomp.push_back(tempdecomp);
  }
  
  // Add up all of the components to get total MC
  TH1* totalmc = (TH1*)mc[0]->Clone("totalmc");
  TH1* totaldecomp = (TH1*)decomp[0]->Clone("totaldecomp");
  for(int i = 1; i < (int)mc.size(); i++) {
    totalmc->Add( mc[i] );
    totaldecomp->Add( decomp[i] );
    std::cout << "Adding component " + std::to_string(i) << endl;
  }
  mc.push_back(totalmc);
  decomp.push_back(totaldecomp);
  std::cout << "Added total MC and Decomp" << endl;

  int tot = mc.size()-1;
  TH1* rdecomp;
  TH1* rmc;
  // make the histograms for the ratios
  rdecomp = MakeRatio(hdata, decomp[tot]);
  rmc = MakeRatio(hdata, mc[tot], 17);


  // Lets start drawing
  // Split plot
  TCanvas* c = new TCanvas("c");
  double ysplit=0.35;
  TPad* p1 =new TPad("p1", "", 0,0,1,1);
  TPad* p2 =new TPad("p2","",0,0,1,1);
  p1->SetBottomMargin(ysplit);
  p2->SetTopMargin(1-ysplit);
  p1->SetFillStyle(0);
  p2->SetFillStyle(0);

  p1->cd();


  std::vector <int> colors = {kMagenta,
                              kGreen+2,
                              kAzure,
                              kMagenta-2,
                              kGreen-3,
                              kRed};
  std::vector<std::string> labels = {"#nu_{e} CC",
                                     "#nu_{#mu} CC",
                                     "NC",
                                     "#bar{#nu}_{e} CC",
                                     "#bar{#nu}_{#mu} CC"};
                                     
                

  decomp[tot]->GetYaxis()->SetTitle(POTstr(datapot));
  decomp[tot]->GetYaxis()->CenterTitle();

  SetAxis(decomp[tot], true, colors[tot]);
  double maxY = decomp[tot]->GetMaximum();
  decomp[tot]->GetYaxis()->SetRangeUser(0, 1.30*maxY);
  decomp[tot]->Draw("hist");
  Legend(decomp[tot]);

  // Makes it look fancy
  Nue2018ThreeBinDivisions(true);
  Nue2018ThreeBinLabels(1.42*maxY, 0.06, kGray+3, true);

  SetAxis(mc[tot], false, colors[tot], 7);
  mc[tot]->Draw("hist same");
  for(int j = 0; j < (int)mc.size()-1 ; j++) {
    SetAxis(decomp[j], false, colors[j]);
    decomp[j]->Draw("hist same");
    SetAxis(mc[j], false, colors[j], 7);
    mc[j]->Draw("hist same");
  }

  hdata->SetMarkerStyle(kFullCircle);
  SetAxis(hdata, false, 1);
  hdata->Draw("ep same");
  
  // Add fancy axes
  Nue2018ThreeBinAxis(decomp[tot], false, true, true); 
  decomp[tot]->GetXaxis()->SetLabelSize(0.);
  decomp[tot]->GetXaxis()->SetTitleSize(0.);
  decomp[tot]->GetYaxis()->SetTitleSize(0.75*decomp[tot]->GetYaxis()->GetTitleSize());
  decomp[tot]->GetYaxis()->SetRangeUser(0.001, 
                            decomp[tot]->GetMaximum());

  // The empty bins of the histograms still get drawn at y=0. Try to 
  // draw over these colored lines with black
  Zero();

  // Draw the ratio
  p2->cd();

  SetAxis(rdecomp, true, 1, 1, false);
  rdecomp->Draw("ep");

  SetAxis(rmc, false, 17, 1, false);
  rmc->Draw("ep same");
  One();
  Nue2018ThreeBinDivisions(true);
  Nue2018ThreeBinAxis(rdecomp, true, true, true);

  rdecomp->GetYaxis()->SetTitleSize(0.75*rdecomp->GetYaxis()->GetTitleSize());
  rmc->GetYaxis()->SetRangeUser(.5,1.5);
  rdecomp->GetYaxis()->SetRangeUser(.5,1.5);
  RatioLegend();
  Preliminary();
  NeutrinoMode(mode);

  p1->Draw();
  p2->Draw("same");
  
  gPad->Print("Energy_Ana2018_Bin_CVNeLow_split.pdf");
}

void SetAxis(TH1* h, bool visible, int color, int style, bool scale) {
  if(scale) h->Scale(.001);
  h->SetLineStyle(style);
  h->SetLineColor(color);
  if(!visible){
    h->SetTickLength(0.);
    h->SetLabelSize(0.);
  }
  double norm = DerivePowScale(h);
}

TString POTstr(double pot) {
  return TString(TString::Format("10^{3} Events / %.3g#times10^{20} POT", pot/1e20)); 
}

double DerivePowScale(TH1* h) {
  // "Next smallest power of 1000 than maximum bin" - from plot_validation_datamc.py
  return 3*(int(log10(h->GetMaximum())))/3;
}

TH1* MakeRatio(TH1* num, TH1* denom, int color) {
  TH1* ret = (TH1*)num->Clone();
  ret->Divide(denom);

  ret->SetMarkerStyle(kFullCircle);
  ret->SetMarkerColor(color);
  ret->SetLineColor(color);
  ret->GetYaxis()->SetRangeUser(.5, 1.5);
  ret->GetYaxis()->SetTitle("Data / MC");
  ret->GetXaxis()->SetNdivisions(1900, kFALSE);
  //  ret->GetXaxis()->SetTickLength(1.);
  return ret;
}

void Zero() {
  TGraph* g = new TGraph();
  g->SetPoint(0, -1e10,0.003);
  g->SetPoint(1, 1e10, 0.003);
  g->SetLineWidth(2);
  g->SetLineStyle(1);
  g->Draw("same");

}

void One(){
  TGraph* g = new TGraph();
  g->SetPoint(0, -1e10,1);
  g->SetPoint(1, 1e10, 1);
  g->SetLineWidth(2);
  g->SetLineStyle(7);
  g->Draw("same");
}

void Legend(TH1* h_to_clone){
  std::vector <std::pair<int, TString>> entries = {{kBlack, "ND data"},
                                                   {kRed, "Total MC"},
                                                   {kAzure, "NC"},
                                                   {kGreen+2, "#nu_{#mu} CC"},
                                                   {kMagenta, "#nu_{e} CC"},
                                                   {kGreen-3, "#bar{#nu}_{#mu} CC"},
                                                   {kMagenta-2, "#bar{#nu}_{e} CC"},
                                                   {kBlack, "Uncorr. MC"},
                                                   {17, "Uncorr. MC"},
                                                   {1, "Decomp MC"}};
                                                       

  TLegend* leg = new TLegend(.78, .4, .95, .89);

  leg->SetFillColor(kWhite);
  leg->SetBorderSize(1);
  leg->SetMargin(.22);

  TH1* h = (TH1F*)h_to_clone->Clone();
  leg->SetTextSize(0.55*h->GetXaxis()->GetTitleSize());
  h->SetLineStyle(1);
  h->SetLineColor(14);
  h->SetMarkerStyle(kFullCircle);
  
  leg->AddEntry(h, entries[0].second, "lep");
  for(int i = 1; i < (int) entries.size()-2; i++){
    h->SetLineColor(entries[i].first);
    if(i == 7) h->SetLineStyle(7); // Uncorr MC
    else h->SetLineStyle(1);
    if(i == 8){
      h->SetMarkerColor(entries[i].first); // Ratio to Uncorr MC
      h->SetMarkerStyle(kFullCircle);
      leg->AddEntry(h->Clone(), entries[i].second, "lep");
      continue;
    }
    if(i == 9) {
      h->SetMarkerColor(entries[i].first);
      h->SetMarkerStyle(kFullCircle);
      leg->AddEntry(h->Clone(), entries[i].second, "lep");
      continue;
    }
    leg->AddEntry(h->Clone(), entries[i].second, "l");                   
  }
  leg->SetTextSize(1.3*leg->GetTextSize());
  leg->Draw();

}

void RatioLegend(){
  TLegend* leg = new TLegend(.12, .13, .47, .18);
  //  leg->SetFillStyle(0);
  leg->SetNColumns(2);
  leg->SetBorderSize(1);
  TH1* h = new TH1F("","", 1, 0, 1);
  h->SetMarkerColor(17);
  h->SetLineColor(17);
  h->SetMarkerStyle(kFullCircle);
  leg->AddEntry(h->Clone(), "Uncorr. MC", "lep");
  h->SetMarkerColor(1);
  h->SetLineColor(1);
  leg->AddEntry(h->Clone(), "Decomp MC", "lep");
  leg->Draw();
}

void NeutrinoMode(std::string mode){
  std::string txt = mode == "fhc"? "Neutrino Mode" : "Antineutrino Mode";
  TLatex* lb = new TLatex(.1, .95, txt.c_str());
  lb->SetNDC();
  lb->SetTextAlign(13);
  lb->SetTextColor(kGray+3);
  lb->SetTextSize(2/45.);
  lb->Draw();

}