datamc.C

Go to the documentation of this file.

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

// analysis cuts/vars
#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncExtra.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Cuts/SpillCuts.h"

#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "TFile.h"
#include "THStack.h"
#include "TH1D.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TPad.h"
#include "TGraph.h"
#include "TArrow.h"


// temporary measure to use prod4 trained muonid
#include "NDAna/muonid/NDXSecMuonPID.h"
namespace nuebarccinc {
  const ana::Cut kMuonIDProd4Cut = ana::kMuonIDProd4 < -0.2;
  void One();
}

void load_libs_muonid()
{
  const int ret = gSystem->Load("libNDAnamuonid");
  if(ret != 0) exit(1);
}




using namespace ana;

struct Bounds_t {
  double min;
  double max = 0;
};

void DrawArrow(TPad * p, double x, TString arrow_opt, int line_color = kGray+1, int line_style = 1, int line_width = 3);

void DrawBounds(TPad * p, double x, int line_color = kGray+1, int line_style = 1, int line_width = 3);

void datamc(std::string input_file_name = "",
            std::string plot_dump = "./plots/")
{
  load_libs_muonid();
  std::map<std::string, const HistAxis*> axes;
  std::map<std::string, Bounds_t> bounds;

  bounds["vtxx"] = {nuebarccinc::vNueCCIncFiducialMin.X(),
                    nuebarccinc::vNueCCIncFiducialMax.X()};
  bounds["vtxy"] = {nuebarccinc::vNueCCIncFiducialMin.Y(),
                    nuebarccinc::vNueCCIncFiducialMax.Y()};
  bounds["vtxz"] = {nuebarccinc::vNueCCIncFiducialMin.Z(),
                    nuebarccinc::vNueCCIncFiducialMax.Z()};
  bounds["top"   ] = {50 };
  bounds["bottom"] = {30 };
  bounds["east"  ] = {50 };
  bounds["west"  ] = {30 };
  bounds["front" ] = {150};
  
  axes["vtxx"] = new HistAxis("Vertex in X Direction (cm)",
                              Binning::Simple(40, -200, 200),
                              kVtxX);
  axes["vtxy"] = new HistAxis("Vertex in Y Direction (cm)",
                              Binning::Simple(40, -200, 200),
                              kVtxY);

  axes["vtxz"] = new HistAxis("Vertex in Z Direction (cm)",
                              Binning::Simple(80, 0, 1200),
                              kVtxZ);
  axes["top"] = new HistAxis("Min. Prong Distance from Top (cm)",
                             Binning::Simple(40, 0, 400),
                             kDistAllTop);
  axes["bottom"] = new HistAxis("Min. Prong Distance from Bottom (cm)",
                                Binning::Simple(40, 0, 400),
                                kDistAllBottom);
  axes["east"] = new HistAxis("Min. Prong Distance from East (cm)",
                              Binning::Simple(40, 0, 400),
                              kDistAllEast);
  axes["west" ] = new HistAxis("Min. Prong Distance from West (cm)",
                               Binning::Simple(40, 0, 400),
                               kDistAllWest);
  axes["front"] = new HistAxis("Min. Prong Distance from Front (cm)",
                               Binning::Simple(80, 0, 1200),
                               kDistAllFront);

  std::map<std::string, const Cut *> cuts;
  cuts["selection"] = new Cut(nuebarccinc::kDQ                          &&
                              nuebarccinc::kFrontPlanes                 &&
                              nuebarccinc::kNHits                       &&
                              nuebarccinc::kMuonIDProd4Cut              &&
                              nuebarccinc::kNuebarCCIncContainmentLoose && 
                              nuebarccinc::kNuebarCCIncFiducialLoose    );
  cuts["selection_and_cvn"] = new Cut(*cuts.at("selection") && nuebarccinc::kCVNNueIDCut);
  if(input_file_name == "") {
    std::map<std::string, SpectrumLoader*> loaders;
    loaders["mc"]   = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_NOMINAL);
    loaders["data"] = new SpectrumLoader(nuebarccinc::PROD5_RHC_DATA);

    // central value weight to be applied to all samples
    const Var * kCVWeight = new Var(kPPFXFluxCVWgt * kXSecCVWgt2020);

    std::map<std::string, Spectrum*> specs;
    std::string key = "";
    for(auto axis_it = axes.begin(); axis_it != axes.end(); axis_it++) {
      for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
        key = "data_" + cut_it->first + "_" + axis_it->first;
        specs[key]  = new Spectrum(*loaders["data"],
                                   *axis_it->second,
                                   *cut_it->second,
                                   kNoShift,
                                   *kCVWeight);
        for(auto ichn = 0; ichn < nuebarccinc::kNumChns; ichn++) {
          key = "mc_" + cut_it->first + "_" + axis_it->first + "_" + nuebarccinc::chns[ichn].name;
          specs[key]  = new Spectrum(*loaders["mc"],
                                     *axis_it->second,
                                     *cut_it->second && nuebarccinc::chns[ichn].cut,
                                     kNoShift,
                                     *kCVWeight);
      
        }
      }
    }

    loaders.at("data")->Go();
    loaders.at("mc"  )->Go();

    TFile * out = new TFile("datamc.root", "recreate");
    for(auto spec_it = specs.begin(); spec_it != specs.end(); spec_it++) {
      spec_it->second->SaveTo(out, spec_it->first.c_str());
    }
  
    out->Close();
  }

  else {

    TFile * input = TFile::Open(input_file_name.c_str());

    std::map<std::string, TH1D*> hists_data;
    std::map<std::string, TH1D*> hists_mc;
    std::string key;
    for(auto axis_it = axes.begin(); axis_it != axes.end(); axis_it++) {
      for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
        key = "data_" + cut_it->first + "_" + axis_it->first;
        auto tmp = Spectrum::LoadFrom(input, key.c_str());
        double POT = tmp->POT();
        hists_data[key] = tmp->ToTH1(POT);
        for(auto ichn = 0; ichn < nuebarccinc::kNumChns; ichn++) {
          key = "mc_" + cut_it->first + "_" + axis_it->first + "_" + nuebarccinc::chns[ichn].name;
          tmp = Spectrum::LoadFrom(input, key.c_str());
          hists_mc[key] = tmp->ToTH1(POT);
        }
      }
    }

    // rebin vtx z and min distance from front
    for(std::string axis : {"front", "vtxz"}) {
      for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
        key = "data_" + cut_it->first + "_" + axis;
        hists_data.at(key)->Rebin(4);
        for(auto ichn = 0; ichn < nuebarccinc::kNumChns; ichn++) {
          key = "mc_" + cut_it->first + "_" + axis + "_" + nuebarccinc::chns[ichn].name;
          hists_mc.at(key)->Rebin(4);
        }
      }
    }
    
    for(auto axis_it  = axes.begin(); axis_it != axes.end(); axis_it++) {
      for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
        TCanvas * c = new TCanvas(UniqueName().c_str());
        TPad * p1 = new TPad("p1", "", 0, 0, 1, 1);
        TPad * p2 = new TPad("p2", "", 0, 0, 1, 1);
        p1->SetLeftMargin(0.15);
        p2->SetLeftMargin(0.15);
        p1->SetBottomMargin(0.35);
        p2->SetTopMargin(0.65);
        p1->SetFillStyle(0);
        p2->SetFillStyle(0);
        TLegend * leg = new TLegend(.65, .89, .89, .65, "", "NDC");
        leg->SetNColumns(2);
        leg->SetFillStyle(0);

        p1->cd();
        key = "data_" + cut_it->first + "_" + axis_it->first;
        auto data = hists_data.at(key);
        
        data->SetMarkerStyle(8);
        data->SetMarkerSize(1);
        data->GetYaxis()->CenterTitle();
        data->GetXaxis()->CenterTitle();
        data->GetYaxis()->SetTitleOffset(1.15);
        auto ratio = (TH1D*) data->Clone(UniqueName().c_str());

        data->GetYaxis()->SetRangeUser(0.0001, data->GetMaximum() * 1.5);
        data->Draw("e0 p");
        //      data->SetTickLength(0);
        data->SetLabelSize(0);

        auto dummy = (TH1D*) data->Clone(UniqueName().c_str());
        THStack * stack = new THStack((key + "_stack").c_str(), "");

        // add nonfiducial to nuebar
        hists_mc.at("mc_" + cut_it->first + "_" + axis_it->first + "_nuebarcc")->Add(hists_mc.at("mc_" + cut_it->first + "_" + axis_it->first + "_nuebarcc_nonfiducial"));

        for(auto ichn = 1; ichn < nuebarccinc::kNumChns; ichn++) {
          if(ichn == 2) continue;       // don't fill MC total and skip nonfiducial
          key = "mc_" + cut_it->first + "_" + axis_it->first + "_" + nuebarccinc::chns[ichn].name;
          hists_mc.at(key)->SetFillColor(nuebarccinc::chns_style[ichn].color);
          auto entry = leg->AddEntry((TObject*)0, nuebarccinc::chns_style[ichn].latex.c_str(), "");
          entry->SetTextColor(nuebarccinc::chns_style[ichn].color);
          stack->Add(hists_mc.at(key));
        }
        leg->AddEntry((TObject*)0, "", "");
        leg->AddEntry(data, "Data", "lp");
        stack->Draw("hist same");
        data->Draw("ep same");

        leg->Draw();

        p2->cd();
        ratio->Divide(hists_mc.at("mc_" + cut_it->first + "_" + axis_it->first + "_mc"));
        ratio->GetYaxis()->SetRangeUser(.5, 1.5);
        ratio->GetYaxis()->SetTitle("Data / MC");
        ratio->Draw("ep");
        nuebarccinc::One();

        c->cd();
        p1->Draw();
        p2->Draw("same");

        DrawArrow(p1, bounds.at(axis_it->first).min, ">");
        if(bounds.at(axis_it->first).max)
          DrawArrow(p1, bounds.at(axis_it->first).max, "<");

        DrawBounds(p2, bounds.at(axis_it->first).min);
        if(bounds.at(axis_it->first).max)
          DrawBounds(p2, bounds.at(axis_it->first).max);

        c->SaveAs((plot_dump + "/datamc_" + cut_it->first + "_" + axis_it->first + ".pdf").c_str());
      }
    }


  }
}

