nue_pid_effs_paper_numu_neweff.C

Go to the documentation of this file.

#ifdef __CINT__
void nue_pid_effs_paper_numu_neweff()
{
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcDumb.h"

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "CAFAna/Cuts/NumuCuts.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Vars/Vars.h"
#include "StandardRecord/Proxy/SRProxy.h"

#include "TArrow.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TGraph.h"
#include "TLegend.h"

#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>

#include "Utilities/func/MathUtil.h"

using namespace ana;
const double kPOTMACRO = 6.754e20;

// 250MeV bins
const Binning kRecoEBinning = Binning::Simple(16, 0, 4);
const Binning kModeBinning  = Binning::Simple(4, -0.5, 3.5);
const Binning kYBinning     = Binning::Simple(50, 0, 1);

/***************************************************************************/

void Print(std::string prefix, std::string name, std::string suffix = "")
{
  name = "plots/assessPaperNuMu/"+prefix+"_"+name;
  if(!suffix.empty()) name += "_"+suffix;

  gPad->Print((name+".eps").c_str());
  gPad->Print((name+".pdf").c_str());
  gPad->Print((name+".png").c_str());
}

const Cut kNumuFAFDContain(
                           //{ "mc.nu.iscc", "mc.nu.pdg", "mc.nu.E",
                            //   "sel.contain.missE", "mc.nnu", "energy.numusimp.trkccE","sel.contain.numucontain","slc.nhit","sel.remid.pid","trk.ncosmic"},
                           [](const caf::SRProxy* sr)
                           {
                             if(sr->mc.nnu == 0) return false;
                             assert(sr->mc.nnu == 1);
                             return (
                                     (
                                      (sr->mc.nu[0].E < 5 && (sr->mc.nu[0].iscc==1 && abs(sr->mc.nu[0].pdg)==14)) || 
                                      (sr->mc.nu[0].iscc==0) || 
                                      (abs(sr->mc.nu[0].pdg)!=14)
                                      ) &&
                                     (
                                      ((sr->sel.contain.missE/sr->mc.nu[0].E) < 0.01 && (sr->mc.nu[0].iscc==1 && abs(sr->mc.nu[0].pdg)==14)) || 
                                      (sr->mc.nu[0].iscc==0) || 
                                      (abs(sr->mc.nu[0].pdg)!=14)
                                      ) &&
                                     sr->slc.nhit>=20 &&
                                     sr->sel.remid.pid>0 &&
                                     sr->trk.cosmic.ntracks>0 &&
                                     sr->mc.nnu>0 &&
                                     sr->energy.numu.trkccE < 5 &&
                                     sr->sel.contain.numucontain
                                     );
                           }
                           );


const Cut kNumuFAFDContainLight(
                                //{ "mc.nu.iscc", "mc.nu.pdg", "mc.nu.E",
                                //    "sel.contain.missE", "mc.nnu", "energy.numusimp.trkccE","sel.contain.numucontain","slc.nhit","sel.remid.pid","trk.ncosmic"},
                                [](const caf::SRProxy* sr)
                                {
                                  if(sr->mc.nnu == 0) return false;
                                  assert(sr->mc.nnu == 1);
                                  return ( 
                                          (
                                           (sr->mc.nu[0].E < 5 && (sr->mc.nu[0].iscc==1 && abs(sr->mc.nu[0].pdg)==14)) || 
                                           (sr->mc.nu[0].iscc==0) || 
                                           (abs(sr->mc.nu[0].pdg)!=14)
                                           ) &&
                                          (
                                           ((sr->sel.contain.missE/sr->mc.nu[0].E) < 0.01 && (sr->mc.nu[0].iscc==1 && abs(sr->mc.nu[0].pdg)==14)) || 
                                           (sr->mc.nu[0].iscc==0) || 
                                           (abs(sr->mc.nu[0].pdg)!=14)
                                           ) &&
                                          sr->mc.nnu>0 &&
                                    sr->energy.numu.trkccE < 5 
                                    );
                                }
  );

// True E, *y for NCs
const Var kTrueEVis(//{"mc.nu.E", "mc.nu.y", "mc.nu.iscc"},
                    [](const caf::SRProxy* sr)
                    {
                      if(sr->mc.nu.empty()) return 0.;
                      double ret = sr->mc.nu[0].E;
                      if(!sr->mc.nu[0].iscc) ret *= sr->mc.nu[0].y;
                      return ret;
                    });

/*
const Var kMode(//{"mc.nu.mode"},
                [](const caf::SRProxy* sr)
                {
                  int tmp = 0;
                  if(sr->mc.nu.empty()) tmp;
                  tmp = sr->mc.nu[0].mode;
                  return tmp;
                });

*/

const Var kY(//{"mc.nu.y"},
             [](const caf::SRProxy* sr)
             {
               float tmp = 0.f;
               if(sr->mc.nu.empty()) return tmp;
               tmp = sr->mc.nu[0].y;
               return tmp;
             });

/***************************************************************************/
// Nominal osc params
void resetCalc(osc::IOscCalcAdjustable* calc)
{
  calc->SetL(810);
  calc->SetRho(2.75);
  calc->SetDmsq21(7.6e-5);
  calc->SetDmsq32(2.35e-3);
  calc->SetTh12(asin(sqrt(.87))/2); // PDG
  calc->SetTh13(asin(sqrt(.095))/2);
  calc->SetTh23(TMath::Pi()/4);
  calc->SetdCP(0);
}

/***************************************************************************/

// Draw cuts on plots
void Arrows(TH1* sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow = false)
{
  sigId->GetYaxis()->SetRangeUser(0,1.05*sigId->GetMaximum());
  
  double maxy = sigId->GetMaximum()/2;
  if(noarrow)
    maxy = maxy*2;
  
  TLine* lin = new TLine(bestCutSoB, 0, bestCutSoB, maxy);
  lin->SetLineWidth(3);
  lin->SetLineColor(kGreen+2);
  lin->Draw();
  if(!noarrow){
    TArrow* arr = new TArrow(bestCutSoB-0.003, maxy, bestCutSoB+.03, maxy, .02, "|>");
    arr->SetLineWidth(3);
    arr->SetLineColor(kGreen+2);
    arr->SetFillColor(kGreen+2);
    arr->Draw();
    maxy *= 1.5;
  }



  lin = new TLine(bestCutSoSB, 0, bestCutSoSB, maxy);
  lin->SetLineWidth(3);
  lin->SetLineStyle(2);
  lin->SetLineColor(kGreen+3);
  lin->Draw();

  if(!noarrow){
    TArrow* arr = new TArrow(bestCutSoSB-0.003, maxy, bestCutSoSB+.04, maxy, .02, "|>");
    arr->SetLineWidth(3);
    arr->SetLineColor(kGreen+3);
    arr->SetFillColor(kGreen+3);
    arr->Draw();
  }


  lin = new TLine(faCut, 0, faCut, maxy);
  lin->SetLineWidth(3);
  lin->SetLineColor(kRed);
  //lin->Draw();

  if(!noarrow){
    TArrow* arr = new TArrow(faCut, maxy, faCut+.04, maxy, .02, "|>");
    arr->SetLineWidth(3);
    arr->SetLineColor(kRed);
    arr->SetFillColor(kRed);
    //arr->Draw();
  }
}

/***************************************************************************/
// Finds best cuts
void FOMPlot(TH1* sig, TH1* bkg, TH1* bkgCC, TH1* bkgBeam, 
             TH1* sig_denom, double* bestCuts, double faCut,
             std::string suffix, std::string eff_suffix)
{
  std::cout << suffix << " " << eff_suffix << ":" << std::endl;
  std::cout << " &\t Cut &\t $s/\\sqrt{b}$ &\t $s/\\sqrt{s+b}$ &\t $\\nu_\\mu$ sig &\t Tot bkg &\t NC &\t Appeared $\\nu_e$ &\t Beam $\\nu_e$ &\t Sig eff &\t Purity\\\\" 
            << std::endl;
  
  TH1F* eff     = new TH1F(("eff_"+suffix+"_"+eff_suffix).c_str(),
                           "Efficiency vs Cut Value; ;Efficiency (%)",
                           sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1), 
                           sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
  
  TH1F* effxpur     = new TH1F(("effxpur_"+suffix+"_"+eff_suffix).c_str(),
                               "Efficiency*Purity vs Cut Value; ;Efficiency*Purity (%)",
                               sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1), 
                               sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
  
  TH1F* pur     = new TH1F(("pur_"+suffix+"_"+eff_suffix).c_str(),
                           "Purity vs Cut Value; ;Purity (%)",
                           sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1), 
                           sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
  TH1F* fomsob  = new TH1F(("sob_"+suffix+"_"+eff_suffix).c_str(),
                           "s/#sqrt{b} vs Cut Value; ;s/#sqrt{b}",
                           sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1), 
                           sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
  TH1F* fomsosb = new TH1F(("sosb_"+suffix+"_"+eff_suffix).c_str(),
                           "s/#sqrt{s+b} vs Cut Value; ;s/#sqrt{s+b}",
                           sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1), 
                           sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));

  eff->GetXaxis()->CenterTitle();
  eff->GetYaxis()->CenterTitle();
  pur->GetXaxis()->CenterTitle();
  pur->GetYaxis()->CenterTitle();
  fomsob->GetXaxis()->CenterTitle();
  fomsob->GetYaxis()->CenterTitle();
  fomsosb->GetXaxis()->CenterTitle();
  fomsosb->GetYaxis()->CenterTitle();
  
  eff->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
  pur->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());

  fomsob->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
  fomsosb->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
  int ngpoints=0;
  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
    if(bkg->Integral(n, -1) == 0) break;
    ngpoints++;
  }
  TGraph* effpurGraph=new TGraph (ngpoints);
  for(int metric = 0; metric < 3; ++metric){
    double bestSoB   = -1;
    double bestSoSB  = -1;
    double bestSig   = 0;
    double bestBkg   = 0;
    double bestBkgCC = 0;
    double bestBeam  = 0;
    double totsig    = sig_denom->Integral(0, -1);
    double fomeff, fompur; 
        
    for(int n = 0; n < sig->GetNbinsX()+2; ++n){
      const double nsig = sig->Integral(n, -1);
      const double nbkg = bkg->Integral(n, -1);
      const double nbkgcc = bkgCC->Integral(n, -1);
      const double nbkgbeam = bkgBeam->Integral(n, -1);

      if(nbkg == 0) break;

      double neff = nsig*100/(totsig);
      double npur = nsig*100/(nsig+nbkg);

      const double fomSoB = nsig/sqrt(nbkg);
      const double fomSoSB = nsig/sqrt(nsig+nbkg);

      if(metric == 0){
        // only fill these hists once.
        pur->Fill(sig->GetXaxis()->GetBinCenter(n), npur);
        eff->Fill(sig->GetXaxis()->GetBinCenter(n), neff);
        effxpur->Fill(sig->GetXaxis()->GetBinCenter(n), 100*(neff/100)*(npur/100));
        effpurGraph->SetPoint(n,neff,npur);
        fomsob->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoB);
        fomsosb->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoSB);
      }

      if((metric == 0 && fomSoB > bestSoB) ||
         (metric == 1 && fomSoSB > bestSoSB) ||
         (metric == 2 && fabs(sig->GetXaxis()->GetBinLowEdge(n) - faCut) < .01)){
        bestSoB   = fomSoB; 
        bestSoSB  = fomSoSB;
        bestSig   = nsig;
        bestBkg   = nbkg;
        bestBkgCC = nbkgcc;
        bestBeam  = nbkgbeam;
        bestCuts[metric] = sig->GetXaxis()->GetBinLowEdge(n);
        fomeff    = neff;
        fompur    = npur;
      }
    }

    if(metric == 0) std::cout << "$s/\\sqrt{b}$ opt";
    if(metric == 1) std::cout << "$s/\\sqrt{s+b}$ opt";
    if(metric == 2) std::cout << "FA cut";

    std::cout << std::fixed << std::setprecision(2);
    std::cout << " &\t " << ((metric == 2) ? faCut : bestCuts[metric])
              << " &\t " << bestSoB << " &\t " << bestSoSB
              << " &\t " << bestSig << " &\t " << bestBkg 
              << " &\t " << bestBkg-bestBkgCC-bestBeam 
              << " &\t " << bestBkgCC << " &\t " << bestBeam 
              << " &\t " << fomeff << "\\% &\t " << fompur << "\\%\\\\"<< std::endl;

  } // end for metric

  new TCanvas;
  eff->Draw("l");
  Arrows(eff, bestCuts[0], bestCuts[1], faCut, true);
  Print("eff", suffix+"_"+eff_suffix);
  
  new TCanvas;
  pur->Draw("l");
  Arrows(pur, bestCuts[0], bestCuts[1], faCut, true);
  Print("pur", suffix);

  new TCanvas;
  eff->SetTitle("");
  eff->GetYaxis()->SetTitle("Percentage");
  eff->SetLineColor(kRed);
  pur->SetLineColor(kBlue);
  effxpur->SetLineColor(kGreen+2);
  eff->SetMaximum(110);
  
  eff->Draw("c");
  pur->Draw("csame");
  effxpur->Draw("csame");
  //Arrows(pur, bestCuts[0], bestCuts[1], faCut, true);

  TLegend* legeffpur = new TLegend(.125, .12, .5, .32);
  legeffpur->AddEntry(eff, "Cumulative Efficiency", "l");
  legeffpur->AddEntry(pur, "Cumulative Purity", "l");
  legeffpur->AddEntry(effxpur, "Cumulative Efficiency #times Purity", "l");
  legeffpur->SetFillStyle(0);
  legeffpur->Draw();
  Print("effpur", suffix+"_"+eff_suffix);

  new TCanvas;
  TH2F* hpx = new TH2F("hpx","hpx",100,0,100,100,0,100); // axis range hpx->SetStats(kFALSE); // no statistics
  hpx->SetStats(kFALSE);
  hpx->SetTitle("");
  hpx->GetXaxis()->SetTitle("Efficiency (%)");
  hpx->GetYaxis()->SetTitle("Purity (%)");
  hpx->GetXaxis()->CenterTitle();
  hpx->GetYaxis()->CenterTitle();
  hpx->Draw();
  effpurGraph->SetLineWidth(2);
  effpurGraph->Draw("C");
  Print("effpurgraph", suffix+"_"+eff_suffix);


  new TCanvas;
  fomsob->Draw("l");
  Arrows(fomsob, bestCuts[0], bestCuts[1], faCut, true);
  Print("fomsob", suffix);

  new TCanvas;
  fomsosb->Draw("l");
  Arrows(fomsosb, bestCuts[0], bestCuts[1], faCut, true);
  Print("fomsosb", suffix);
}

