plot_nueFD_Signal_REWvsNOM.C

Go to the documentation of this file.

// Macro to make table of nominal predictions at oscillation values

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"

#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Analysis/Style.h"

//#include "CAFAna/nue/SecondAna/draw_plots_util.h"

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

using namespace ana;

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->SetLineColor(kNueSignalColor);
  leg->AddEntry(dummy->Clone(), "#nu_{e} signal NOM.", "l");
  dummy->SetLineStyle(3);
  leg->AddEntry(dummy->Clone(), "trueQ2 REW.", "l");
  dummy->SetLineStyle(8);
  leg->AddEntry(dummy->Clone(), "PtP REW.", "l");
  dummy->SetLineStyle(9);
  leg->AddEntry(dummy->Clone(), "CosTheta REW.", "l");  

  leg->Draw();
}

void plot_nueFD_Signal_REWvsNOM(const std::string filename_no, const std::string filename_re
                                      )

{
  const int kNumSels = 1;
  const int kNumVars = 1;
    
  struct HistDef
  {
    std::string name;
  };
  
  const HistDef defs[kNumVars] = {
    "recoE",
  };
  
  const std::string selNames[kNumSels] = {
      "AllSamples" 
      };  
  
  TFile * file_no = new TFile (filename_no.c_str(),"READ");
  if(file_no->IsZombie()){std::cerr << "Can't find " << filename_no.c_str() << std::endl; return;}
  TFile * file_re = new TFile (filename_re.c_str(),"READ");
  if(file_re->IsZombie()){std::cerr << "Can't find " << filename_re.c_str() << std::endl; return;}
    
  osc::OscCalcPMNSOpt mycalc;
  ResetOscCalcToDefault(&mycalc);
  
  mycalc.SetDmsq32(+fabs(mycalc.GetDmsq32()));
  mycalc.SetdCP(3*TMath::Pi()/2);
  mycalc.SetTh23(asin(sqrt(0.404)));
  
  std::ofstream tfile;
  const int kNumTabs = 1;
  const std::string tabNames[kNumTabs] = {
      "AllSamples"
  };
  
    
    tfile.open("plots/FDextrap-TableOut.tex", ios::out);
    
    tfile << "\\begin{tabular}{lllll} \n";
    tfile << "\\multicolumn{5}{c}{prop FD $\nu_e$ signal ALL SAMPLES} \\\\ \n";
    tfile << "Sample & Sig.+Bkg.  & /NOMINAL  & Signal & /NOMINAL  \\\\ \n";
    tfile << "\\hline \n";  
  
  double intTotalFull, intAllFull, intSigFull, intNCFull, intCCFull, intBeamFull, intTauFull;
  
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const std::string selName = selNames[selIdx];
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      bool pidaxes = varIdx;
      const std::string varName = defs[varIdx].name;
      TDirectory * dpred_no = file_no->GetDirectory("prediction/nue_pred_AllSamples_recoE");
      TDirectory * dpred_re = file_re->GetDirectory("prediction/nue_pred_trueQ2_AllSamples_recoE");
      TDirectory * dpred_re_PtP = file_re->GetDirectory("prediction/nue_pred_PtP_AllSamples_recoE");
      TDirectory * dpred_re_Cos = file_re->GetDirectory("prediction/nue_pred_Cos_AllSamples_recoE");
        
      std::cout << "HELLO" << endl;
      //FD Plots
      auto pred_no = ana::LoadFrom<IPrediction> (dpred_no);      
      auto pred_re = ana::LoadFrom<IPrediction> (dpred_re);  
      auto pred_re_PtP = ana::LoadFrom<IPrediction> (dpred_re_PtP);
      auto pred_re_Cos = ana::LoadFrom<IPrediction> (dpred_re_Cos); 
   
      std::cout << "HEllO" << endl;
      if(selIdx == 0 && varIdx == 0)
          std::cout << "MC POT: "
                  << pred_no->Predict(&mycalc).POT() << std::endl;
                  
      //play with the binning choice for given variables:  
      int rebinFactor=1;
    
      //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();
    
      // POT set by hand from total POT of ND sample
      const double kPOTFD = 8.09e20;
        
      //Backgrounds                                         
      TH1* hNC = pred_no.get()->PredictComponent(&mycalc,
                                             Flavors::kAll,
                                             Current::kNC,
                                             Sign::kBoth).ToTH1(kPOTFD);
                                             
      TH1* hCC = pred_no.get()->PredictComponent(&mycalc,
                                             Flavors::kAllNuMu,
                                             Current::kCC,
                                             Sign::kBoth).ToTH1(kPOTFD);
                                             
      TH1* hBeam = pred_no.get()->PredictComponent(&mycalc,
                                               Flavors::kNuEToNuE,
                                               Current::kCC,
                                               Sign::kBoth).ToTH1(kPOTFD);
        
      TH1* hAll = (TH1*)hNC->Clone(UniqueName().c_str());
      hAll->Add(hCC);
      hAll->Add(hBeam);
        
      //Signals
      TH1* hSig_no = pred_no.get()->PredictComponent(&mycalc,
                                              Flavors::kNuMuToNuE,
                                              Current::kCC,
                                              Sign::kBoth).ToTH1(kPOTFD);
      hSig_no->SetLineColor(kNueSignalColor);
      hSig_no->Rebin(rebinFactor);
        
      TH1* hSig_re = pred_re.get()->PredictComponent(&mycalc,
                                              Flavors::kNuMuToNuE,
                                              Current::kCC,
                                              Sign::kBoth).ToTH1(kPOTFD);
      hSig_re->SetLineColor(kNueSignalColor);
      hSig_re->SetLineStyle(3);
      hSig_re->Rebin(rebinFactor);
        
      TH1* hSig_re_PtP = pred_re_PtP.get()->PredictComponent(&mycalc,
                                              Flavors::kNuMuToNuE,
                                              Current::kCC,
                                              Sign::kBoth).ToTH1(kPOTFD);
      hSig_re_PtP->SetLineColor(kNueSignalColor);
      hSig_re_PtP->SetLineStyle(8);
      hSig_re_PtP->Rebin(rebinFactor);
                                        
      TH1* hSig_re_Cos = pred_re_Cos.get()->PredictComponent(&mycalc,
                                              Flavors::kNuMuToNuE,
                                              Current::kCC,
                                              Sign::kBoth).ToTH1(kPOTFD);
      hSig_re_Cos->SetLineColor(kNueSignalColor);
      hSig_re_Cos->SetLineStyle(9);
      hSig_re_Cos->Rebin(rebinFactor);
      
      double intAll = hAll->Integral();
      double intSig_no = hSig_no->Integral();
      double intSig_re = hSig_re->Integral();
      double intSig_re_PtP = hSig_re_PtP->Integral();
      double intSig_re_Cos = hSig_re_Cos->Integral();
        
      tfile << std::fixed;
      tfile.precision(2);
      tfile << "NOMINAL" << " & " << round(100.0*(intSig_no+intAll))/100.0 << " & 0.00 & ";
      tfile << round(100.0*intSig_no)/100.0 << " & 0.00 \\\\ \n";
      tfile << "trueQ2 REW." << " & " << round(100.0*(intSig_re+intAll))/100.0 << " & " << round(1000.0*((intSig_re+intAll)/(intSig_no+intAll)-1.0))/10.0 << " \\% & ";
      tfile << round(100.0*intSig_re)/100.0 << " & "<< round(1000.0*((intSig_re/intSig_no)-1.0))/10.0 << " \\% \\\\ \n";
      tfile << "PtP REW." << " & " << round(100.0*(intSig_re_PtP+intAll))/100.0 << " & " << round(1000.0*((intSig_re_PtP+intAll)/(intSig_no+intAll)-1.0))/10.0 << " \\% & ";
      tfile << round(100.0*intSig_re_PtP)/100.0 << " & " << round(1000.0*((intSig_re_PtP/intSig_no)-1.0))/10.0 << " \\% \\\\ \n";
      tfile << "Cos REW." << " & " << round(100.0*(intSig_re_Cos+intAll))/100.0 << " & " << round(1000.0*((intSig_re_Cos+intAll)/(intSig_no+intAll)-1.0))/10.0 << " \\% & ";
      tfile << round(100.0*intSig_re_Cos)/100.0 << " & " << round(1000.0*((intSig_re_Cos/intSig_no))-1.0)/10.0 << " \\% \\\\ \n";
      tfile << "\\hline \n";
      tfile << "\\end{tabular}";
      tfile.close();
        
      hSig_re->GetXaxis()->CenterTitle();
      hSig_re->GetYaxis()->CenterTitle();
      hSig_re->GetXaxis()->SetDecimals();
      hSig_re->GetYaxis()->SetDecimals();
      hSig_re->GetYaxis()->SetTitleOffset(1.15);
      hSig_re->GetYaxis()->SetRangeUser(.001,1.005*TMath::Max(hSig_re->GetMaximum(),hSig_no->GetMaximum()));
      hSig_re->GetYaxis()->SetLabelSize(3.7/70.);
      hSig_re->GetXaxis()->SetLabelSize(.0);
      hSig_re->SetMaximum(1.1*TMath::Max(hSig_re->GetMaximum(),hSig_no->GetMaximum()));
      hSig_re->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{20} POT", kPOTFD/1E20));
        
      if(varName=="recoE"){
        hSig_re->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
      }
        
      hSig_re->GetXaxis()->SetNdivisions(405,kFALSE);
    
      hSig_re->DrawCopy("hist");
      hSig_no->Draw("same hist");
      hSig_re_PtP->Draw("same hist");
      hSig_re_Cos->Draw("same hist");
        
      //Fix up the axis
      padAbs->RedrawAxis();
      padAbs->Update();
        
      double legendxlow=0.5;
      double legendxhigh=0.95;
        
      Legend(legendxlow, .55, legendxhigh, .85);
      
      TLatex* selTitle = new TLatex(.935, .95, "FD #nu_{e} signal reweighted");
      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* hRatio = (TH1*)hSig_re -> Clone("hRatio");
      TH1* hRatio_PtP = (TH1*)hSig_re_PtP -> Clone("hRatio_PtP");
      TH1* hRatio_Cos = (TH1*)hSig_re_Cos -> Clone("hRatio_Cos");
      TLine* lOne = new TLine(0,1.0,5.0,1.0);
      lOne -> SetLineStyle(2);
      
      hRatio -> Divide(hSig_no);
      hRatio_PtP -> Divide(hSig_no);
      hRatio_Cos -> Divide(hSig_no);
      
      hRatio->GetXaxis()->CenterTitle();
      hRatio->GetYaxis()->CenterTitle();
      hRatio->GetYaxis()->SetLabelSize(3.7/30.);
      hRatio->GetXaxis()->SetLabelSize(3.7/30.);
      hRatio->GetXaxis()->SetLabelOffset(.04);
      hRatio->GetXaxis()->SetTickSize(.06);
      hRatio->GetYaxis()->SetTitleSize(4.5/30.);
      hRatio->GetXaxis()->SetTitleSize(4.5/30.);
      hRatio->GetYaxis()->SetRangeUser(0.7,1.3);
      hRatio->SetMarkerSize(.1);
      hRatio->GetXaxis()->SetDecimals();
      hRatio->GetYaxis()->SetDecimals();
      hRatio->GetYaxis()->SetTitleOffset(0.4);
      hRatio->GetYaxis()->SetTitle("REW./NOM.");
      
      if(varName=="recoE"){
        hRatio->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
      }
      
      hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
      hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
      
      hRatio -> Draw("hist");
      hRatio_PtP -> Draw("hist same");
      hRatio_Cos -> Draw("hist same");
      lOne -> Draw("same");
        
      c2->Print(("plots/FDnue_pred_recoE_abs.png"));
    } // end for varIdx
  } // end for selIdx
} // end