numu_cut_flow.C

Go to the documentation of this file.

// Fills spectra for systematic error band
#ifdef __CINT__
void numu_cut_flow()
{
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else



#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/NumuCuts.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/NumuSysts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Vars/HistAxes.h"
#include "CAFAna/Core/OscillatableSpectrum.h"
#include "CAFAna/Core/Spectrum.h"
#include "OscLib/IOscCalc.h"
#include "OscLib/OscCalcDumb.h"


#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.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 <cmath>
#include <iostream>
#include <vector>
#include <list>
#include <sstream>
#include <string>
#include <fstream>
#include <iomanip>
#include <stdexcept>


using namespace ana;

class CanMan {

public:
  CanMan():
  fCans()
  {}

  TCanvas* NewCan(std::string name, std::string title)
  {
    fCans.push_back(new TCanvas(name.c_str(), title.c_str()));
    return fCans.back();
  }

  TCanvas* NewCan(std::string nameTitle)
  {
    return NewCan(nameTitle, nameTitle);
  }

  TCanvas* SaveCans(std::string prefix="", std::string suffix=".pdf")
  {
    for(const auto& can:fCans)
    {
      std::string name = prefix + std::string(can->GetName()) + suffix;
      can->SaveAs(name.c_str());
    }
  }

  std::vector<TCanvas*> fCans;

};


void addSame(std::string& opt)
{
  if (opt.find("same") == std::string::npos &&
      opt.find("SAME") == std::string::npos)
    opt += " SAME ";

}


class Sample
{
public:

  Sample(std::string shortName, std::string longName, Cut cut)
  :
  fShortName(shortName),
  fLongName(longName),
  fCut(cut)
  {};

  std::string fShortName;
  std::string fLongName;
  Cut         fCut;

};

class CutLevel
{
public:

  CutLevel(std::string shortName, std::string longName, Cut cut,
           SpectrumLoader& loader, const HistAxis& axis, int color)
  :
  fSpec(loader, axis, cut),
  fCut(cut),
  fShortName(shortName),
  fLongName(longName),
  fColor(color)
  {}

  void Draw(float pot, std::string opt,
            osc::IOscCalc* calc, int from, int to, TLegend* leg,
            int style=1)
  {
    TH1* h = fSpec.Oscillated(calc, from, to).ToTH1(pot);
    h->SetLineColor(fColor);
    h->SetLineStyle(style);
    h->Draw(opt.c_str());
    leg->AddEntry(h, fLongName.c_str(), "l");

  };

  void DrawEfficiency(float pot, std::string opt,
            osc::IOscCalc* calc, int from, int to, TH1* denominator,
            TLegend* leg)
  {
    TH1* h = fSpec.Oscillated(calc, from, to).ToTH1(pot);
    h->GetYaxis()->SetTitle("Efficiency");
    h->Divide(denominator);
    h->SetLineColor(fColor);
    h->SetLineStyle(fStyle);
    h->Draw(opt.c_str());
    leg->AddEntry(h, fLongName.c_str(), "l");
  };

  TH1* ToTH1(float pot, osc::IOscCalc* calc, int from, int to)
  {
    return fSpec.Oscillated(calc, from, to).ToTH1(pot);
  }



  OscillatableSpectrum fSpec;
  Cut fCut;
  std::string fShortName;
  std::string fLongName;
  int fColor;
  int fStyle;
};


class CutFlow
{
public:
  CutFlow(std::string name, SpectrumLoader& loader, const HistAxis& axis,
          std::vector<Sample> samples):
  fName(name),
  fLoader(loader),
  fAxis(axis)
  {

    fSamples.reserve(samples.size());
    for(const auto& sample:samples)
    {
     fSamples.emplace_back(sample.fShortName, sample.fLongName,
                                    sample.fCut,
                                    fLoader, fAxis, kBlack);
    }
    fLevels.resize(fSamples.size());
  };

  void AddLevel(std::string shortName, std::string longName, Cut cut, int color)
  {
    for(size_t iSample = 0; iSample < fSamples.size(); ++iSample)
    {
      fLevels[iSample].emplace_back(shortName, longName,
                                    fSamples[iSample].fCut && cut,
                                    fLoader, fAxis, color);
    }

  };


  void Reserve(size_t n)
  {

    for(auto& samp:fLevels)
      samp.reserve(n);
  };


  int findSample(std::string sampleName)
  {
    int i = -1;
    for(size_t iSample = 0; iSample < fSamples.size(); ++iSample){
      if(!sampleName.compare(fSamples[iSample].fShortName))
      {
        i = iSample;
      }
    }
    if(i == -1) throw std::invalid_argument("Sample not found");
    return i;
  }

  void Draw(float pot, std::string sample, std::string opt,
            osc::IOscCalc* calc, int from, int to, int style=1)
  {
    int i = findSample(sample);

    std::string name = fName + "_" + fSamples[i].fShortName;
    std::string currOpt(opt);
    TLegend* leg = new TLegend(0.55, 0.55, 0.85, 0.85);
    fSamples[i].Draw(pot, currOpt, calc, from, to, leg, style);
    addSame(currOpt);
    for(size_t iLevel = 0; iLevel < fLevels[i].size(); ++iLevel){
      fLevels[i][iLevel].Draw(pot, currOpt, calc, from, to, leg, style);
    }
    leg->Draw();

  };


  void DrawEfficiency(float pot, std::string sample, std::string opt,
            osc::IOscCalc* calc, int from, int to)
  {

    int i = findSample(sample);
    std::string name = fName + "_" + fSamples[i].fShortName;
    std::string currOpt(opt);
    TLegend* leg = AutoPlaceLegend(0.3, 0.23, 0.33);

    TH1* denominator = fSamples[i].ToTH1(pot, calc, from, to);

    for(size_t iLevel = 0; iLevel < fLevels[i].size(); ++iLevel){
      fLevels[i][iLevel].DrawEfficiency(pot, currOpt,
                                        calc, from, to, denominator, leg);
      addSame(currOpt);
    }
    leg->Draw();

  };



  std::string fName;
  SpectrumLoader& fLoader;
  HistAxis fAxis;
  std::vector<CutLevel> fSamples;

  std::vector<std::vector<CutLevel> > fLevels;


};



void numu_cut_flow()
{


  Cut truthContain({"sel.contain.missE", "mc.nu.E", "mc.nnu"},
                    [](const caf::StandardRecord* sr)
                    {
                      if(sr->mc.nnu == 0) return false;
                      return sr->sel.contain.missE / sr->mc.nu[0].E < 0.01;
                    });

  Cut cvn({"sel.cvnnumu.output", "sel.cvnnumu.noutput"},
                    [](const caf::StandardRecord* sr)
                    {
                      if(!sr->sel.cvnnumu.noutput) return false;
                      return sr->sel.cvnnumu.output[0] > 0.35;
                    });

  Cut cvn2({"sel.cvnnumu.output", "sel.cvnnumu.noutput"},
                    [](const caf::StandardRecord* sr)
                    {
                      if(!sr->sel.cvnnumu.noutput) return false;
                      return sr->sel.cvnnumu.output[0] > 0.5;
                    });


  Var cvnOut({"sel.cvnnumu.output", "sel.cvnnumu.noutput"},
                    [](const caf::StandardRecord* sr)
                    {
                      if(!sr->sel.cvnnumu.noutput) return -5.f;
                      return sr->sel.cvnnumu.output[0];
                    });

  HistAxis cvnAxis("; CVN #nu_{#mu};Events", Binning::Simple(40, 0,1), cvnOut);

  SpectrumLoader nonswap(
    "/cloud/login/rocco/development/results/2015-10-20_caf_nonswap_numu/fardet_genie_fhc_nonswap_ideal_3000_r01000001_s*");
  nonswap.SetSpillCut(kStandardSpillCuts);
  CutFlow remid("remid", nonswap,
                  kNumuCCAxis,
                  {Sample("sig", "Truth Containment",
                           kIsNumuCC && truthContain),
                   Sample("bkg", "Truth Containment",
                            !kIsNumuCC && truthContain)});


  remid.Reserve(10);


  remid.AddLevel("quality", "Quality", truthContain && kNumuQuality,
                   kGreen+2);
  remid.AddLevel("containment", "Containment",
                   truthContain && kNumuQuality && kNumuContainFD, kBlue);
  remid.AddLevel("NCRej", "CVN Loose",
                   truthContain && kNumuQuality && kNumuContainFD && cvn,
                   kMagenta-3);

  remid.AddLevel("NCRej", "CVN Tight",
                   truthContain && kNumuQuality && kNumuContainFD && cvn2,
                   kCyan-3);
  remid.AddLevel("NCRej", "ReMId",
                   truthContain && kNumuQuality && kNumuContainFD && kNumuNCRej,
                   kRed);


  OscillatableSpectrum cvnSig(nonswap, cvnAxis, kIsNumuCC && truthContain);
  OscillatableSpectrum cvnBkg(nonswap, cvnAxis, !kIsNumuCC && truthContain);



  nonswap.Go();

  CanMan canMan;

  osc::IOscCalc* calc = new osc::NoOscillations();
  canMan.NewCan("sig");
  remid.Draw(3e20, "sig","H", calc, 14, 14);

  canMan.NewCan("sig_eff");
  remid.DrawEfficiency(3e20, "sig","H", calc, 14, 14);

  canMan.NewCan("bkg");
  remid.Draw(3e20, "bkg","H", calc, 14, 14);

  canMan.NewCan("bkg_eff");
  remid.DrawEfficiency(3e20, "bkg","H", calc, 14, 14);

  canMan.NewCan("cvn_dist");
  TH1* hSig = cvnSig.Oscillated(calc, 14, 14).ToTH1(3e20);
  TH1* hBkg = cvnBkg.Oscillated(calc, 14, 14).ToTH1(3e20);
  hSig->SetLineColor(kRed);
  hBkg->SetLineColor(kBlue);
  hSig->Draw("HIST ");
  hBkg->Draw("HIST  SAME");

  TLegend* leg = AutoPlaceLegend(0.35, 0.15);
  leg->AddEntry(hSig, "Contained #nu_{#mu} Signal", "l");
  leg->AddEntry(hBkg, "Contained Background", "l");
  leg->Draw();



  delete calc;

  osc::IOscCalc* calcOsc = new osc::OscCalcDumb();
  canMan.NewCan("sig_osc");
  remid.Draw(3e20, "sig","H", calcOsc, 14, 14);

  canMan.NewCan("sig_eff_osc");
  remid.DrawEfficiency(3e20, "sig","H", calcOsc, 14, 14);

  canMan.NewCan("bkg_osc");
  remid.Draw(3e20, "bkg","H", calcOsc, 14, 14);

  canMan.NewCan("bkg_eff_osc");
  remid.DrawEfficiency(3e20, "bkg","H", calcOsc, 14, 14);

  canMan.NewCan("cvn_dist_osc");
  TH1* hSigOsc = cvnSig.Oscillated(calcOsc, 14, 14).ToTH1(3e20);
  TH1* hBkgOsc = cvnBkg.Oscillated(calcOsc, 14, 14).ToTH1(3e20);
  hSigOsc->SetLineColor(kRed);
  hBkgOsc->SetLineColor(kBlue);
  hBkgOsc->Draw("HIST");
  hSigOsc->Draw("HIST  SAME ");

  TLegend* legOsc = AutoPlaceLegend(0.35, 0.15);
  legOsc->AddEntry(hSigOsc, "Contained #nu_{#mu} Signal", "l");
  legOsc->AddEntry(hBkgOsc, "Contained Background", "l");
  legOsc->Draw();

  canMan.SaveCans("cvn_plots/", ".pdf");

}


#endif