/***************************************************************************/

void GetSpectra(IPrediction* pred, osc::IOscCalc* calc,
                TH1*& hSig, TH1*& hNC, TH1*& hCC, TH1*& hBeam, TH1*&hTotBkg)
{
  hSig = pred->PredictComponent(calc,
                                Flavors::kNuMuToNuMu,
                                Current::kCC,
                                Sign::kBoth).ToTH1(kPOTMACRO);
  hSig->SetLineColor(kNumuBackgroundColor);

  hNC = pred->PredictComponent(calc,
                               Flavors::kAll,
                               Current::kNC,
                               Sign::kBoth).ToTH1(kPOTMACRO);
  hNC->SetLineColor(kNCBackgroundColor);

  hCC = pred->PredictComponent(calc,
                               Flavors::kNuMuToNuE,
                               Current::kCC,
                               Sign::kBoth).ToTH1(kPOTMACRO);
  hCC->SetLineColor(kNueSignalColor);

  hBeam = pred->PredictComponent(calc,
                                 Flavors::kNuEToNuE,
                                 Current::kCC,
                                 Sign::kBoth).ToTH1(kPOTMACRO);
  hBeam->SetLineColor(kBeamNueBackgroundColor);

  hTotBkg = (TH1*)hNC->Clone(UniqueName().c_str());
  hTotBkg->Add(hCC);
  hTotBkg->Add(hBeam);

  hTotBkg->SetLineColor(kTotalMCColor);
}

