plot_pid.C

Go to the documentation of this file.

///////////////////////////////////////////////////////////
//
// Script to plot data/mc comparisons for prong CVN with
// different cuts.
//
//////////////////////////////////////////////////////////

#include <sstream>
#include <iomanip>
#include <fstream>

#include "TCanvas.h"
#include "TH1.h"
#include "TF1.h"
#include "TFile.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TArrayD.h"
#include "TAxis.h"
#include "TLatex.h"
#include "TStyle.h"
#include "THStack.h"

#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/GenieWeights.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "CAFAna/reco/prongID/make_pid.h"

using namespace ana;

void HCtag(TString tag,double pot)
{

  std::stringstream out;
  out << std::setprecision(3);
  out << pot;

  std::string potstr = out.str();
  std::string coef = out.str().substr(0,4);
  std::string power = out.str().substr(6,2);
  potstr = coef+"x10^{"+power+"}";

  TLatex* tex = new TLatex(.625, .85, "#scale[0.8]{"+potstr+" POT-equiv "+tag+"}");
  //if ( tag == "FHC") tex->SetTextColor(kPink+3);
  //if ( tag == "RHC") tex->SetTextColor(kAzure+3);
  tex->SetNDC();
  tex->SetTextSize(0.043);
  tex->SetTextFont(42);
  tex->SetTextAlign(12);
  tex->Draw();
}

std::string ToString(double value) {
  std::stringstream out;
  out << std::setprecision(3);
  out << value;
  return out.str();
}

void Preliminary()
{
  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
  prelim->SetTextColor(kAzure);
  prelim->SetNDC();
  prelim->SetTextSize(2/30.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}

void plot_pid(const std::string inFile,bool isFHC = 1){

  const std::string cStr = isFHC ? "FHC" : "RHC";

  Spectrum* spects[kNSels][kNPlots];

  for(unsigned int i = 0;i < kNSels;++i){
    for(unsigned int j = 0;j < kNPlots;++j){
      spects[i][j] = LoadFromFile<Spectrum>(inFile,(sels[i].suffix+"_"
                                                    +plots[j].suffix).c_str()).release();
    }
  }

  // Get spectra from file 

  TFile fout(("plot_pid_output_"+cStr+".root").c_str(),"RECREATE");

  for ( unsigned int i = 0 ; i < kNSels ; ++i ){
    for ( unsigned int j = 0 ; j < 5 ; ++j ){

      std::string theSuffix = sels[i].suffix+"_"+plots[j].suffix;
      TCanvas *c = new TCanvas(theSuffix.c_str(),theSuffix.c_str());
      TH1* h[kNPlots];
      
      h[j] = spects[i][j]->ToTH1(spects[i][j]->POT());
      
      THStack *hs = new THStack("hs","");
      for ( int k = j+5;k < kNPlots;k = k+5 ){
        h[k] = spects[i][k]->ToTH1(spects[i][j]->POT());
        hs->Add(h[k]);
      }
      
      h[j]->SetMarkerStyle(kFullCircle);
      h[j+5]->SetLineColor(PDGcolor["Electron"]);
      h[j+10]->SetLineColor(PDGcolor["Muon"]);
      h[j+15]->SetLineColor(PDGcolor["Gamma"]);
      h[j+20]->SetLineColor(PDGcolor["Pion"]);
      h[j+25]->SetLineColor(PDGcolor["Proton"]);
      
      for(int k = j+5;k < kNPlots;k = k+5) FillWithDimColor(h[k]);
      
      hs->Draw("hist");
      h[j]->Draw("ep same");
      
      hs->GetYaxis()->SetTitle(("Prongs / "+ToString(h[i]->GetBinWidth(1))+" PID").c_str());
      hs->GetXaxis()->SetTitle(plots[j].label.c_str());
      hs->GetXaxis()->CenterTitle();
      hs->GetYaxis()->CenterTitle();
      
      HCtag(cStr,spects[i][j]->POT());
      Preliminary();
      
      hs->SetMaximum(1.05*h[j]->GetMaximum());
      
      gPad->Modified();
      gPad->RedrawAxis();
      gPad->Update();
      
      TLegend *l = new TLegend(0.35,0.5,0.7,0.8);
      l->AddEntry(h[j],"Data","lpe");
      l->AddEntry(h[j+25],"True Proton","f");
      l->AddEntry(h[j+20],"True Pion","f");
      l->AddEntry(h[j+15],"True Photon","f");
      l->AddEntry(h[j+10],"True Muon","f");
      l->AddEntry(h[j+5],"True Electron","f");
      l->Draw("same");
      
      fout.cd();
      c->Write();
      c->Print(("pidPlot_"+cStr+"_"+sels[i].suffix+"_"+plots[j].suffix+".root").c_str());
      c->Print(("pidPlot_"+cStr+"_"+sels[i].suffix+"_"+plots[j].suffix+".pdf").c_str());
      gPad->SetLogy();
      c->Print(("pidPlot_"+cStr+"_"+sels[i].suffix+"_"+plots[j].suffix+"_log.pdf").c_str());
    } 
  }
}