drawLongTerm_2019.C

Go to the documentation of this file.

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Cuts/NumuCuts2018.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Vars/NumuVars.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Plots.h"

#include "TStyle.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TProfile.h"
#include "TSpectrum.h"

#include "VarsAndCuts_2019.h"
#include "Periods_2019.h"
#include "CustomFunctions_2019.h"

using namespace ana;

void drawLongTerm_2019(std::string samDef, unsigned int endTime, unsigned int nDays, std::string outDir)
{
  //SET SIG FIGS.
  gStyle->SetStatFormat("6.3g");
  gStyle->SetPaintTextFormat("6.3g");
  gStyle->SetFitFormat("6.3g");

  //DEFINE A TFILE TO WRITE ALL CANVASES TO.
  TFile *f_Out = new TFile((TString)(outDir+"BeamMonitoring_LongTerm.root"),"RECREATE");
  TDirectory *dir_Profile = f_Out->mkdir("Profiles");
  TDirectory *dir_1D      = f_Out->mkdir("1Ds");
  TDirectory *dir_Ratio   = f_Out->mkdir("Ratios");

  //APPLY TIME CUTS TO DATASET, DEFINE THE LOADER AND SET SPILL CUTS.
  unsigned int startTime = endTime - 86400*nDays;
  const std::string fFiles = "dataset_def_name_newest_snapshot " + samDef + Form(" and online.subrunendtime >= '%u' and online.subrunendtime <= '%u'", startTime, endTime);
  SpectrumLoader loader(fFiles);
  loader.SetSpillCut(kStandardSpillCuts);

  const Binning kTimeBinning = Binning::Simple(nDays, startTime, endTime);
  const Cut     kNumuNDCut = kNumuCutND2018;

  //PLOT SPILL TIME FOR A PARTICULAR PERIOD SO THAT WE HAVE A POT TO NORMALISE BY.
  //CUT USING getYearAndMonth.
  std::vector<TH1D*> vec_hPOT;

  for(int i = 0; i < kNumPeriods; i++)
  {
    TH1D *hPOT = new TH1D(Form("hPOT_Period%i", i),";Central Time;POT (10^{15})", nDays, startTime, endTime);
    vec_hPOT.push_back(hPOT);

    loader.AddSpillHistogram(vec_hPOT.at(i), kSpillTime, getYearAndMonth(datasets[i].fYear, datasets[i].fMonth), kPOTe15);
  }

  //DEFINE PLOTS OF THE VARIABLES WE WANT TO SHOW.
  std::vector<Plot> vec_Plots = {
    {"Central Time", "Reconstructed Neutrino Energy, (GeV)", "ccE",      Binning::Simple(10,0,5),   kCCE,       kNumuNDCut,  1.5,  2.2 },
    {"Central Time", "Reconstructed Hadronic Energy, (GeV)", "hadE",     Binning::Simple(12,0,3),   kHadE,      kNumuNDCut,  0.2,  0.6 },
    {"Central Time", "Reconstructed Muon Energy, (GeV)",     "muonE",    Binning::Simple(15,0,4.5), kMuE,       kNumuNDCut,  0.9,  1.7 },
    {"Central Time", "Hadronic Energy Fraction",             "hadEfrac", Binning::Simple(10,0,1.0), kHadEFrac,  kNumuNDCut,  0.05, 0.35},
    {"Central Time", "No. of Hits",                          "nhit",     Binning::Simple(10,0,250), kNHit,      kNumuNDCut,  85,   105 },
    {"Central Time", "Vertex X, (m)",                        "startX",   Binning::Simple(10,-2,2),  kTrkStartX, kNumuNDCut, -0.40, 0.0 },
    {"Central Time", "Vertex Y, (m)",                        "startY",   Binning::Simple(10,-2,2),  kTrkStartY, kNumuNDCut, -0.25, 0.15},
    {"Central Time", "Vertex Z, (m)",                        "startZ",   Binning::Simple(12,0,12),  kTrkStartZ, kNumuNDCut,  4.5,  5.8 }
  };

  std::map<std::pair<unsigned int,unsigned int>,Spectrum*> map_PeriodAndVarToSpectrum;
  std::map<std::string,TCanvas*> map_VarToCanvas;
  std::map<std::string,TH1D*>    map_VarToMeanHist;

  for(unsigned int i = 0; i < vec_hPOT.size(); i++)
  {
    for(unsigned int j = 0; j < vec_Plots.size(); j++)
    {
      Spectrum *h_Spectrum = new Spectrum(vec_Plots.at(j).fLabelX, vec_Plots.at(j).fLabelY, loader, kTimeBinning, kTime, vec_Plots.at(j).fBins, vec_Plots.at(j).fVar, 
                                         kYear == datasets[i].fYear && kMonth == datasets[i].fMonth && vec_Plots.at(j).fCut);
      map_PeriodAndVarToSpectrum[{i,j}] = h_Spectrum;
      if(i==0)
      {
        map_VarToCanvas[vec_Plots.at(j).fName] = new TCanvas((TString)("c_"+vec_Plots.at(j).fName+"_Profile"),(TString)("c_"+vec_Plots.at(j).fName+"_Profile"));
        map_VarToCanvas[vec_Plots.at(j).fName]->SetBottomMargin(0.15);
        map_VarToCanvas[vec_Plots.at(j).fName]->SetLeftMargin  (0.15);
        map_VarToCanvas[vec_Plots.at(j).fName]->SetRightMargin (0.15);
        map_VarToCanvas[vec_Plots.at(j).fName+"_1D"] = new TCanvas((TString)("c_"+vec_Plots.at(j).fName+"_1D"),(TString)("c_"+vec_Plots.at(j).fName+"_1D"),2*map_VarToCanvas[vec_Plots.at(j).fName]->GetWw(), 
                                                                                                                                                             map_VarToCanvas[vec_Plots.at(j).fName]->GetWh());
        map_VarToCanvas[vec_Plots.at(j).fName+"_1D"]->SetBottomMargin(0.15);
        map_VarToCanvas[vec_Plots.at(j).fName+"_1D"]->SetLeftMargin  (0.15);
        map_VarToCanvas[vec_Plots.at(j).fName+"_1D"]->SetRightMargin (0.15);
        map_VarToCanvas[vec_Plots.at(j).fName+"_Ratio"] = new TCanvas((TString)("c_"+vec_Plots.at(j).fName+"_Ratio"),(TString)("c_"+vec_Plots.at(j).fName+"_Ratio"),2*map_VarToCanvas[vec_Plots.at(j).fName]->GetWw(), 
                                                                                                                                                             map_VarToCanvas[vec_Plots.at(j).fName]->GetWh());
        map_VarToCanvas[vec_Plots.at(j).fName+"_Ratio"]->SetBottomMargin(0.15);
        map_VarToCanvas[vec_Plots.at(j).fName+"_Ratio"]->SetLeftMargin  (0.15);
        map_VarToCanvas[vec_Plots.at(j).fName+"_Ratio"]->SetRightMargin (0.15);

        TH1D *h_Mean = new TH1D();
        h_Mean->SetName((TString)vec_Plots.at(j).fName);
        h_Mean->SetLineWidth(1);
        h_Mean->SetLineColor(kBlack);
        map_VarToMeanHist[vec_Plots.at(j).fName] = h_Mean;
      }
    }
  } 

  loader.Go();

  //PROFILES.
  TLegend* leg = new TLegend(0.17,0.17,0.83,0.4);
  leg->SetFillColor(0);
  leg->SetNColumns(3);
  TLegend* leg_1D_1 = new TLegend(0.6,0.6,0.85,0.85);
  leg_1D_1->SetFillColor(0);
  leg_1D_1->SetNColumns(3);
  TLegend* leg_1D_2 = new TLegend(0.35,0.20,0.62,0.48);
  leg_1D_2->SetFillColor(0);
  leg_1D_2->SetNColumns(3);

  for(auto spectrum : map_PeriodAndVarToSpectrum)
  {
    TH2      *h_2D   = spectrum.second->ToTH2(spectrum.second->POT());
    TProfile *h_Prof = h_2D->ProfileX((TString)("prof_"+vec_Plots.at(spectrum.first.second).fName+"_"+Form("%i_%i", datasets[spectrum.first.first].fYear, datasets[spectrum.first.first].fMonth)));
    h_Prof->SetMarkerStyle(datasets[spectrum.first.first].fStyle);
    h_Prof->SetMarkerColor(datasets[spectrum.first.first].fColor);
    h_Prof->SetLineColor  (datasets[spectrum.first.first].fColor);
    h_Prof->SetMarkerSize (1.4);
    if(spectrum.first.first==0)
    {
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
      h_Prof->SetTitle("");
      h_Prof->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(0));
      h_Prof->GetYaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
      h_Prof->GetXaxis()->CenterTitle();
      h_Prof->GetYaxis()->CenterTitle();
      h_Prof->GetXaxis()->SetTimeDisplay(1);
      h_Prof->GetXaxis()->SetTimeFormat("%y-%m-%d%F1970-01-01 00:00:00");
      h_Prof->Draw();

      unsigned int minBin  = h_Prof->GetXaxis()->FindBin(startTime);
      unsigned int maxBin  = h_Prof->GetXaxis()->FindBin(endTime);
      double       minZoom = h_Prof->GetXaxis()->GetBinLowEdge(minBin-1);
      double       maxZoom = h_Prof->GetXaxis()->GetBinLowEdge(maxBin+2);
      h_Prof->GetXaxis()->SetRangeUser(minZoom, maxZoom);
      h_Prof->GetYaxis()->SetRangeUser(vec_Plots.at(spectrum.first.second).fYRangeLow, vec_Plots.at(spectrum.first.second).fYRangeHigh);
    }
    else
    {
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
      h_Prof->Draw("SAME");
    }

    if(spectrum.first.second==0)
    {
      leg->AddEntry(h_Prof, (TString)datasets[spectrum.first.first].fLabel, "P");
      leg_1D_1->AddEntry(h_Prof, (TString)datasets[spectrum.first.first].fLabel, "P");
      leg_1D_2->AddEntry(h_Prof, (TString)datasets[spectrum.first.first].fLabel, "P");
    }
    if(spectrum.first.first==datasets.size()-1)
    {
      leg->Draw();
      lastUpdated(endTime);
    }
  }
  TH1D *h_Dummy = new TH1D("h_Dummy","h_Dummy",1,0,1);
  h_Dummy->SetLineColor(kBlack);
  leg_1D_1->AddEntry(h_Dummy, "Average (Mean)", "L");
  leg_1D_2->AddEntry(h_Dummy, "Average (Mean)", "L");

  //1D.
  for(auto spectrum : map_PeriodAndVarToSpectrum)
  {
    TH2  *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
    TH1D *h_1D = h_2D->ProjectionY((TString)("1D_"+vec_Plots.at(spectrum.first.second).fName+"_"+Form("%i_%i", datasets[spectrum.first.first].fYear, datasets[spectrum.first.first].fMonth)));
    if(spectrum.first.first==0)
    {
      map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName] = (TH1D*)h_1D->Clone((TString)("h_Mean"+vec_Plots.at(spectrum.first.second).fName));
    }
    else
    {
      map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Add(h_1D);
    }

    if(vec_hPOT.at(spectrum.first.first)->Integral()!=0.){h_1D->Scale(1.0/vec_hPOT.at(spectrum.first.first)->Integral());}
    TGraphErrors *gr       = new TGraphErrors(h_1D);                 
    gr->SetMarkerStyle(datasets[spectrum.first.first].fStyle); 
    gr->SetMarkerColor(datasets[spectrum.first.first].fColor); 
    gr->SetLineColor  (datasets[spectrum.first.first].fColor); 

    for(int j = 0; j < gr->GetN(); j++)
    {
      double x, y, dy;
      double shift = (h_1D->GetBinWidth(1)/(double)datasets.size())*((spectrum.first.first+1)*1. - datasets.size()*0.5);
      gr->GetPoint(j, x, y);
      dy = gr->GetErrorY(j);
      gr->SetPoint(j, x + shift, y);
      gr->SetPointError(j, 0, dy);
    }
    if(spectrum.first.first==0)
    {
      gr->SetTitle("");
      gr->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
      gr->GetYaxis()->SetTitle("Slices / 10^{15} POT");
      gr->GetXaxis()->CenterTitle();
      gr->GetYaxis()->CenterTitle();
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_1D"]->cd();
      gr->Draw("AP");
    }
    else
    {
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_1D"]->cd();
      gr->Draw("P");
    }

    if(spectrum.first.first==datasets.size()-1)
    {
      map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Scale(1.0e15/map_PeriodAndVarToSpectrum[{0, spectrum.first.second}]->POT());
      map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->SetTitle("");
      map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Draw("HIST SAME");
      if(vec_Plots.at(spectrum.first.second).fName=="ccE"      ||
         vec_Plots.at(spectrum.first.second).fName=="muonE"    ||
         vec_Plots.at(spectrum.first.second).fName=="hadE"     ||
         vec_Plots.at(spectrum.first.second).fName=="hadEfrac" ||
         vec_Plots.at(spectrum.first.second).fName=="nhit") {
        leg_1D_1->Draw();
      }
      else{
        leg_1D_2->Draw();
      }                 

      lastUpdated(endTime);
    }
  }
  //RATIOS.
  for(auto spectrum : map_PeriodAndVarToSpectrum)
  {
    TH2  *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
    TH1D *h_1D = h_2D->ProjectionY((TString)("1D_"+vec_Plots.at(spectrum.first.second).fName+"_"+Form("%i_%i", datasets[spectrum.first.first].fYear, datasets[spectrum.first.first].fMonth)));
    if(vec_hPOT.at(spectrum.first.first)->Integral()!=0.){h_1D->Scale(1.0/vec_hPOT.at(spectrum.first.first)->Integral());}
    TGraphErrors *gr = new TGraphErrors(h_1D);                 
    gr->SetMarkerStyle(datasets[spectrum.first.first].fStyle); 
    gr->SetMarkerColor(datasets[spectrum.first.first].fColor); 
    gr->SetLineColor  (datasets[spectrum.first.first].fColor); 

    for(int j = 0; j < gr->GetN(); j++)
    {
      double x, y, my, dy;
      double shift = (h_1D->GetBinWidth(1)/(double)datasets.size())*((spectrum.first.first+1)*1. - datasets.size()*0.5);
      my = map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->GetBinContent(j+1);
      gr->GetPoint(j, x, y);
      if(my>0)
      {
        y /= my;
        dy = gr->GetErrorY(j)/my;
        gr->SetPoint(j, x + shift, y);
        gr->SetPointError(j, 0, dy);
      }
    }
    double xMin, xMax;
    if(spectrum.first.first==0)
    {
      map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName] = (TH1D*)h_1D->Clone((TString)("h_Mean"+vec_Plots.at(spectrum.first.second).fName));
      gr->SetTitle("");
      gr->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
      gr->GetYaxis()->SetTitle("Slices per 10^{15} POT / Average");
      gr->GetXaxis()->CenterTitle();
      gr->GetYaxis()->CenterTitle();
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_Ratio"]->cd();
      gr->Draw("AP");
      gr->GetYaxis()->SetRangeUser(0.8,1.2);

      double binsArrary[h_1D->GetNbinsX()];
      h_1D->GetLowEdge(binsArrary);
      for(int j = 0; j < h_1D->GetNbinsX(); j++)
      {
        TLine *l = new TLine(binsArrary[j], 0.8, binsArrary[j], 1.2);
        l->SetLineStyle(kDashed);
        l->Draw();
      }
    }
    else
    {
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_Ratio"]->cd();
      gr->Draw("P");
    }

    if(spectrum.first.first==datasets.size()-1)
    {
      lastUpdated(endTime);
      TGraph *xAxis = new TGraph;
      xAxis->SetPoint(0, -1.e99, 1);
      xAxis->SetPoint(1,  1.e99, 1);
      xAxis->Draw("L");
    }
  }

  for(auto it : map_VarToCanvas)
  {
    it.second->SaveAs((TString)(outDir+it.first+"_LongTerm.png"));
    if     (((std::string)(it.second->GetName())).find("Profile")!=std::string::npos){dir_Profile->cd();}
    else if(((std::string)(it.second->GetName())).find("1D")     !=std::string::npos){dir_1D     ->cd();}
    else if(((std::string)(it.second->GetName())).find("Ratio")  !=std::string::npos){dir_Ratio  ->cd();}
    it.second->Write();
    it.second->Clear();
  }

  std::cout <<  "ALL PLOTS MADE" << std::endl;

}