plot_nominal_spectra.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#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 "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "3FlavorAna/Prediction/PredictionAddRock.h"
#include "3FlavorAna/Ana2017/nue/joint_fit_2017_tools.h"
#include "OscLib/IOscCalc.h"

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

using namespace ana;

void PlotNuePredictionFourBins (const IPrediction* pred, osc::IOscCalc* calc, double fdpot);
void TableNuePredictionFourBins(const IPrediction* pred, osc::IOscCalc* calc, double fdpot);
void PlotNueDecompFourBins     (MichelDecomp* combo = 0, ProportionalDecomp * prop = 0, 
                                bool mccomponly = false, bool overlay = false);
void TableNueNDDataMC    (MichelDecomp* combo);
void TableNueNDComponents(MichelDecomp* combo, ProportionalDecomp * prop = 0);

void plot_nominal_spectra(){


  for(std::string decomp :{"noextrap","prop","combo"}){
    auto pred = GetNuePrediction2017(decomp, DefaultOscCalc());
    
    // FD plots
    auto calc = DefaultOscCalc();
    SetFakeCalc(calc, 0.402, 2.7e-3, 1.5);
    
    double fdpot = kAna2017POT;
    
    PlotNuePredictionFourBins(pred, calc, fdpot);
    gPad->Print((TString)"plot_nue_fd_"+decomp+"_rock_spectrum.pdf");
    std::cout << decomp ;    
    TableNuePredictionFourBins(pred, calc, fdpot);
  }
  // ND Plots
  //Add prop to the mix
  auto pred_file = new TFile("/nova/ana/nu_e_ana/Ana2017/Predictions/provisional/pred_nue_reduced_v4.root", "read");
  auto xp_dir = pred_file->GetDirectory((TString) "pred_" + "xp_combo" + 
                                          "_Nominal/noShift/");
  auto xp_dir2 = pred_file->GetDirectory((TString) "pred_" + "xp_prop" + 
                                          "_Nominal/noShift/");
  TString deco_dir  ="extrap/EEextrap/Decomp";
  auto combo = MichelDecomp::LoadFrom (xp_dir, deco_dir).release();
  auto prop  = ProportionalDecomp::LoadFrom(xp_dir2, deco_dir).release();


  PlotNueDecompFourBins (combo,0,true,false);
  gPad->Print("plot_nue_nd_nominal_mconly.pdf");

  PlotNueDecompFourBins (combo,0,false,true);
  gPad->Print("plot_nue_nd_nominal_combo_vs_mc.pdf");

  PlotNueDecompFourBins (combo,prop,false,true);
  gPad->Print("plot_nue_nd_nominal_prop_vs_mc.pdf");


  TableNueNDDataMC(combo);

  TableNueNDComponents(combo,0);

  TableNueNDComponents(combo,prop);

}
void PrintTableFourBins(std::vector <TH1D*> mcnom, std::vector <TString> labels){
  
  assert(mcnom.size()==labels.size());  
  
  int N = 4;
  int Nbins = mcnom[0]->GetNbinsX();
  int binNbins = Nbins / N;
  for(TString b:{"Bin1","Bin2","Bin3","Bin4","AllBins"}){
    std::cout << "\t & \t " << b;  
  }
  std::cout << "\t \\\\ \\\\hline \n";
  for(int i = 0; i < (int) mcnom.size(); ++ i){
    std::cout << labels[i] << "\t & \t";
    for(int n = 0; n < 5; ++n){    
      int bin1 = 1 + (n * binNbins);
      int bin2 = (n+1) * binNbins;
      if(n==4) {bin1 = 1; bin2 = Nbins;}
      std::cout << mcnom[i]->Integral(bin1,bin2)
                << (n == 4 ? "\t \\\\\n" : "\t & \t ");
    }//end bins
  }//end rows
  std::cout << "\n\n";
}



