plot_nueFD_Signal_REWvsNOM_pTExtrap_RHC.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 "3FlavorAna/Plotting/NuePlotStyle.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/3flavorCuts.h"
#include "3FlavorAna/Cuts/NueCuts2020.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 "3FlavorAna/Prediction/PredictionExtendToPeripheral.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 "3FlavorAna/Plotting/NuePlotStyle.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 "TStyle.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);
  leg->SetTextSize(2/30.);

  //  leg->SetHeader("#splitline{#bf{Neutrino Beam}}{#splitline{FD prediction}{#nu_{e} signal+WS}}");
  leg->SetHeader("#splitline{#bf{Neutrino Beam}}{#nu_{e} signal+WS}");
  TH1* dummy = new TH1F("", "", 1, 0, 1);
  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
  leg->AddEntry(dummy->Clone(), "","");
  dummyFill->SetLineColor(kTotalMCColor);
  dummyFill->SetFillColor(kTotalMCErrorBandColor);

  dummy->SetLineColor(kNueSignalColor);
  leg->AddEntry(dummy->Clone(), "Nominal", "l");
  dummy->SetLineStyle(9);
  dummy->SetLineColor(kOrange+1);
  leg->AddEntry(dummy->Clone(), "Q^{2} reweighted", "l");
  dummy->SetLineStyle(8);
  dummy->SetLineColor(kGray+2);
  leg->AddEntry(dummy->Clone(), "p_{T}/p reweighted", "l");
  dummy->SetLineStyle(3);
  //  dummy->SetLineColor(kGray+2);  
  leg->AddEntry(dummy->Clone(), "cos #theta reweighted", "l");  

  leg->Draw();
}

TH1* create_coreperi_TH1(const TH1* coreTH1, const TH1* periTH1)
{
  TH1* coreperiTH1 = (TH1*)coreTH1->Clone("coreperi");
  for(int i = 18; i < 28; ++i){
    coreperiTH1->SetBinContent(i, periTH1->GetBinContent(i));
  }
  
  return coreperiTH1;
}

void plot_nueFD_Signal_REWvsNOM_pTExtrap_RHC(const std::string filename_no, const std::string filename_re_core, const std::string filename_re_peri)

