plot_datamc_ND_numu.C

Go to the documentation of this file.

#include "CAFAna/Core/LoadFromFile.h"

#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"

#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Analysis/Style.h"
using namespace ana;

#include "OscLib/IOscCalc.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TLatex.h"
#include "TLine.h"
#include "TLegend.h"
#include "TPad.h"
#include "TStyle.h"

#include <iostream>
#include <fstream>

//make a pretty legend 
void Legend(double x0, double y0, double x1, double y1)
{
  // Doesn't handle log-y well
  //  TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
  TLegend* leg = new TLegend(x0, y0, x1, y1);
  leg->SetFillStyle(0);

  TH1* dummy = new TH1F("", "", 1, 0, 1);
  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
  dummyFill->SetLineColor(kTotalMCColor);
  dummyFill->SetFillColor(kTotalMCErrorBandColor);

  dummy->SetMarkerStyle(kFullCircle);

  leg->AddEntry(dummy->Clone(), "ND data", "lep");
  dummy->SetLineColor(kGray+2);
  dummy->SetMarkerColor(kGray+2);
  //  leg->AddEntry(dummy->Clone(), "ND data stagger", "lep");
  dummy->SetLineColor(kTotalMCColor);
  leg->AddEntry(dummy->Clone(), "Total MC", "l");
  //leg->AddEntry(dummyFill->Clone(), "Flux Uncert.", "fl");
  dummy->SetLineColor(kNCBackgroundColor);
  leg->AddEntry(dummy->Clone(), "NC", "l");
  dummy->SetLineColor(kBeamNueBackgroundColor);
  leg->AddEntry(dummy->Clone(), "Beam #nu_{e} CC", "l");
  dummy->SetLineColor(kNumuBackgroundColor);
  leg->AddEntry(dummy->Clone(), "#nu_{#mu} CC", "l");

  leg->Draw();
}