/***************************************************************************/

void PlotSpectra(IPrediction* pred, osc::IOscCalc* calc, std::string title, double maxy = 0)
{
  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
  GetSpectra(pred, calc, hSig, hNC, hCC, hBeam, hTotBkg);

  new TCanvas;
  hSig->SetTitle("");//title.c_str());
  hSig->GetYaxis()->SetTitle("Events / 18#times10^{20} POT");
  //hSig->GetYaxis()->SetTitle("Arbitary Units");
  hSig->GetXaxis()->CenterTitle();
  hSig->GetYaxis()->CenterTitle();
  hSig->Draw("hist");
  if(maxy) hSig->GetYaxis()->SetRangeUser(0, maxy);
  hNC->Draw("hist same");
  hCC->Draw("hist same");
  hBeam->Draw("hist same");
  TLegend* leg2 = new TLegend(.125, .55, .5, .85);
  leg2->AddEntry(hCC, "Appeared #nu_{e}", "l");
  leg2->AddEntry(hSig, "Survived #nu_{#mu}", "l");
  leg2->AddEntry(hNC, "NC background", "l");
  leg2->AddEntry(hBeam, "Beam #nu_{e} background", "l");
  leg2->SetFillStyle(0);
  leg2->Draw();

  static bool once = true;
  if(once){
    once = false;
    TVirtualPad* bak = gPad;
    new TCanvas;
    TLegend* leg = new TLegend(.1, .1, .9, .9);
    leg->AddEntry(hSig, "#nu_{#mu} signal", "l");
    leg->AddEntry(hNC, "NC background", "l");
    leg->AddEntry(hCC, "Appeared #nu_{e} background", "l");
    leg->AddEntry(hBeam, "Beam #nu_{e} background", "l");
    leg->Draw();
    gPad->Print("plots/assessPaperNuMu/eff_leg.pdf");
    gPad->Print("plots/assessPaperNuMu/eff_leg.eps");
    gPad->Print("plots/assessPaperNuMu/eff_leg.png");
    bak->cd();
  }
}