void nuebarccinc::One()
{
  TGraph* one = new TGraph();
  one->SetPoint(0, -1e6, 1);
  one->SetPoint(1,  1e6, 1);
  one->SetLineStyle(7);
  one->SetLineColor(kGray+3);
  one->Draw("same");
}

void DrawBounds(TPad * p, double x, int line_color, int line_style, int line_width)
{
  p->cd();
  p->Update();

  TGraph* g = new TGraph();
  g->SetPoint(0, x, p->GetY1());
  g->SetPoint(1, x, p->GetY2());
  g->SetLineStyle(line_style);
  g->SetLineColor(line_color);
  g->SetLineWidth(line_width);
  g->Draw("same");

  g->Draw();
  p->Update();
}

void DrawArrow(TPad * p, double x, TString arrow_opt, int line_color, int line_style, int line_width)
{
  p->cd();
  p->Update();
  bool left = arrow_opt.Contains("<");
  double offset = 0.07 * (p->GetX2() - p->GetX1());
  if(left) offset *= -1;

  double y = p->GetY1() + 0.8 * (p->GetY2() - p->GetY1());

  std::cout << p->GetY1() << "\t" << p->GetY2() << std::endl;

  TLine * l1 = new TLine(x, 0, x, y);
  TArrow * arrow = new TArrow(x, y, x+offset, y, 0.04, arrow_opt.Data());
  if(left) arrow = new TArrow(x+offset, y, x, y, 0.04, arrow_opt.Data());  
  
  arrow->SetLineColor(line_color);
  l1   ->SetLineColor(line_color);

  arrow->SetLineStyle(line_style);
  l1   ->SetLineStyle(line_style);

  arrow->SetLineWidth(line_width);
  l1   ->SetLineWidth(line_width);

  arrow->Draw();
  l1->Draw();
  p->Update();
}