drawLongTerm.C

Go to the documentation of this file.

#ifdef __CINT__
void drawLongTerm(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming){
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#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.h"
#include "Periods.h"
#include "CustomFunctions.h"

using namespace ana;

void drawLongTerm(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming)
{

  std::cout << "Drawing plots..." << std::endl;

  gStyle->SetStatFormat("6.3g"); // Significant figures. Fix if this is not satisfactory
  gStyle->SetPaintTextFormat("6.3g");
  gStyle->SetFitFormat("6.3g");

  // Output directory
  const std::string outDir = "/nova/ana/users/novadq/BeamMonitoring/Plots/stability/";
  // Loader for data
  const std::string fFiles = "dataset_def_name_newest_snapshot " + samdef + Form(" and online.subrunendtime >= '%u' and online.subrunendtime <= '%u'", timestamp - 86400*nDays, timestamp);

  // int nBins = nDays/7 + 1; // weekly looks better than daily
  int nBins = nDays;

  const Binning kTimeBinning = Binning::Simple(nBins, timestamp - 86400*nDays, timestamp);

  const Cut kNumuNDCut = kNumuCutND2018;

  SpectrumLoader loader(fFiles);
  loader.SetSpillCut(kStandardSpillCuts);

  TH1D* hPOT[kNumPeriods];

  for(int j = 0; j < kNumPeriods; ++j)
  {
    hPOT[j] = new TH1D(Form("pot_%i",j+1),";Central Time;POT (10^{15})", nBins, timestamp - 86400*nDays, timestamp);

    loader.AddSpillHistogram(hPOT[j], kSpillTime, getYearAndMonth(datasets[j].year, datasets[j].month), kPOTe15);
  }

  const int kNumPlots = 8; 

  Plot plots[kNumPlots] = {
    {";Central Time;Reconstructed neutrino energy (GeV)","ccE", Binning::Simple(10,0,5), kCCE, kNumuNDCut, 1.5, 2.0},
    {";Central Time;Reconstructed hadronic energy (GeV)","hadE", Binning::Simple(12,0,3), kHadE, kNumuNDCut, 0.15, 0.50},
    {";Central Time;Reconstructed muon energy (GeV)","muonE", Binning::Simple(15,0,4.5), kMuE, kNumuNDCut, 1.0, 2.0},
    {";Central Time;Hadronic energy fraction","hadEfrac", Binning::Simple(10,0,1.0), kHadEFrac, kNumuNDCut, 0, 0.25},
    {";Central Time;No. of hits","nhit", Binning::Simple(10,0,250), kNHit, kNumuNDCut, 90, 110},
    {";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, 5.0, 6.25}
  };

  /// ====================================================================== ///

  Spectrum* hData[kNumPlots][kNumPeriods];

  for(int i = 0; i < kNumPlots; ++i){
    for(int j = 0; j < kNumPeriods; ++j){
      hData[i][j] = new Spectrum(plots[i].label, loader, kTimeBinning, kTime, plots[i].bins, plots[i].var, kYear == datasets[j].year && kMonth == datasets[j].month && plots[i].cut);
    } // end loop over j
  } // end loop over i

  // GO!
  loader.Go();

for(int i = 0; i < kNumPlots; ++i){ 

  std::string fullName = plots[i].fname+"_"+naming;
  TCanvas* c = new TCanvas("c",(fullName).c_str());
  c->SetBottomMargin(0.15);
  c->SetLeftMargin(0.15);
  c->SetRightMargin(0.15);

  TLegend* leg = new TLegend(0.17,0.17,0.83,0.4);
  leg->SetFillColor(0);
  leg->SetNColumns(3);

  // First the profile
  for(int j = 0; j < kNumPeriods; ++j){

    std::string monthName = plots[i].fname+Form("_period%i_",j+1)+naming;
    
    TH2* Histogram = hData[i][j]->ToTH2(hData[i][j]->POT());

    TProfile* Profile = Histogram->ProfileX((monthName+"Profile").c_str());
    Profile->SetMarkerStyle(datasets[j].style);
    Profile->SetMarkerColor(datasets[j].color);
    Profile->SetLineColor(datasets[j].color);
    Profile->SetMarkerSize(1.4);

    Profile->GetXaxis()->CenterTitle();
    Profile->GetYaxis()->CenterTitle();
    Profile->GetYaxis()->SetTitle(Histogram->GetYaxis()->GetTitle());
    Profile->GetXaxis()->SetTimeDisplay(1);
    Profile->GetXaxis()->SetTimeFormat("%y-%m-%d%F1970-01-01 00:00:00");

    unsigned int minBin = Profile->GetXaxis()->FindBin(timestamp - 86400*nDays);
    unsigned int maxBin = Profile->GetXaxis()->FindBin(timestamp);

    double minZoom = Profile->GetXaxis()->GetBinLowEdge(minBin-1);
    double maxZoom = Profile->GetXaxis()->GetBinLowEdge(maxBin+2);

    Profile->GetXaxis()->SetRangeUser(minZoom, maxZoom);

    if(j == 0) 
      {
        Profile->GetYaxis()->SetRangeUser(plots[i].yRangeLow, plots[i].yRangeHigh);
        Profile->Draw();
      }

    else Profile->Draw("same");

    leg->AddEntry(Profile, datasets[j].label.c_str(), "p");

    } // End loop over j

    leg->Draw();
    lastUpdated(timestamp);

    //c->SaveAs(("stability/"+fullName+".png").c_str());
    //c->SaveAs(("stability/"+fullName+".pdf").c_str());
    //c->SaveAs(("stability/"+fullName+".root").c_str());
    c->SaveAs((outDir+fullName+".png").c_str());
    c->SaveAs((outDir+fullName+".pdf").c_str());
    c->SaveAs((outDir+fullName+".root").c_str());

    // And finally the 1D plots
    double ww = c->GetWw();
    double wh = c->GetWh();
    c->Close();
    TCanvas* c2 = new TCanvas("c2",(fullName).c_str(), 2*ww, wh);
    c2->SetBottomMargin(0.15);
    c2->SetLeftMargin(0.15);
    c2->SetRightMargin(0.15);

    c2->SetWindowSize(2*ww, wh); // wide for spectra

    TLegend* leg2;
    if(i>4) leg2 = new TLegend(0.35,0.20,0.62,0.48);
    else leg2 = new TLegend(0.55,0.60,0.82,0.88);
    leg2->SetFillColor(0);
    leg2->SetNColumns(3);

    std::string dim1Name = plots[i].fname+"_1D_"+naming;

    TH1* hMean1D = new TH1();

    for(int j = 0; j < kNumPeriods; ++j){
      std::string dim1NameMonth = plots[i].fname+Form("_period%i_1D_",j+1)+naming;
      
      TH2* Histogram = hData[i][j]->ToTH2(hData[i][j]->POT());

      TH1* Hist1D = Histogram->ProjectionY((dim1NameMonth).c_str());

      if(j == 0) hMean1D = (TH1*) Hist1D->Clone("Mean_1D");
      else hMean1D->Add(Hist1D);

      Hist1D->Scale(1.0/hPOT[j]->Integral());

      TGraphErrors* Graph = new TGraphErrors(Hist1D);
      Graph->SetName((dim1NameMonth+"_gr").c_str());
      Graph->SetMarkerStyle(datasets[j].style);
      Graph->SetMarkerColor(datasets[j].color);
      Graph->SetLineColor(datasets[j].color);
      Graph->GetXaxis()->SetTitle(Hist1D->GetXaxis()->GetTitle());
      Graph->GetYaxis()->SetTitle("Slices / 10^{15} POT");
      Graph->GetXaxis()->CenterTitle();
      Graph->GetYaxis()->CenterTitle();

      for(int fPoint = 0; fPoint<Graph->GetN(); ++fPoint)
      {
        double x, y, dy;
        double shift = (Hist1D->GetBinWidth(1) / kNumPeriods) * ((j+1)*1. - kNumPeriods*0.5);
        Graph->GetPoint(fPoint, x, y);
        dy = Graph->GetErrorY(fPoint);
        Graph->SetPoint(fPoint, x + shift, y);
        Graph->SetPointError(fPoint,0,dy);
      }

      if(j == 0)
      {
        Hist1D->SetLineWidth(0);
        Hist1D->SetLineColor(10);
        Hist1D->SetMarkerStyle(0);
        Hist1D->SetMarkerColor(0);
        Hist1D->Draw("hist");
      }
      
      Graph->Draw("P");

      leg2->AddEntry(Graph, datasets[j].label.c_str(), "p");

    } // End loop over j

    hMean1D->Scale(1.0e15/hData[i][0]->POT());
    hMean1D->SetLineWidth(1);
    hMean1D->SetLineColor(kBlack);
    hMean1D->Draw("hist same");

    leg2->AddEntry(hMean1D,"Average","l");
    leg2->Draw();
    lastUpdated(timestamp);

    //c2->SaveAs(("stability/"+dim1Name+".png").c_str());
    //c2->SaveAs(("stability/"+dim1Name+".pdf").c_str());
    //c2->SaveAs(("stability/"+dim1Name+".root").c_str());
    c2->SaveAs((outDir+dim1Name+".png").c_str());
    c2->SaveAs((outDir+dim1Name+".pdf").c_str());
    c2->SaveAs((outDir+dim1Name+".root").c_str());

    // Now the ratio
    c2->Clear();

    std::string ratioName = plots[i].fname+"_ratio_"+naming;

    for(int j = 0; j < kNumPeriods; ++j){
      std::string ratioNameMonth = plots[i].fname+Form("_period%i_ratio_",j+1)+naming;
      
      TH2* Histogram = hData[i][j]->ToTH2(hData[i][j]->POT());

      TH1* Hist1D = Histogram->ProjectionY((ratioNameMonth).c_str());

      Hist1D->Scale(1.0/hPOT[j]->Integral());

      TGraphErrors* Graph = new TGraphErrors(Hist1D);
      Hist1D->Reset();
      Graph->SetName((ratioNameMonth+"_gr").c_str());
      Graph->SetMarkerStyle(datasets[j].style);
      Graph->SetMarkerColor(datasets[j].color);
      Graph->SetLineColor(datasets[j].color);
      Graph->GetXaxis()->SetTitle(Hist1D->GetXaxis()->GetTitle());
      Graph->GetYaxis()->SetTitle("Slices per 10^{15} POT / Average");
      Graph->GetXaxis()->CenterTitle();
      Graph->GetYaxis()->CenterTitle();

      for(int fPoint = 0; fPoint<Graph->GetN(); ++fPoint)
      {
        double x, y, my, dy;
        double shift = (Hist1D->GetBinWidth(1) / kNumPeriods) * ((j+1)*1. - kNumPeriods*0.5);
        my = hMean1D->GetBinContent(fPoint+1);
        Graph->GetPoint(fPoint, x, y);
        if(my > 0)
        {
          y = y/my;
          dy = Graph->GetErrorY(fPoint)/my;
          Graph->SetPoint(fPoint, x + shift, y);
          Graph->SetPointError(fPoint,0,dy);
        }      
      }

      if(j == 0)
      {
        Hist1D->GetYaxis()->SetRangeUser(0.8,1.2);
        Hist1D->SetLineWidth(0);
        Hist1D->SetLineColor(10);
        Hist1D->SetMarkerStyle(0);
        Hist1D->SetMarkerColor(0);
        Hist1D->Draw("hist");
      }
      
      Graph->Draw("P");
    } // End loop over j

    TGraph* gr1 = new TGraph;
    gr1->SetPoint(0,-1.e99,1.);
    gr1->SetPoint(1, 1.e99,1.);
    gr1->SetLineColor(kBlack);
    gr1->SetLineWidth(1);
    gr1->SetFillColor(10);
    gr1->Draw("L");

    lastUpdated(timestamp);

    //c2->SaveAs(("stability/"+ratioName+".png").c_str());
    //c2->SaveAs(("stability/"+ratioName+".pdf").c_str());
    //c2->SaveAs(("stability/"+ratioName+".root").c_str());
    c2->SaveAs((outDir+ratioName+".png").c_str());
    c2->SaveAs((outDir+ratioName+".pdf").c_str());
    c2->SaveAs((outDir+ratioName+".root").c_str());

} // End loop over i

  std::cout << "All plots finished!" << std::endl;

} // End of function

#endif