nue_nd_datastability.C

Go to the documentation of this file.

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/LoadFromFile.h"

#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"

#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/HistAxes.h"

#include "TCanvas.h"
#include "TH2.h"
#include "TFile.h"
#include "TLegend.h"
#include "TPad.h"
#include "TLine.h"

#include <vector>
#include <string>

using namespace ana;

const Var kRun = SIMPLEVAR(hdr.run);

const Var kEPH([](const caf::SRProxy* sr){
                 if (sr->slc.nhit == 0) return float(-1);
                 return float(sr->slc.calE/sr->slc.nhit);});
const Var kHadCalE([](const caf::SRProxy* sr){
    if (!sr->vtx.elastic.IsValid) return float(-1);
    if (sr->vtx.elastic.fuzzyk.nshwlid < 1) return float(sr->slc.calE);
    return float(sr->slc.calE -
                 sr->vtx.elastic.fuzzyk.png[0].shwlid.calE);
  });
const Var kShwLIDE([](const caf::SRProxy* sr){
    if ( !sr->vtx.elastic.IsValid ) return float(-1);
    if (sr->vtx.elastic.fuzzyk.nshwlid < 1) return float(-1);
    return float(sr->vtx.elastic.fuzzyk.png[0].shwlid.calE);
  });
const Var kShwLIDEPH([](const caf::SRProxy* sr){
    if (!sr->vtx.elastic.IsValid ) return float(-1);
    if (sr->vtx.elastic.fuzzyk.nshwlid < 1) return float(-1);
    float cale = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
    int   nhit = sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;
    float eph = (nhit > 0) ? cale/nhit : -1;
    return eph;
  });





std::vector<Spectrum*> MakeCutFlowSpectra(SpectrumLoader *l,Var wei)
{
  std::vector<Spectrum*> ret;
  HistAxis axis("CalE",Binning::Simple(10,0,5),kCaloE);
  //HistAxis axis = kNueSAEnergyAxis;


  std::vector<Cut> cuts = {kNoCut};
  cuts.push_back(kNueDQ2017CVN);
  cuts.push_back(cuts.back() && kNue2017NDFiducial);
  cuts.push_back(cuts.back() && kNue2017NDContain && kNue2017NDFrontPlanes);
  cuts.push_back(cuts.back() && kNue2017NDNHits &&
                 kNue2017NDEnergy && kNue2017NDProngLength);
  cuts.push_back(cuts.back() && kCVNe>0.75);

  for (Cut cut : cuts)
    ret.push_back(new Spectrum(*l,axis,cut,kNoShift,wei));

  return ret;
}

std::vector<HistAxis> InitAxes()
{
  std::vector<HistAxis> axes;

  axes.push_back(HistAxis("CalE",Binning::Simple(200,0,5),kCaloE));
  axes.push_back(HistAxis("NHit",Binning::Simple(200,0,200),kNHit));
  axes.push_back(HistAxis("NHit",Binning::Simple(200,0,0.05),kEPH));
  axes.push_back(HistAxis("ShwE",Binning::Simple(200,0,5),kShwLIDE));
  axes.push_back(HistAxis("ShwEPH",Binning::Simple(200,0,0.05),kShwLIDEPH));
  axes.push_back(HistAxis("HadE",Binning::Simple(200,0,2),kHadCalE));
  axes.push_back(HistAxis("pTp",Binning::Simple(100,0,1),kPtP));

  return axes;
}

std::vector<Spectrum*> MakeKinSpectra(SpectrumLoader *l, Var wei, Cut sel)
{
  std::vector<Spectrum*> ret;

  std::vector<HistAxis> axes = InitAxes();
  for (HistAxis axis : axes)
    ret.push_back(new Spectrum(*l,axis,sel,kNoShift,wei));

  return ret;
}

std::vector<Spectrum*> MakeKin2DSpectra(SpectrumLoader *l, Var wei, Cut sel)
{
  std::vector<Spectrum*> ret;

  HistAxis runaxis("Run",Binning::Simple(200,10300,12300),kRun);

  std::vector<HistAxis> axes = InitAxes();
  for (HistAxis axis : axes)
    ret.push_back(new Spectrum(*l,runaxis,axis,sel,kNoShift,wei));

  return ret;
}