void PrintEffs(TH1* h, const char* title)
{
  std::cout << title;
  for(int i = 0; i < h->GetNbinsX()+2; ++i){
    if(h->GetBinCenter(i) > .75 &&
       h->GetBinCenter(i) < 3.5){
      std::cout << " &\t " << h->GetBinContent(i);
    }
  }
  std::cout << "\\\\" << std::endl;
}

/***************************************************************************/

void PlotEffs(IPrediction* predNum, IPrediction* predDenom,
              osc::IOscCalc* calc, std::string title, 
              std::string label, std::string id, double maxy = 0)
{
  TH1 *hSigNum, *hNCNum, *hCCNum, *hBeamNum, *hTotBkgNum;
  GetSpectra(predNum, calc, hSigNum, hNCNum, hCCNum, hBeamNum, hTotBkgNum);

  TH1 *hSigDenom, *hNCDenom, *hCCDenom, *hBeamDenom, *hTotBkgDenom;
  GetSpectra(predDenom, calc, hSigDenom, hNCDenom, hCCDenom, hBeamDenom, hTotBkgDenom);

  hSigDenom->SetLineStyle(2);
  hNCDenom->SetLineStyle(2);
  hBeamDenom->SetLineStyle(2);
  hTotBkgDenom->SetLineStyle(2);
  hCCDenom->SetLineStyle(2);

  new TCanvas();

  if( label == "mode" ){
    hNCDenom->GetXaxis()->SetLabelSize(0.06);
    hNCDenom->GetXaxis()->SetBinLabel(1, "QE");
    hNCDenom->GetXaxis()->SetBinLabel(2, "Res");
    hNCDenom->GetXaxis()->SetBinLabel(3, "DIS");
    hNCDenom->GetXaxis()->SetBinLabel(4, "Coh");

    hSigNum->GetXaxis()->SetLabelSize(0.06);
    hSigNum->GetXaxis()->SetBinLabel(1, "QE");
    hSigNum->GetXaxis()->SetBinLabel(2, "Res");
    hSigNum->GetXaxis()->SetBinLabel(3, "DIS");
    hSigNum->GetXaxis()->SetBinLabel(4, "Coh");

  }

  hNCDenom->Draw("hist");
  hCCDenom->Draw("hist same");
  hBeamDenom->Draw("hist same");
  hSigDenom->Draw("hist same");
  Print(label, id, "nocut");

  new TCanvas();
  if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
  //hSigNum->GetYaxis()->SetTitle("Arbitary Units");
  hSigNum->GetYaxis()->CenterTitle();
  hSigNum->Draw("hist");
  hNCNum->Draw("hist same");
  hBeamNum->Draw("hist same");
  hCCNum->Draw("hist same");
  Print(label, id, "cut");

  new TCanvas;

  std::cout << "Total efficiency"
            << " sig: " << 100*hSigNum->Integral(0, -1)/hSigDenom->Integral(0, -1)
            << " NC: " << 100*hNCNum->Integral(0, -1)/hNCDenom->Integral(0, -1)
            << " CC: " << 100*hCCNum->Integral(0, -1)/hCCDenom->Integral(0, -1)
            << " Beam: " << 100*hBeamNum->Integral(0, -1)/hBeamDenom->Integral(0, -1) << std::endl;

  hSigNum->Divide(hSigDenom);
  hSigNum->Scale(100);
  hSigNum->SetTitle(title.c_str());
  hSigNum->GetYaxis()->SetTitle("Efficiency (%)");
  hSigNum->GetYaxis()->SetRangeUser(0, 100);
  //if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
  hSigNum->GetXaxis()->CenterTitle();
  hSigNum->GetYaxis()->CenterTitle();
  hSigNum->Draw("hist");

  PrintEffs(hSigNum, "$\\nu_e$ sig");

  hNCNum->Divide(hNCDenom);
  hNCNum->Scale(100);
  hNCNum->Draw("hist same");

  PrintEffs(hNCNum, "NC");

  hCCNum->Divide(hCCDenom);
  hCCNum->Scale(100);
  hCCNum->Draw("hist same");

  PrintEffs(hCCNum, "$\\nu_\\mu$ CC");

  hBeamNum->Divide(hBeamDenom);
  hBeamNum->Scale(100);
  hBeamNum->Draw("hist same");

  PrintEffs(hBeamNum, "beam $\\nu_e$");
}

