energyResolution.C

Go to the documentation of this file.

// to get histograms of cosmics from cosmic trigger data and numi stream, and GENIE prediction, by period and by quantile
// argument is period: 1, 2, 3, or 5
// author: Kirk Bays (bays@caltech.edu) Sep 7 2017
// Histograms made: cosrej BDT inputs, neutrino E, and optimized neutrino E
// suitable to run on the grid (set bool grid to true to keep track of pot/livetime and not scale per job)
// another important option is 'quantiletype': 0 means use quantiles determined per period, 1 means use quantiles determined from all periods combined

#include "3FlavorAna/Ana2018/numu/includes.h"
#include "3FlavorAna/Ana2018/numu/varsandcuts.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Analysis/Plots.h"

#include <cmath>
#include <vector>
#include <sstream>
#include <string>
#include <iomanip>


using namespace ana;

void energyResolution(bool fhc = true, bool isFirstRun = false)
{

  gStyle->SetStatFormat("6.3g");
  gStyle->SetPaintTextFormat("6.3g");
  gStyle->SetFitFormat("6.3g");
  gStyle->SetTitleOffset(1.15,"x");
  gStyle->SetTitleOffset(0.95,"y");
  gStyle->SetPadLeftMargin(0.12);
  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetHistLineWidth(3);

  double livetime, pot;
  std::string numifiles, cosfiles, geniefiles, outfile, quantpath, horn, prodgen;
  if(fhc){
    livetime = kAna2018FHCLivetime;
    pot      = kAna2018FHCPOT;
    horn     = "fhc";
    prodgen  = ".d";
  }
  if(!fhc){
    livetime = kAna2018RHCLivetime;
    pot      = kAna2018RHCPOT;
    horn     = "rhc";
    prodgen  = ".e";
  }

  geniefiles = "prod_sumdecaf_R17-11-14-prod4reco" + prodgen + "_fd_genie_nonswap_" + horn + "_nova_v08_full_v1_numu2018";
  outfile = horn + "_energyResolution_numu2018.root";

  const Var kRMTOT([](const caf::SRProxy *sr) 
                   {
                     if(sr->mc.nnu == 0) return -10.0;
                     return (kCCE(sr)-kTrueE(sr))/kTrueE(sr);
                   });

  std::cout << "\n\nGet quantile cuts from file" << std::endl;
  quantpath = "/nova/ana/nu_mu_ana/Ana2018/Quantiles/quantiles__" + horn + "_full__numu2018.root";
  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
  TH2* quanthist = (TH2*)quantfile->FindObjectAny("FDSpec2D");
  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
  quantfile->Close();
  delete quantfile;


  if(isFirstRun){
    SpectrumLoader loaderGenie(geniefiles);
    loaderGenie.SetSpillCut(kStandardSpillCuts);

    HistAxis axisRMTOT("(E_{reco}-E_{true})/E_{true}", Binning::Simple(200,-1,1), kRMTOT);

    Spectrum* specTotal = new Spectrum(loaderGenie, axisRMTOT, kNumuCutFD2018, 
                                       kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2018);
    std::map< TString, Spectrum* > specMap;
    //  for(auto thisCut : QuantCuts){
    for(std::vector<Cut>::size_type i = 0; i != QuantCuts.size(); i++) {
      std::cout<< "i: " << i <<std::endl;
      specMap[Form("rmtot_q%d",(i+1))] = new Spectrum(loaderGenie, axisRMTOT, QuantCuts[i] && kNumuCutFD2018, 
                                                      kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2018);
    }  


    loaderGenie.Go();


    TFile fout(outfile.c_str(), "recreate");
    TH1D* histTotal = specTotal -> ToTH1(pot);
    histTotal -> Write("rmtot_comb");
    for(std::vector<Cut>::size_type i = 0; i != QuantCuts.size(); i++) {
      std::cout<< "i: " << i <<std::endl;
      TH1D* hist = specMap[Form("rmtot_q%d",(i+1))] -> ToTH1(pot);
      hist -> Write(Form("rmtot_q%d",(i+1)));    
    }
    fout.Close();

  } 


  //////////////////////////////
  // Draw hists
  DontAddDirectory guard;
  TFile fout(outfile.c_str());

  TH1D* histTotal = (TH1D*)fout.Get("rmtot_comb");

  TCanvas* c = new TCanvas();
  TLegend *leg = new TLegend(0.15,0.45,0.48,0.85);
  histTotal -> GetYaxis() -> SetRangeUser(0, 1.3 * histTotal->GetMaximum());
  histTotal -> Draw("hist");
  histTotal -> GetXaxis() -> CenterTitle();  
  double rms = histTotal -> GetRMS();
  double mean = histTotal -> GetMean();
  leg -> AddEntry(histTotal, "Combined", "l");
  leg -> AddEntry(histTotal, Form("Mean:%.03lf, RMS:%.03lf",mean,rms), "l");

  for(std::vector<Cut>::size_type i = 0; i != QuantCuts.size(); i++) {
    TH1D* hist = (TH1D*) fout.Get(Form("rmtot_q%d",(i+1)));    
    hist -> SetLineColor(i+2);
    if((i+2) == 5) hist -> SetLineColor(6);
    hist -> Draw("histsame");
    double rms = hist -> GetRMS();
    double mean = hist -> GetMean();
    leg -> AddEntry(hist, Form("Q%d",(i+1)), "l");
    leg -> AddEntry(hist, Form("Mean:%.03lf, RMS:%.03lf",mean,rms), "l");
  }
  leg -> SetFillStyle(0);
  leg -> SetFillColor(0);
  leg -> Draw();
  c -> SaveAs((horn+"_energyResSpec.pdf").c_str());

  fout.Close();

  

}