// We want to measure the cutflow for events, and look at simple reco vars
// as a function of run number
void CreateAndFillSpectra(std::string fname)
{
  SpectrumLoader *lDa2(new SpectrumLoader("prod_caf_R16-11-12-feature_caf_size.b_nd_numi_fhc_epochs1-3c_v1_goodruns"));
  SpectrumLoader *lDa3(new SpectrumLoader("prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns"));

  SpectrumLoader *lMC2(new SpectrumLoader("prod_caf_R16-11-12-feature_caf_size.b_nd_genie_nonswap_genierw_fhc_nova_v08_full_v1"));
  SpectrumLoader *lMC3(new SpectrumLoader("prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1"));

  lDa2->SetSpillCut(kStandardSpillCuts);
  lDa3->SetSpillCut(kStandardSpillCuts);
  lMC2->SetSpillCut(kStandardSpillCuts);
  lMC3->SetSpillCut(kStandardSpillCuts);

  Var wei = kXSecCVWgt2017 * kPPFXFluxCVWgt;
  Cut kNueND = kNue2017NDPresel && kCVNe > 0.75;

  std::vector<Spectrum*> mc_cutflow2 = MakeCutFlowSpectra(lMC2,wei);
  std::vector<Spectrum*> mc_cutflow3 = MakeCutFlowSpectra(lMC3,wei);
  std::vector<Spectrum*> da_cutflow2 = MakeCutFlowSpectra(lDa2,wei);
  std::vector<Spectrum*> da_cutflow3 = MakeCutFlowSpectra(lDa3,wei);

  std::vector<Spectrum*> mc_kin2 = MakeKinSpectra(lMC2,wei,kNueND);
  std::vector<Spectrum*> mc_kin3 = MakeKinSpectra(lMC3,wei,kNueND);
  std::vector<Spectrum*> da_kin2 = MakeKinSpectra(lDa2,wei,kNueND);
  std::vector<Spectrum*> da_kin3 = MakeKinSpectra(lDa3,wei,kNueND);

  std::vector<Spectrum*> mc_2D2 = MakeKin2DSpectra(lMC2,wei,kNueND);
  std::vector<Spectrum*> mc_2D3 = MakeKin2DSpectra(lMC3,wei,kNueND);
  std::vector<Spectrum*> da_2D2 = MakeKin2DSpectra(lDa2,wei,kNueND);
  std::vector<Spectrum*> da_2D3 = MakeKin2DSpectra(lDa3,wei,kNueND);

  lDa2->Go();
  lDa3->Go();
  lMC2->Go();
  lMC3->Go();

  std::cout << "fname " << fname << std::endl;
  TFile *out = new TFile(fname.c_str(),"recreate");

  std::vector<std::string> cutname = {"nocut","dq","fiducial","contain",
                                       "presel","ptp","cvn"};
  assert(cutname.size() == da_cutflow2.size());
  for (int cIdx = 0; cIdx < int(cutname.size()); cIdx++){
    mc_cutflow2[cIdx]->SaveTo(out, ("prod2_mc_"+cutname[cIdx]).c_str());
    mc_cutflow3[cIdx]->SaveTo(out, ("prod3_mc_"+cutname[cIdx]).c_str());

    da_cutflow2[cIdx]->SaveTo(out, ("prod2_da_"+cutname[cIdx]).c_str());
    da_cutflow3[cIdx]->SaveTo(out, ("prod3_da_"+cutname[cIdx]).c_str());
  }

  std::vector<std::string> varname = {"calE","NHit","EPH","shwE",
                                       "shwEPH","HadE","pTp"};
  assert(varname.size() == da_kin2.size());
  for (int kIdx = 0; kIdx < int(varname.size()); kIdx++){
    mc_kin2[kIdx]->SaveTo(out, ("prod2_mc_"+varname[kIdx]).c_str());
    mc_kin3[kIdx]->SaveTo(out, ("prod3_mc_"+varname[kIdx]).c_str());
    mc_2D2[kIdx]->SaveTo(out, ("prod2_mc_runvs"+varname[kIdx]).c_str());
    mc_2D3[kIdx]->SaveTo(out, ("prod3_mc_runvs"+varname[kIdx]).c_str());

    da_kin2[kIdx]->SaveTo(out, ("prod2_da_"+varname[kIdx]).c_str());
    da_kin3[kIdx]->SaveTo(out, ("prod3_da_"+varname[kIdx]).c_str());
    da_2D2[kIdx]->SaveTo(out, ("prod2_da_runvs"+varname[kIdx]).c_str());
    da_2D3[kIdx]->SaveTo(out, ("prod3_da_runvs"+varname[kIdx]).c_str());
  }
}

void CenterTitles(TH1 *h)
{
  h->GetXaxis()->CenterTitle();
  h->GetYaxis()->CenterTitle();
}

