efficiency.C

Go to the documentation of this file.

// Fills spectra for systematic error band

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Systs/Systs.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "3FlavorAna/Vars/Binnings.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Vars/HistAxes.h"

#include "OscLib/OscCalcPMNSOpt.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"
#include "TGaxis.h"
#include "TMultiGraph.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TStyle.h"
#include "THStack.h"
#include "TPaveText.h"
#include "TList.h"
#include "TGaxis.h"
#include "TAttLine.h"
#include "TAttMarker.h"

#include <cmath>
#include <iostream>
#include <vector>
#include <list>
#include <sstream>
#include <string>
#include <sstream>
#include <fstream>
#include <iomanip>

const double pot = 6e20;

using namespace ana;

void efficiency()
{

  std::string fd_nonswap = "prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1";

  Loaders loaders;
  loaders.SetLoaderPath( fd_nonswap,  caf::kFARDET,  Loaders::kMC,   ana::kBeam, Loaders::kNonSwap );

  loaders.SetSpillCut(kStandardSpillCuts);

  osc::OscCalcPMNSOpt calc;
  ana::ResetOscCalcToDefault(&calc);


  //// ------ Now the cuts ------ /////

  const Cut kAny(// Any mode
                  [](const caf::SRProxy* sr)
                  {
                    return (sr->mc.nnu > 0);
                  });
  const Cut kQE( // QE
                  [](const caf::SRProxy* sr)
                  {
                    return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 0;
                  });
  const Cut kRes( // Res
                  [](const caf::SRProxy* sr)
                  {
                    return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 1;
                  });
  const Cut kDIS( // DIS
                  [](const caf::SRProxy* sr)
                  {
                    return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 2;
                  });
  const Cut kCoh( // Coh
                  [](const caf::SRProxy* sr)
                  {
                    return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 3;
                  });
  const Cut kMEC( // MEC
                  [](const caf::SRProxy* sr)
                  {
                    return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 10;
                  });

  const Cut kTruthContainFD( // only approx
                            [](const caf::SRProxy* sr)
                            {
                              return sr->mc.nnu > 0 && sr->mc.nu[0].E > 0 && sr->sel.contain.missE/sr->mc.nu[0].E < 0.01;
                                });

  struct Modes
  {
    std::string fname;
    Cut cut;
  };

  const int kNumModes = 6;

  Modes modes[kNumModes] = {
    {"All", kAny},
    {"QE", kQE},
    {"Res", kRes},
    {"DIS", kDIS},
    {"Coh", kCoh},
    {"MEC", kMEC}
  };

  struct Plot
  {
    std::string label;
    std::string fname;
    Binning bins;
    Var var;
  };

  const int kNumPlots = 2;

  Plot plots[kNumPlots] = {
    {";Reconstructed Neutrino Energy (GeV);Efficiency", "kCCE", kNumuCCEOptimisedBinning, kCCE},
    {";True Neutrino Energy (GeV);Efficiency", "kTrueE", kNumuCCEOptimisedBinning, kTrueE}
  };

  struct Selection
  {
    std::string label;
    std::string fname;
    Cut cut;
  };

  const int kNumCuts = 5;


  const Cut kBasicCut = kIsNumuCC && kTruthContainFD;

  Selection cuts[kNumCuts] = {
    {"No cuts","NoCut", kBasicCut},
    {"Data quality","DQ", kBasicCut && kNumuQuality},
    {"+ containment","DQ_contain", kBasicCut && kNumuQuality && kNumuContainFD2017},
    {"+ PIDs","DQ_contain_NCRej", kBasicCut && kNumuQuality && kNumuContainFD2017 && kNumuPID2017},
    {"+ Cosmic rej.","NumuFD", kBasicCut && kNumuQuality && kNumuContainFD2017 && kNumuPID2017 && kNumuCosmicRej2017}
  };

  int colours[kNumCuts] = {kBlack, kBlue, kGreen+2, kOrange-3, kRed};

  PredictionNoExtrap* sSpect[kNumPlots][kNumCuts][kNumModes];

  for(int k = 0; k < kNumModes; k++)
  {
    Cut mode = modes[k].cut;

    for(int i = 0; i < kNumPlots; i++)
    {
      Plot p = plots[i];

      for(int j = 0; j < kNumCuts; j++)
      {
        Cut cut = cuts[j].cut;

        sSpect[i][j][k] = new PredictionNoExtrap(loaders, p.label, p.bins, p.var, cut && mode, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017 );
      } // loop over cuts
    } // loop over plots
  } // loop over modes

  loaders.Go();

  TH1* hHist[kNumPlots][kNumCuts][kNumModes];

  for(int k = 0; k < kNumModes; k++)
  {
    for(int i = 0; i < kNumPlots; i++)
    {
      for(int j = 0; j < kNumCuts; j++)
      {
        std::string histName = plots[i].fname+"_"+cuts[j].fname;

        // hHist[i][j][k] = sSpect[i][j][k]->Predict(&calc).ToTH1(pot);
        hHist[i][j][k] = sSpect[i][j][k]->PredictUnoscillated().ToTH1(pot);
        hHist[i][j][k]->SetName(histName.c_str());

        hHist[i][j][k]->SetMarkerColor(colours[j]);
        hHist[i][j][k]->SetLineColor(colours[j]);      
      } // end loop over cuts
    } // end loop over plots
  } // end loop over modes

  for(int k = 0; k < kNumModes; k++)
  {
    for(int i = 0; i < kNumPlots; i++)
    {
      std::string canvasName = plots[i].fname+"_"+modes[k].fname;
      TCanvas* c1 = new TCanvas("",canvasName.c_str());

      TLegend* leg = new TLegend(0.65,0.18,0.85,0.40);
      leg->SetFillColor(10);
      leg->SetLineWidth(0);
      leg->SetLineColor(10);

      for (int j = 1; j < kNumCuts; j++)
      {
        hHist[i][j][k]->Divide(hHist[i][0][k]);

            std::cout << modes[k].fname.c_str() << " "
                  << plots[i].fname.c_str() << " " 
                  << cuts[j].fname.c_str() << " : "
                  << hHist[i][j][k]->Integral() << std::endl;

        if(j == 1) hHist[i][j][k]->Draw("hist");
        else hHist[i][j][k]->Draw("hist same");
        leg->AddEntry(hHist[i][j][k], cuts[j].label.c_str(), "l");
      } // loop over cuts

      leg->Draw();
      c1->SaveAs(("eff/"+canvasName+".pdf").c_str());
      c1->SaveAs(("eff/"+canvasName+".root").c_str());
    } // loop over plots
  } // loop over modes

}