std::vector <TH1D*> GetNueNDComponents(IDecomp* decomp, double pot)
{
  std::vector <TH1D*> ret;
  ret.push_back( (decomp->NueComponent()  + decomp->AntiNueComponent() +
                  decomp->NumuComponent() + decomp->AntiNumuComponent() +
                  decomp->NCTotalComponent() ). ToTH1(pot,kTotalMCColor) );  
  ret.push_back( (decomp->NueComponent() + decomp->AntiNueComponent()).
                 ToTH1(pot,kBeamNueBackgroundColor) );
  ret.push_back( (decomp->NCTotalComponent()).
                 ToTH1(pot,kNCBackgroundColor) );
  ret.push_back( (decomp->NumuComponent() + decomp->AntiNumuComponent()).
                 ToTH1(pot,kNumuBackgroundColor) );
  return ret;
}

std::vector <TH1D*> GetNueNDComponentsMC(MichelDecomp* decomp, double pot)
{
  std::vector <TH1D*> ret;
  ret.push_back( (decomp->MC_NueComponent()  + decomp->MC_AntiNueComponent() +
                  decomp->MC_NumuComponent() + decomp->MC_AntiNumuComponent() +
                  decomp->MC_NCTotalComponent() ). ToTH1(pot,kTotalMCColor) );  
  ret.push_back( (decomp->MC_NueComponent() + decomp->MC_AntiNueComponent()).
                 ToTH1(pot,kBeamNueBackgroundColor) );
  ret.push_back( (decomp->MC_NCTotalComponent()).
                 ToTH1(pot,kNCBackgroundColor) );
  ret.push_back( (decomp->MC_NumuComponent() + decomp->MC_AntiNumuComponent()).
                 ToTH1(pot,kNumuBackgroundColor) );
  return ret;
}


void TableNueNDDataMC(MichelDecomp* combo)
{
  auto data = combo->Data_Component();
  auto ndpot = data.POT();
  auto hdata = data.ToTH1(ndpot);
  auto comp = GetNueNDComponentsMC (combo, ndpot);

  std::cout << "Data vs uncorrected MC\n";
  comp.insert(comp.begin(),hdata);
  std::vector <TString> labels = {"Data","Tot MC","Beam","NC","Numu"};
  std::cout << std::setprecision(0) << std::fixed; 
  PrintTableFourBins(comp,labels);
}

void TableNueNDComponents(MichelDecomp* combo, ProportionalDecomp * prop)
{
  auto data = combo->Data_Component();
  auto ndpot = data.POT();
  std::vector <TH1D*> comp;
  if(prop)  comp = GetNueNDComponents (prop, ndpot);
  else      comp = GetNueNDComponents (combo, ndpot);
  std::vector <TString> labels = {"Tot MC","Beam","NC","Numu"};
  std::cout << std::setprecision(0) << std::fixed; 
  std::cout << (prop?"prop":"combo");
  PrintTableFourBins(comp,labels);
}

void PlotNueDecompFourBins (MichelDecomp* combo, ProportionalDecomp * prop, 
                            bool mccomponly, bool overlay)
{
  auto data = combo->Data_Component();
  auto ndpot = data.POT();
  auto hdata = data.ToTH1(ndpot);
  hdata->SetMarkerStyle(kFullCircle);
  hdata->Draw("e1");

  std::vector <TH1D*> comp_decomp;
  if(prop) comp_decomp =  GetNueNDComponents(prop, ndpot);
  else comp_decomp = GetNueNDComponents(combo, ndpot); 

  if(!mccomponly) for (auto h:comp_decomp) h->Draw("hist same");

  std::vector <TH1D*> comp_mc;
  if( mccomponly || overlay) comp_mc =  GetNueNDComponentsMC(combo, ndpot); 
  if( overlay ) for (auto h:comp_mc) h->SetLineStyle(7);
  if( mccomponly || overlay) for (auto h:comp_mc) h->Draw("hist same");
  hdata->Draw("e1 same");

  hdata->GetYaxis()->SetTitle("Events / 0.5 GeV bin");
  hdata->SetMaximum(1.2*hdata->GetMaximum());
  Nue2017FourBinLabels(0.98, 0.8 * hdata->GetXaxis()->GetLabelSize());
  Nue2017FourBinDivisions();
  Nue2017FourBinAxis(hdata);
    
  auto leg = new TLegend (0.70,0.2,0.89,0.75);
  leg->SetTextSize(hdata->GetYaxis()->GetLabelSize());
  leg->SetHeader(TString::Format("#splitline{NOvA ND}{%.2f#times10^{20} POT}",
                                 ndpot/1.E20));
  leg->AddEntry(hdata, "ND Data","epl");
  leg->AddEntry(comp_decomp[0], "Total","l");
  leg->AddEntry(comp_decomp[1], "Beam #nu_{e} CC","l");
  leg->AddEntry(comp_decomp[2], "NC","l");
  leg->AddEntry(comp_decomp[3], "#nu_{#mu} CC","l");
  if(overlay){
    auto dummy = new TGraph();
    dummy->SetLineWidth(2);
    dummy->SetLineColor(kGray+2);
    dummy->SetLineStyle(7);
    leg->AddEntry(dummy, "Uncorr. MC", "l");
  }
  leg->Draw();
}