void plot_datamc_ND_numu(std::string fname = "datamc_ND_numu_kinematics.root")
{
  struct HistDef
  {
    std::string name;
  };
  
  //which variables do you actually want to look at?
  const int kNumVars =4;
  const HistDef defs[kNumVars] = {
    "recoE",
    "recoQ2",
    "PtP",
    "CosNumi"
  };

  //which cut tiers do you want to look at?
  const int kNumSels = 1;
  const std::string selNames[kNumSels] = {
      "numuND" 
      };

  Spectrum* spects[kNumSels][kNumVars];
  IPrediction* preds[kNumSels][kNumVars];

  osc::NoOscillations noosc;
  
  //pick up spectra/prediction objects from the input root file
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const char* selName = selNames[selIdx].c_str();
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const char* varName = defs[varIdx].name.c_str();

      std::cout<<TString::Format("%s/spect_%s", selName, varName).Data()<<std::endl;
      spects[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/spect_%s", selName, varName).Data()).release();
      
      preds[selIdx][varIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/pred_%s", selName, varName).Data()).release();
      
      if(selIdx == 0 && varIdx == 0)
        std::cout << "data, MC POT: "
                  << spects[selIdx][varIdx]->POT() << " "
                  << preds[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;

    } // end for varIdx
  } // end for selIdx
    
  std::ofstream tfile;

  tfile.open("plots/TableOut-numuND.tex", ios::out);
    
  tfile << "\\begin{tabular}{llllll} \n";
  tfile << "\\multicolumn{6}{c}{numuND} \\\\ \n";
  tfile << "Sample & \\%Data/MC   & Data  & MC & NC & $\\nu_\\mu$ CC \\\\ \n";
  tfile << "\\hline \n";

  //loop over the different selector/cut levels to plot each one
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const std::string selName = selNames[selIdx];
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      
      const std::string varName = defs[varIdx].name;

      //play with the binning choice for given variables:  
      int rebinFactor=1;

      if(varName=="recoE"){
        rebinFactor=2;
      } else if(varName=="recoQ2"){
        rebinFactor=10;
      } else {
        rebinFactor=20;
      }

      //create a canvas   
      TCanvas *c2 =new TCanvas("c2","c2",500,400);      
      
      
      TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
      padAbs -> SetBottomMargin(0); 
      padAbs->SetFillStyle(0);
      padAbs -> SetLeftMargin(0.15);
      
      padAbs -> Draw();
      padAbs -> cd();

      double legendxlow=0.6;
      double legendxhigh=0.95;
      const double kPOTFD = 8.09e20;

      //Pull down the variable we want
      Spectrum* spect = spects[selIdx][varIdx];

      IPrediction* pred = preds[selIdx][varIdx];

      //Get components of the prediction and store them in histograms
      TH1* hMC = pred->Predict(&noosc).ToTH1(kPOTFD);
      hMC->SetLineColor(kTotalMCColor);
      hMC->Rebin(rebinFactor);
      
      TH1* hNC = pred->PredictComponent(&noosc,
                                        Flavors::kAll,
                                        Current::kNC,
                                        Sign::kBoth).ToTH1(kPOTFD);
      hNC->SetLineColor(kNCBackgroundColor);
      hNC->Rebin(rebinFactor);
      
      TH1* hCC = pred->PredictComponent(&noosc,
                                        Flavors::kAllNuMu,
                                        Current::kCC,
                                        Sign::kBoth).ToTH1(kPOTFD);
      hCC->SetLineColor(kNumuBackgroundColor);
      hCC->Rebin(rebinFactor);
      
      TH1* hBeam = pred->PredictComponent(&noosc,
                                          Flavors::kNuEToNuE,
                                          Current::kCC,
                                          Sign::kBoth).ToTH1(kPOTFD);
      hBeam->SetLineColor(kBeamNueBackgroundColor);
      hBeam->Rebin(rebinFactor);
      
      TH1* hData = spect->ToTH1(kPOTFD);
      hData->SetMarkerStyle(kFullCircle);
      hData->Rebin(rebinFactor);
      
      double intData = round(100.0*hData->Integral())/100.0;
      double intMC = round(100.0*hMC->Integral())/100.0;
      double intNC = round(100.0*hNC->Integral())/100.0;
      double intCC = round(100.0*hCC->Integral())/100.0;
      double intBeam = round(100.0*hBeam->Integral())/100.0;
      
      tfile.precision(4);
      tfile << std::fixed;
      tfile << "NOMINAL" << " & " << (intData/intMC) << " & ";
      tfile.precision(0);
      tfile << intData << " & ";
      tfile.precision(0);
      tfile << intMC << " & ";
      tfile.precision(0);
      tfile << intNC << " & ";
      tfile.precision(0);
      tfile << intCC << " \\\\ \n";
      tfile << "\\hline \n";
      tfile << "\\end{tabular}";
      tfile.close();
      
      hMC->GetXaxis()->CenterTitle();
      hMC->GetYaxis()->CenterTitle();
      hMC->GetXaxis()->SetDecimals();
      hMC->GetYaxis()->SetDecimals();
      hMC->GetYaxis()->SetTitleOffset(1.15);
      hMC->GetXaxis()->SetLabelSize(.0);
      hMC->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{20} POT", kPOTFD/1E20));

      if(varName=="recoE"){
        hMC->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
      }
      
      if(varName=="recoQ2" || varName=="trueQ2"){
          hMC->GetXaxis()->SetRangeUser(0,2);
      }
      
      //Draw the error band, MC, Data, MC components, and then Data one last time

      if(varName=="cvne_zoom"||varName=="shw_len"){
        legendxlow=0.2;
        legendxhigh=0.55;
      }
            
      hMC->GetXaxis()->SetNdivisions(405,kFALSE);

      hMC ->SetMinimum(.001);
      hMC ->SetMaximum(1.1*TMath::Max(hData->GetMaximum(),hMC->GetMaximum()));
      
      hMC ->DrawCopy("hist");
      hData->Draw("ep same");   
  
      hNC->Draw("hist same");
      hCC->Draw("hist same");
      hBeam->Draw("hist same");
      hMC->Draw("hist same");
      hData->Draw("ep same");   

      //Fix up the axis
      padAbs->RedrawAxis();
      padAbs->Update();
      
      //Add the legend and preliminary
      Legend(legendxlow, .5, legendxhigh, .85);
      
      TLatex* selTitle = new TLatex(.935, .95, "numuND");
      selTitle->SetNDC();
      selTitle->SetTextSize(2/30.);
      selTitle->SetTextAlign(32);
      selTitle->Draw();
      
      c2->cd();
      TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
      padRatio -> SetTopMargin(0);
      padRatio -> SetBottomMargin(0.25);
      padRatio -> SetFillStyle(0);
      padRatio -> SetLeftMargin(0.15);
      
      padRatio -> Draw();
      padRatio -> cd();
    
      TH1* hRatioMC = (TH1*)hData -> Clone("hRatioMC");
      double b = 5.0;
      
      if(varName=="recoQ2"){b=2.0;}
      if(varName=="CosNumi" || varName=="PtP"){b=1.0;}
      
      TLine* lOne = new TLine(0,1.0,b,1.0);
      lOne -> SetLineStyle(2);
      
      hRatioMC -> Divide(hMC);
      
      TH1* hRatio = (TH1*)hRatioMC -> Clone("hRatio");
      
      hRatioMC->GetXaxis()->CenterTitle();
      hRatioMC->GetYaxis()->CenterTitle();
      hRatioMC->GetYaxis()->SetLabelSize(3.7/30.);
      hRatioMC->GetXaxis()->SetLabelSize(3.7/30.);
      hRatioMC->GetXaxis()->SetLabelOffset(.04);
      hRatioMC->GetXaxis()->SetTickSize(.06);
      hRatioMC->GetYaxis()->SetTitleSize(4.5/30.);
      hRatioMC->GetXaxis()->SetTitleSize(4.5/30.);
      hRatioMC->GetYaxis()->SetRangeUser(0.5,1.5);
      hRatioMC->SetMarkerSize(.1);
      hRatioMC->GetXaxis()->SetDecimals();
      hRatioMC->GetYaxis()->SetDecimals();
      hRatioMC->GetYaxis()->SetTitleOffset(0.4);
      hRatioMC->GetYaxis()->SetTitle("Data/MC");
      
      if(varName=="recoQ2" || varName=="trueQ2"){
          hRatioMC->GetXaxis()->SetRangeUser(0,2);
      }

      if(varName=="recoE"){
        hRatioMC->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
      }
      
      hRatioMC->GetYaxis()->SetNdivisions(502,kFALSE);
      hRatioMC->GetXaxis()->SetNdivisions(405,kFALSE);

      hRatioMC -> Draw("ep");
      lOne -> Draw("same");
      
      c2->Print(("plots/"+selName+"_"+varName+".png").c_str());
      
    } // end for varIdx
  } // end for selIdx
}