void MakePlot(std::string fname, std::string name, std::string xtitle)
{
  TCanvas *can = new TCanvas();

  Spectrum *da_prod2=LoadFromFile<Spectrum>(fname,"prod2_da_"+name).release();
  Spectrum *da_prod3=LoadFromFile<Spectrum>(fname,"prod3_da_"+name).release();
  double da_pot2 = da_prod2->POT();
  double da_pot3 = da_prod3->POT();
  Spectrum *mc_prod2=LoadFromFile<Spectrum>(fname,"prod2_mc_"+name).release();
  Spectrum *mc_prod3=LoadFromFile<Spectrum>(fname,"prod3_mc_"+name).release();
  double mc_pot2 = mc_prod2->POT();
  double mc_pot3 = mc_prod3->POT();

  TH1 *da2 = da_prod2->ToTH1(da_pot3);
  TH1 *da3 = da_prod3->ToTH1(da_pot3);
  TH1 *mc2 = mc_prod2->ToTH1(da_pot3);
  TH1 *mc3 = mc_prod3->ToTH1(da_pot3);

  TH1 *rat_da2 = da_prod2->ToTH1(da_pot3);
  TH1 *rat_da3 = da_prod3->ToTH1(da_pot3);
  TH1 *rat_mc3 = mc_prod3->ToTH1(da_pot3);

  int nrebin = 1;
  if (name == "EPH" || name == "shwEPH") nrebin = 10;
  if (name == "calE" || name == "shwE" ||
      name == "HadE" || name == "NHit") nrebin = 20;
  da2->Rebin(nrebin);
  da3->Rebin(nrebin);
  mc2->Rebin(nrebin);
  mc3->Rebin(nrebin);
  rat_da2->Rebin(nrebin);
  rat_da3->Rebin(nrebin);
  rat_mc3->Rebin(nrebin);

  rat_da2->Divide(mc2);
  rat_da3->Divide(mc2);
  rat_mc3->Divide(mc2);

  TPad *pRat = new TPad(("rat"+name) .c_str(),"",0,0,1,0.35);
  TPad *pSpec= new TPad(("spec"+name).c_str(),"",0,0.35,1,1);

  pRat->Draw();
  pSpec->Draw();

  pSpec->cd();
  gPad->SetBottomMargin(0);
  mc2->SetLineColor(kBlue);
  mc3->SetLineColor(kRed);
  da2->SetLineColor(kGreen+2);
  da3->SetLineColor(kBlack);
  da2->SetMinimum(1e-3);

  double max = std::max(da2->GetMaximum(),da3->GetMaximum());
  max = std::max(max,mc2->GetMaximum());
  max = std::max(max,mc3->GetMaximum());

  da2->SetMaximum(1.2*max);
  da2->GetXaxis()->SetTitle("");
  da2->GetYaxis()->SetTitle("Events");
  da2->SetTitle("");
  CenterTitles(da2);
  da2->GetXaxis()->SetTitleSize(0.06);
  da2->Draw("E");
  da3->Draw("E,same");
  mc2->Draw("E,same");
  mc3->Draw("E,same");
  TLegend *leg = new TLegend(0.6,0.65,0.85,0.85);
  leg->AddEntry(da2, "Data prod2", "le");
  leg->AddEntry(mc2, "MC prod2", "le");
  leg->AddEntry(da3, "Data prod3", "le");
  leg->AddEntry(mc3, "MC prod3", "le");
  leg->Draw();

  rat_mc3->SetLineColor(kRed);
  rat_da2->SetLineColor(kGreen+2);
  rat_da3->SetLineColor(kBlack);

  pRat->cd();
  gPad->SetTopMargin(0);
  gPad->SetBottomMargin(0.2);
  rat_mc3->GetXaxis()->SetTitle(xtitle.c_str());
  rat_mc3->GetYaxis()->SetTitle("Various Ratios");
  rat_mc3->SetTitle("");
  CenterTitles(rat_mc3);
  rat_mc3->GetXaxis()->SetTitleSize(0.10);
  rat_mc3->GetXaxis()->SetLabelSize(0.08);
  rat_mc3->GetXaxis()->SetTitleOffset(0.8);
  rat_mc3->GetYaxis()->SetTitleSize(0.10);
  rat_mc3->GetYaxis()->SetLabelSize(0.08);
  rat_mc3->GetYaxis()->SetTitleOffset(0.5);
  rat_mc3->GetYaxis()->SetRangeUser(0.5001,1.4999);
  rat_mc3->Draw();
  rat_da2->Draw("same");
  rat_da3->Draw("same");

  can->SaveAs(("Fig_nddatastability_"+name+".png").c_str());
}