void PlotNuePredictionFourBins(const IPrediction* pred, osc::IOscCalc * calc, double fdpot){
    
  auto htot = pred->Predict(calc).ToTH1(fdpot, kTotalMCColor);
  auto hsig = pred->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).
    ToTH1(fdpot, kNueSignalColor);
  auto hbeam = pred->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).
    ToTH1(fdpot, kBeamNueBackgroundColor);
  auto hnc   = pred->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).
    ToTH1(fdpot, kNCBackgroundColor);
  auto hnumu = pred->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).
    ToTH1(fdpot, kNumuBackgroundColor);

  htot->Draw("hist");
  for(auto h:{hsig,hbeam,hnc,hnumu})h->Draw("hist same");
    
  htot->GetYaxis()->SetTitle("Events / 0.5 GeV bin");
  htot->SetMaximum(1.2*htot->GetMaximum());
  Nue2017FourBinLabels(0.98, 0.8 * htot->GetXaxis()->GetLabelSize());
  Nue2017FourBinDivisions();
  Nue2017FourBinAxis(htot);
    
  auto leg = new TLegend (0.12,0.4,0.35,0.75);
  leg->SetTextSize(htot->GetYaxis()->GetLabelSize());
  leg->SetHeader(TString::Format("NOvA FD, %.2f #times10^{20} POT eq.",
                                 kAna2017FullDetEquivPOT/1.E20));
  leg->AddEntry(htot, "Nominal Osc. Pred.","l");
  leg->AddEntry(hsig, "Signal","l");
  leg->AddEntry(hbeam, "Beam #nu_{e} CC","l");
  leg->AddEntry(hnc, "NC","l");
  leg->AddEntry(hnumu, "#nu_{#mu} CC","l");
  leg->Draw();
}

void TableNuePredictionFourBins(const IPrediction* pred, osc::IOscCalc *calc, double fdpot)
{
  auto htot = pred->Predict(calc).ToTH1(fdpot, kTotalMCColor);
  auto hsig = pred->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).
    ToTH1(fdpot, kNueSignalColor);
  auto hbkg = (TH1D*) htot->Clone(UniqueName().c_str());
  hbkg->Add(hsig,-1);
  auto hbeam = pred->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).
    ToTH1(fdpot, kBeamNueBackgroundColor);
  auto hnc   = pred->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).
    ToTH1(fdpot, kNCBackgroundColor);
  auto hnumu = pred->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).
    ToTH1(fdpot, kNumuBackgroundColor);
  auto htau = pred->PredictComponent(calc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).
    ToTH1(fdpot, kNumuBackgroundColor);


  std::vector <TH1D*> mcnom = {htot,hsig,hbkg,hbeam,hnc,hnumu, htau};
  std::vector <TString> labels = {"Total", "Signal","Bkg.","Beam","NC","Numu","Tau"};

  std::cout << std::setprecision(3) ;
  PrintTableFourBins(mcnom,labels);
}