{
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(0);

  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_core = new TFile (filename_re_core.c_str(),"READ");
  if(file_re_core->IsZombie()){std::cerr << "Can't find " << filename_re_core.c_str() << std::endl; return;}
  TFile * file_re_peri = new TFile (filename_re_peri.c_str(), "READ");
  if(file_re_peri->IsZombie()){std::cerr << "Can't find " << filename_re_peri.c_str() << std::endl; return;}   

  osc::OscCalcDumb mycalc;

  TFile* file = new TFile("FD_KinematicsCorrection_pTExtrap_RHC.root","RECREATE");

  TDirectory*  dpred_noextrap = file_no->GetDirectory("NOReweight/nueAxis_NoExtrap");
  PredictionNoExtrap* pred_noextrap  = ana::LoadFrom<PredictionNoExtrap>(dpred_noextrap).release();       
 
  //Save the predictions for each quant for nom, core and peri    

  std::vector<const PredictionExtendToPeripheral* > pred_no_quant;
  std::vector<const PredictionExtendToPeripheral* > pred_re_core_q2_quant;
  std::vector<const PredictionExtendToPeripheral* > pred_re_core_PtP_quant;
  std::vector<const PredictionExtendToPeripheral* > pred_re_core_Cos_quant;
  std::vector<const PredictionExtendToPeripheral* > pred_re_peri_q2_quant;
  std::vector<const PredictionExtendToPeripheral* > pred_re_peri_PtP_quant;
  std::vector<const PredictionExtendToPeripheral* > pred_re_peri_Cos_quant;

  const unsigned int nQuants = 3;

  for ( unsigned int quant = 0; quant < nQuants; ++quant ){
    //nominal
    TDirectory* dpred_no_quant = file_no->GetDirectory(Form("NOReweight/nueAxis_NueSignalExtrap_Quant%i",quant+1));
    auto extend_nom = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_no_quant).release()), pred_noextrap, kNue2020Binning, false);
    pred_no_quant.push_back(extend_nom);

    //core
    TDirectory* dpred_re_core_q2_quant = file_re_core->GetDirectory(Form("CORE/nueAxis_NueSignalExtrap_Q2Weight_Quant%i", quant+1));
    auto extend_core_q2 = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_core_q2_quant).release()), pred_noextrap, kNue2020Binning, false);
    pred_re_core_q2_quant.push_back(extend_core_q2);

    TDirectory* dpred_re_core_PtP_quant = file_re_core->GetDirectory(Form("CORE/nueAxis_NueSignalExtrap_PtPWeight_Quant%i", quant+1));
    auto extend_core_PtP = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_core_PtP_quant).release()), pred_noextrap, kNue2020Binning, false);
    pred_re_core_PtP_quant.push_back(extend_core_PtP);

    TDirectory* dpred_re_core_Cos_quant = file_re_core->GetDirectory(Form("CORE/nueAxis_NueSignalExtrap_CosWeight_Quant%i", quant+1));
    auto extend_core_Cos = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_core_Cos_quant).release()), pred_noextrap, kNue2020Binning, false);
    pred_re_core_Cos_quant.push_back(extend_core_Cos);

    //peri
    TDirectory* dpred_re_peri_q2_quant = file_re_peri->GetDirectory(Form("PERI/nueAxis_NueSignalExtrap_Q2Weight_Quant%i", quant+1));
    auto extend_peri_q2 = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_peri_q2_quant).release()), pred_noextrap, kNue2020Binning, false);
    pred_re_peri_q2_quant.push_back(extend_peri_q2);

    TDirectory* dpred_re_peri_PtP_quant = file_re_peri->GetDirectory(Form("PERI/nueAxis_NueSignalExtrap_PtPWeight_Quant%i", quant+1));
    auto extend_peri_PtP = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_peri_PtP_quant).release()), pred_noextrap, kNue2020Binning, false);
    pred_re_peri_PtP_quant.push_back(extend_peri_PtP);

    TDirectory* dpred_re_peri_Cos_quant = file_re_peri->GetDirectory(Form("PERI/nueAxis_NueSignalExtrap_CosWeight_Quant%i", quant+1));
    auto extend_peri_Cos = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_peri_Cos_quant).release()), pred_noextrap, kNue2020Binning, false);
    pred_re_peri_Cos_quant.push_back(extend_peri_Cos);    
  }

  //Add quants to form a single prediction for nom, core and peri
  //add to form Spectrums (so mycalc already implemented)

  Spectrum spec_nom      = pred_no_quant.front()->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);
  spec_nom.Clear();

  Spectrum spec_core_q2  = pred_re_core_q2_quant.front()->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);
  spec_core_q2.Clear();

  Spectrum spec_core_PtP = pred_re_core_PtP_quant.front()->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);
  spec_core_PtP.Clear();

  Spectrum spec_core_Cos = pred_re_core_Cos_quant.front()->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);
  spec_core_Cos.Clear();

  Spectrum spec_peri_q2  = pred_re_peri_q2_quant.front()->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);
  spec_peri_q2.Clear();

  Spectrum spec_peri_PtP = pred_re_peri_PtP_quant.front()->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);
  spec_peri_PtP.Clear();

  Spectrum spec_peri_Cos = pred_re_peri_Cos_quant.front()->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);
  spec_peri_Cos.Clear();



  for (unsigned int quant = 0; quant < nQuants; ++quant ){
    spec_nom      += pred_no_quant[quant]->PredictComponent(&mycalc,
                                            Flavors::kNuMuToNuE,
                                            Current::kCC,
                                            Sign::kBoth);
    
    spec_core_q2  += pred_re_core_q2_quant[quant]->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);

    spec_core_PtP += pred_re_core_PtP_quant[quant]->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);

    spec_core_Cos += pred_re_core_Cos_quant[quant]->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);

    spec_peri_q2  += pred_re_peri_q2_quant[quant]->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);

    spec_peri_PtP += pred_re_peri_PtP_quant[quant]->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);

    spec_peri_Cos += pred_re_peri_Cos_quant[quant]->PredictComponent(&mycalc,
                                                 Flavors::kNuMuToNuE,
                                                 Current::kCC,
                                                 Sign::kBoth);

  }
  
  //form histograms and merge core+peri to form a single prediction for FD
  // POT set by hand from total POT of ND sample
  const double kPOTFD = kAna2019RHCPOT;

  TH1* hAll_no       = spec_nom.ToTH1(kPOTFD);

  TH1* hAll_core_q2  = spec_core_q2.ToTH1(kPOTFD);
  TH1* hAll_peri_q2  = spec_peri_q2.ToTH1(kPOTFD);
  TH1* hAll_re_q2    = (TH1*)(create_coreperi_TH1(hAll_core_q2, hAll_peri_q2))->Clone("hAll_re_q2");  

  TH1* hAll_core_PtP = spec_core_PtP.ToTH1(kPOTFD);
  TH1* hAll_peri_PtP = spec_peri_PtP.ToTH1(kPOTFD);
  TH1* hAll_re_PtP   = (TH1*)(create_coreperi_TH1(hAll_core_PtP, hAll_peri_PtP))->Clone("hAll_re_PtP");  

  TH1* hAll_core_Cos = spec_core_Cos.ToTH1(kPOTFD);
  TH1* hAll_peri_Cos = spec_peri_Cos.ToTH1(kPOTFD);
  TH1* hAll_re_Cos   = (TH1*)(create_coreperi_TH1(hAll_core_Cos, hAll_peri_Cos))->Clone("hAll_re_Cos");  

  //Set hist styles
  hAll_re_q2->SetLineColor(kOrange+1);
  hAll_re_q2->SetLineStyle(9);
  //      hAll_re_q2->SetLineWidth(1);

  hAll_re_PtP->SetLineColor(kGray+2);
  hAll_re_PtP->SetLineStyle(8);
  hAll_re_PtP->SetLineWidth(1);

  hAll_re_Cos->SetLineColor(kGray+2);
  hAll_re_Cos->SetLineStyle(3); 
  hAll_re_Cos->SetLineWidth(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();
        
  hAll_re_q2->GetXaxis()->CenterTitle();
  hAll_re_q2->GetYaxis()->CenterTitle();
  hAll_re_q2->GetXaxis()->SetDecimals();
  hAll_re_q2->GetYaxis()->SetDecimals();
  hAll_re_q2->GetYaxis()->SetTitleOffset(.8);
  hAll_re_q2->GetYaxis()->SetTitleSize(.06);
  hAll_re_q2->GetYaxis()->SetLabelSize(.06);
  hAll_re_q2->GetXaxis()->SetLabelSize(.0);
  hAll_re_q2->SetMaximum(1.1*hAll_re_q2->GetMaximum());
  hAll_re_q2->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{21} POT", kPOTFD/1E21));
    
  hAll_no->Draw("hist");
  hAll_re_PtP->Draw("same hist");
  hAll_re_Cos->Draw("same hist");
  hAll_re_q2->Draw("same hist");

  Nue2018ThreeBinAxis(hAll_re_q2, true, false, false);
  hAll_re_q2->GetYaxis()->SetRangeUser(.001,20);
  
  //Fix up the axis
  padAbs->RedrawAxis();
  padAbs->Update();
        
  double legendxlow=0.55;
  double legendxhigh=0.95;
        
  Legend(legendxlow, .5, legendxhigh, .85);
      
  TLatex* selTitle = new TLatex(.9, .95, "NOvA Simulation");
  selTitle->SetNDC();
  selTitle->SetTextColor(kGray+1);
  selTitle->SetTextSize(2/30.);
  selTitle->SetTextAlign(32);
  selTitle->Draw();
  
  c2->cd();
  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.0, 1, 0.33);
  padRatio -> SetTopMargin(0.11);
  padRatio -> SetBottomMargin(0.33);
  padRatio -> SetFillStyle(0);
  padRatio -> SetLeftMargin(0.15);
        
  padRatio -> Draw();
  padRatio -> cd();
    
  TH1* hRatio_q2  = (TH1*)hAll_re_q2 -> Clone("hRatio_q2");
  TH1* hRatio_PtP = (TH1*)hAll_re_PtP -> Clone("hRatio_PtP");
  TH1* hRatio_Cos = (TH1*)hAll_re_Cos -> Clone("hRatio_Cos");
  TLine* lOne = new TLine(0,1.0,27.0,1.0);
  lOne -> SetLineStyle(2);
  
  hRatio_q2 -> Divide(hAll_no);
  hRatio_PtP -> Divide(hAll_no);
  hRatio_Cos -> Divide(hAll_no);
      
  hRatio_q2->GetXaxis()->CenterTitle();
  hRatio_q2->GetYaxis()->CenterTitle();
  hRatio_q2->GetYaxis()->SetLabelSize(.12);
  hRatio_q2->GetXaxis()->SetLabelSize(.12);      
  hRatio_q2->GetXaxis()->SetTitleOffset(1.0);
  hRatio_q2->GetXaxis()->SetTickSize(.1);
  hRatio_q2->GetYaxis()->SetTitleSize(.12);
  hRatio_q2->GetXaxis()->SetTitleSize(.12);
  hRatio_q2->SetMarkerSize(.1);
  hRatio_q2->GetXaxis()->SetDecimals();
  hRatio_q2->GetYaxis()->SetDecimals();
  hRatio_q2->GetYaxis()->SetTitleOffset(0.4);
  hRatio_q2->GetYaxis()->SetTitle("#splitline{Ratio to}{nominal}");
      
  hRatio_q2->GetYaxis()->SetNdivisions(502,kFALSE);
  hRatio_q2 -> Draw("hist");
  hRatio_PtP -> Draw("same hist");
  hRatio_Cos -> Draw("same hist");
  lOne -> Draw("same");
  hRatio_q2 -> SetMinimum(0.8);
  Nue2018ThreeBinAxis(hRatio_q2, true, false, false);
  hRatio_q2->GetYaxis()->SetRangeUser(0.95,1.05);
        
  c2->Print(("plots/FDnue_pred_recoE_pTExtrap_RHC.png"));
  c2->SaveAs("plots/FDnue_pred_recoE_pTExtrap_RHC.pdf","pdf");
      
  hRatio_q2-> Write("trueQ2_FD_KinematicsCorrection_RHC");
  hRatio_PtP-> Write("PtP_FD_KinematicsCorrection_RHC");
  hRatio_Cos-> Write("Cos_FD_KinematicsCorrection_RHC");
  //  file->Close();

  //Save the integral stats in a tex file
  double intAll_no = hAll_no->Integral();
  double intAll_re_q2 = hAll_re_q2->Integral();
  double intAll_re_PtP = hAll_re_PtP->Integral();
  double intAll_re_Cos = hAll_re_Cos->Integral();
  
  std::cout << intAll_no << std::endl;
  std::cout << intAll_re_q2 << std::endl;
  std::cout << intAll_re_PtP << std::endl;
  std::cout << intAll_re_Cos << std::endl;
    
  std::ofstream tfile;
  
  tfile.open("plots/FDextrap-TableOut-pTExtrap-RHC.tex", ios::out);
    
  tfile << "\\begin{tabular}{lll} \n";
  tfile << "\\multicolumn{3}{c}{prop FD $\nu_e$ signal + WS} \\\\ \n";
  tfile << "Sample & Signal+WS & /NOMINAL  \\\\ \n";
  tfile << "\\hline \n";  
        
  tfile << std::fixed;
  tfile.precision(2);
  tfile << "NOMINAL"     << " & " << round(100.0*intAll_no)/100.0 << " & \\\\ \n";
  tfile << "trueQ2 REW." << " & " << round(100.0*intAll_re_q2)/100.0  << " & " << round(1000.0*((intAll_re_q2/intAll_no)-1.0))/10.0  << " \\% \\\\ \n";
  tfile << "PtP REW."    << " & " << round(100.0*intAll_re_PtP)/100.0 << " & " << round(1000.0*((intAll_re_PtP/intAll_no)-1.0))/10.0 << " \\% \\\\ \n";
  tfile << "Cos REW."    << " & " << round(100.0*intAll_re_Cos)/100.0 << " & " << round(1000.0*((intAll_re_Cos/intAll_no)-1.0))/10.0 << " \\% \\\\ \n";
  tfile << "\\hline \n";
  tfile << "\\end{tabular}";
  tfile.close();

} // end