void Make2DPlot(std::string fname, std::string name, std::string ytitle)
{
  TCanvas *can = new TCanvas();

  Spectrum *mc_prod2 = LoadFromFile<Spectrum>(fname,"prod2_mc_runvs"+name).release();
  Spectrum *mc_prod3 = LoadFromFile<Spectrum>(fname,"prod3_mc_runvs"+name).release();
  Spectrum *da_prod2 = LoadFromFile<Spectrum>(fname,"prod2_da_runvs"+name).release();
  Spectrum *da_prod3 = LoadFromFile<Spectrum>(fname,"prod3_da_runvs"+name).release();
  double pot = da_prod3->POT();

  TH2F *mc2 = (TH2F*)mc_prod2->ToTH2(pot);
  TH2F *mc3 = (TH2F*)mc_prod3->ToTH2(pot);
  TH2F *da2 = (TH2F*)da_prod2->ToTH2(pot);
  TH2F *da3 = (TH2F*)da_prod3->ToTH2(pot);

  int nrebin = 2;
  if (name == "HadE" || name == "shwE")
    nrebin = 4;
  mc2->RebinX(nrebin);
  mc3->RebinX(nrebin);
  da2->RebinX(nrebin);
  da3->RebinX(nrebin);

  TH1D *pmc2 = (TH1D*)mc2->ProfileX((name+"pmc2").c_str());
  TH1D *pmc3 = (TH1D*)mc3->ProfileX((name+"pmc3").c_str());
  TH1D *pda2 = (TH1D*)da2->ProfileX((name+"pda2").c_str());
  TH1D *pda3 = (TH1D*)da3->ProfileX((name+"pda3").c_str());

  pmc2->SetLineColor(kBlue);
  pmc3->SetLineColor(kRed);
  pda2->SetLineColor(kGreen+2);
  pda3->SetLineColor(kBlack);

  pda2->GetXaxis()->SetTitle("Run");
  pda2->GetYaxis()->SetTitle(ytitle.c_str());
  CenterTitles(pda2);

  int nnonzero = 0;
  double avgval = 0;
  for (int i = 1; i <= pda3->GetNbinsX(); i++){
    double curcont = pda3->GetBinContent(i);
    if (curcont != 0){
      nnonzero++;
      avgval += curcont;
    }
  }
  avgval /= nnonzero;
  pda2->GetYaxis()->SetRangeUser(0.5*avgval,1.5*avgval);

  pda2->Draw();
  pda3->Draw("same");
  pmc3->Draw("same");
  pmc2->Draw("same");

  TLegend *leg = new TLegend(0.6,0.65,0.85,0.85);
  leg->AddEntry(pda2, "Data prod2", "le");
  leg->AddEntry(pmc2, "MC prod2", "le");
  leg->AddEntry(pda3, "Data prod3", "le");
  leg->AddEntry(pmc3, "MC prod3", "le");
  leg->Draw();


  can->SaveAs(("Fig_nddatastability_"+name+"_vsrun.png").c_str());
}


void MakePlots(std::string fname)
{
  Make2DPlot(fname,"calE","CalE [GeV]");
  Make2DPlot(fname,"NHit","NHit");
  Make2DPlot(fname,"EPH","CalE / NHit [GeV]");
  Make2DPlot(fname,"HadE","HadE [GeV]");
  Make2DPlot(fname,"shwE","ShwE [GeV]");
  Make2DPlot(fname,"shwEPH","ShwE / ShwNHit [GeV]");



  MakePlot(fname,"calE","CalE [GeV]");
  MakePlot(fname,"NHit","NHit");
  MakePlot(fname,"EPH","CalE / NHit [GeV]");
  MakePlot(fname,"HadE","HadE [GeV]");
  MakePlot(fname,"shwE","ShwE [GeV]");
  MakePlot(fname,"shwEPH","ShwE / ShwNHit [GeV]");

  MakePlot(fname,"nocut","CalE [GeV]");
  MakePlot(fname,"dq","CalE [GeV]");
  MakePlot(fname,"fiducial","CalE [GeV]");
  MakePlot(fname,"contain","CalE [GeV]");
  MakePlot(fname,"presel","CalE [GeV]");
  MakePlot(fname,"ptp","CalE [GeV]");
  MakePlot(fname,"cvn","CalE [GeV]");
}


void nd_datastability(bool fillSpectra, std::string fname)
{
  if (fillSpectra)
    CreateAndFillSpectra(fname);
  else
    MakePlots(fname);

}