/***************************************************************************/

std::string PIDLongName(int pidIdx, bool rhc)
{
  std::string ret;
  if(pidIdx == -1) ret = "No selection";
  if(pidIdx == 0) ret = "remID";
  if(pidIdx == 1) ret = "CVN#mu";

  if(rhc) ret += " (RHC)"; else ret += " (FHC)";

  return ret;
}

/***************************************************************************/

std::string PIDFileTag(int pidIdx, bool rhc)
{
  std::string ret;
  if(pidIdx == -1) ret = "nosel";
  if(pidIdx == 0) ret = "lem";
  if(pidIdx == 1) ret = "cvn";

  if(rhc) ret += "_rhc"; else ret += "_fhc";

  return ret;
}

/***************************************************************************/

void nue_pid_effs_paper_numu_neweff()
{

  osc::OscCalcPMNSOpt calc;
  calc.SetL(810);
  calc.SetRho(2.84);
  calc.SetDmsq21(7.53e-5);
  calc.SetDmsq32(2.44e-3);
  calc.SetTh12(asin(sqrt(0.846))/2);
  calc.SetTh13(asin(sqrt(.085))/2);
  calc.SetTh23(asin(sqrt(0.999))/2);
  calc.SetdCP(0);

  for(int rhc = false; rhc <= true; ++rhc){
    if(rhc) continue; // Don't exist yet for S13-06-18

    /*
    const std::string fnameMC = 
      TString::Format("defname: prod_caf_S15-05-22_fd_genie_fhc_nonswap_ideal with stride 1",
                      rhc ? "rhc" : "fhc").Data();
    
    const std::string fnameSwap = 
      TString::Format("defname: prod_caf_S15-05-22_fd_genie_fhc_fluxswap_ideal with stride 1",
                      rhc ? "rhc" : "fhc").Data();
    */

    //SpectrumLoader loader("defname: prod_caf_R17-03-01-prod3reco.g_fd_genie_nonswap_fhc_nova_v08_full_v1 with limit 3");
    //SpectrumLoader loaderSwap("defname: prod_caf_R17-03-01-prod3reco.g_fd_genie_fluxswap_fhc_nova_v08_full_v1 with limit 3");

    //SpectrumLoader loader("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/fardet_genie_nonswap_genierw_fhc_v08_1000_r00014457_s20_c000_S17-05-31_v1_20170330_113817_sim.caf.root");
    //SpectrumLoader loaderSwap("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/fardet_genie_fluxswap_genierw_fhc_v08_1000_r00014031_s02_c000_S17-05-31_v1_20170330_145648_sim.caf.root");

    SpectrumLoader loader("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/semfe_nonS_big.root");
    SpectrumLoader loaderSwap("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/semfe_S_big.root");

    //SpectrumLoader loader("/pnfs/nova/scratch/users/atsaris/cvnFS_eval_fullInt_afterColl_swap/00014*/*/*.root");
    //SpectrumLoader loaderSwap("/pnfs/nova/scratch/users/atsaris/cvnFS_eval_fullInt_afterColl_nonSwap/00014*/*/*.root");

    const int kNumPIDs = 2;

    IPrediction* preds[kNumPIDs];
    //kNumuFAFDContain
    preds[0] = new PredictionNoExtrap(loader, loaderSwap,
                                      "remID", Binning::Simple(100, 0.0, 1.0),
                                      kRemID, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);

    preds[1] = new PredictionNoExtrap(loader, loaderSwap,
                                      "#nu_{#mu} CC Classifier Output", Binning::Simple(100, 0.0, 1.0),
                                      //kANGNuMu, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);
                                      kCVNm, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);

    PredictionNoExtrap predDenomNoCutRecoE(loader, loaderSwap,
                                           "Reco E (GeV)", kRecoEBinning, 
                                           kCaloE, kNumuFAFDContainLight);//kNumuQuality && kNumuContainFD);

    IPrediction* predDenomPreselRecoE[kNumPIDs];
    predDenomPreselRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
                                                     "Reco E (GeV)",
                                                     kRecoEBinning, 
                                                     kCaloE,
                                                     kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);
    
    predDenomPreselRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
                                                     "Reco E (GeV)",
                                                     kRecoEBinning, 
                                                     kCaloE,
                                                     kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);

    PredictionNoExtrap predDenomNoCutTrueE(loader, loaderSwap,
                                           "True Visible E (GeV)",
                                           kRecoEBinning, 
                                           kTrueEVis, kNoCut);

    IPrediction* predDenomPreselTrueE[kNumPIDs];
    predDenomPreselTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
                                                     "True Visible E (GeV)",
                                                     kRecoEBinning, 
                                                     kTrueEVis,
                                                     kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);
    predDenomPreselTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
                                                     "True Visible E (GeV)",
                                                     kRecoEBinning, 
                                                     kTrueEVis,
                                                     kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);

    IPrediction* predsRecoE[kNumPIDs];
    IPrediction* predsTrueE[kNumPIDs];

    PredictionNoExtrap predDenomMode(loader, loaderSwap,
                                     "", kModeBinning, 
                                     kMode, kNoCut);

    IPrediction* predsMode[kNumPIDs];

    PredictionNoExtrap predDenomY(loader, loaderSwap,
                                  "", kYBinning, 
                                  kY, kNoCut);

    IPrediction* predsY[kNumPIDs];

    predsRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
                                           "Reco E (GeV)", kRecoEBinning,
                                           kCaloE, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaLEMOpt);

    predsRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
                                           "Reco E (GeV)", kRecoEBinning,
                                           //kCaloE, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaANGOpt);
                                           kCaloE, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);

    predsTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
                                           "True Visible E (GeV)", kRecoEBinning,
                                           kTrueEVis, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaLEMOpt);

    predsTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
                                           "True Visible E (GeV)", kRecoEBinning,
                                           //kTrueEVis, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaANGOpt);
                                           kTrueEVis, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);

    predsMode[0] = new PredictionNoExtrap(loader, loaderSwap,
                                           "", kModeBinning,
                                          kMode, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaLEMOpt);

    predsMode[1] = new PredictionNoExtrap(loader, loaderSwap,
                                           "", kModeBinning,
                                          //kMode, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaANGOpt);
                                           kMode, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);

    predsY[0] = new PredictionNoExtrap(loader, loaderSwap,
                                       "True y", kYBinning,
                                       kY, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaLEMOpt);
    
    predsY[1] = new PredictionNoExtrap(loader, loaderSwap,
                                       "True y", kYBinning,
                                       //kY, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaANGOpt);
                                       kY, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);

    loader.Go();
    loaderSwap.Go();

    for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx){
      TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
      GetSpectra(preds[pidIdx], &calc, hSig, hNC, hCC, hBeam, hTotBkg);

      double cuts[2];
      const double faCut = 0.75;
      
      TH1* backoscnueNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
                                                                Flavors::kNuMuToNuE,
                                                                Current::kCC,
                                                                Sign::kBoth).ToTH1(kPOTMACRO);
      
      TH1* sigdenomNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
                                                                Flavors::kNuMuToNuMu,
                                                                Current::kCC,
                                                                Sign::kBoth).ToTH1(kPOTMACRO);

      TH1* backnueNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
                                                                Flavors::kNuEToNuE,
                                                                Current::kCC,
                                                                Sign::kBoth).ToTH1(kPOTMACRO);

      TH1* backncNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
                                                              Flavors::kAll,
                                                                Current::kNC,
                                                                Sign::kBoth).ToTH1(kPOTMACRO);


      double totalBckOscNuECC=backoscnueNoCut->Integral(0, -1);
      double totalSig=sigdenomNoCut->Integral(0, -1);
      double totalBckNuECC=backnueNoCut->Integral(0, -1);
      double totalBckNC=backncNoCut->Integral(0, -1);
      double totalBackground=totalBckOscNuECC+totalBckNuECC+totalBckNC;
      double totalsob=totalSig/sqrt(totalBackground);
      double totalsosb=totalSig/sqrt(totalSig+totalBackground);
      double totalefficiency=100;
      double totalpurity=(totalSig*100)/(totalSig+totalBackground);

      std::cout << std::fixed << std::setprecision(2);
      std::cout << " &\t " << totalsob << " &\t " << totalsosb
                << " &\t " << totalSig << " &\t " << totalBackground
                << " &\t " << totalBckNC
                << " &\t " << totalBckOscNuECC << " &\t " << totalBckNuECC
                << " &\t " << totalefficiency << "\\% &\t " << totalpurity << "\\%\\\\"<< std::endl;
      // & 2.47 & 2.14 & 18.31 & 54.93 & 18.31 & 18.31  18.31 & 100.00\% & 75.00\%\\

      sigdenomNoCut->SetLineColor(kNueSignalColor);

      TH1* sigdenomPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
                                                                           Flavors::kNuMuToNuMu,
                                                                           Current::kCC,
                                                                           Sign::kBoth).ToTH1(kPOTMACRO);
      sigdenomPresel->SetLineColor(kNueSignalColor);
      
      FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomNoCut, cuts, faCut, PIDFileTag(pidIdx, rhc), "nocut");
      //FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomPresel, cuts, faCut, PIDFileTag(pidIdx, rhc), "presel");

      PlotSpectra(preds[pidIdx], &calc, PIDLongName(pidIdx, rhc));//,16.0);

      Arrows(hSig, cuts[0], cuts[1], faCut);
      Print("pid", PIDFileTag(pidIdx, rhc));
    } // end for pidIdx

    PlotSpectra(&predDenomNoCutRecoE, &calc, PIDLongName(-1, rhc), 40);
    Print("recoE", PIDFileTag(-1, rhc));
    PlotSpectra(&predDenomNoCutTrueE, &calc, PIDLongName(-1, rhc), 40);
    Print("trueE", PIDFileTag(-1, rhc));

    return;
  } // end for rhc
}

#endif