drawIntensityEffect_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/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 "TTimeStamp.h"
#include "TPaveLabel.h"

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

using namespace ana;

void drawIntensityEffect_2019(std::string samdef, unsigned int endTime, unsigned int nDays, std::string outDir)
{
  gStyle->SetOptFit(1111);

  //DEFINE A FILE TO WRITE ALL OF THE CANVASES TO.
  TFile *f_Out = new TFile((TString)(outDir+"BeamMonitoring_IntensityEffect.root"),"RECREATE");

  //APPLY TIME CUTS TO DATAASET, 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 kSingleBin = Binning::Simple(1,0.5,1.5);
  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;
  std::vector<TH1D*> vec_hNumPOT;

  for(unsigned int i = 0; i < datasets.size(); i++)
  {
    TH1D *hPOT    = new TH1D(Form("hPOT_Period%i",    i),";POT (10^{15});POT (10^{15})", 1000,   0, 100);
    TH1D *hNumPOT = new TH1D(Form("hNumPOT_Period%i", i),";POT (10^{15});POT (10^{15})", 1,    0.5, 1.5);
    vec_hPOT   .push_back(hPOT);
    vec_hNumPOT.push_back(hNumPOT);

    loader.AddSpillHistogram(vec_hPOT.at(i),    kPOTe12,          getYearAndMonth(datasets[i].fYear, datasets[i].fMonth)         );
    loader.AddSpillHistogram(vec_hNumPOT.at(i), kSpillUnweighted, getYearAndMonth(datasets[i].fYear, datasets[i].fMonth), kPOTe12);
  }

  
  //DEFINE PLOTS OF THE VARIABLES WE WANT TO SHOW.
  std::vector<Plot> vec_Plots = {
    {"Central Time", "Number of slices / 10^{15} POT", "slices", kSingleBin, kUnweighted, kNumuNDCut, 0, 2.0}
  };

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

  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).fLabelX, loader, kSingleBin, kUnweighted, 
                                          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+"_IntensityEffect"),(TString)("c_"+vec_Plots.at(j).fName+"_IntensityEffect"));
        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+"_AllPeriods"] = new TCanvas((TString)("c_"+vec_Plots.at(j).fName+"_AllPeriods_IntensityEffect"),
                                                                           (TString)("c_"+vec_Plots.at(j).fName+"_AllPeriods_IntensityEffect"));
        map_VarToCanvas[vec_Plots.at(j).fName+"_AllPeriods"]->SetBottomMargin(0.15);
        map_VarToCanvas[vec_Plots.at(j).fName+"_AllPeriods"]->SetLeftMargin  (0.15);
        map_VarToCanvas[vec_Plots.at(j).fName+"_AllPeriods"]->SetRightMargin (0.15);
      }
    }
  } 

  loader.Go();

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

  TGraphErrors *g_AllPeriods = new TGraphErrors;
  g_AllPeriods->SetTitle("");
  g_AllPeriods->SetMarkerStyle(20);
  g_AllPeriods->SetMarkerColor(kBlack);
  g_AllPeriods->SetLineColor(kBlack);
  g_AllPeriods->GetXaxis()->CenterTitle();
  g_AllPeriods->GetYaxis()->CenterTitle();
  g_AllPeriods->GetXaxis()->SetTitle("Slices / 10^{9} POT");
  g_AllPeriods->GetYaxis()->SetTitle("Average POT / 10^{12}");

  f_Out->cd();
  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->Draw();
    }

    TGraphErrors *g = new TGraphErrors;
    g->SetName((TString)("1D_"+vec_Plots.at(spectrum.first.second).fName+"_"+Form("%i_%i", datasets[spectrum.first.first].fYear, datasets[spectrum.first.first].fMonth)));
    double y    = h_2D->Integral();
    double dy   = std::sqrt(y);
    double x    = vec_hPOT.at(spectrum.first.first)   ->GetMean();
    double npot = vec_hNumPOT.at(spectrum.first.first)->Integral()/1000.;
    g->SetTitle("");
    g->SetMarkerStyle(datasets[spectrum.first.first].fStyle);
    g->SetMarkerColor(datasets[spectrum.first.first].fColor);
    g->SetLineColor  (datasets[spectrum.first.first].fColor);
    if(x > 0 && npot > 0)
    {
      g->SetPoint     (0, x, y/npot);
      g->SetPointError(0, 0, dy/npot);

      g_AllPeriods->SetPoint     (g_AllPeriods->GetN(),   x, y/npot );
      g_AllPeriods->SetPointError(g_AllPeriods->GetN()-1, 0, dy/npot);
    }

    if(spectrum.first.first==0) 
    {
      TH2D* h_Canvas = new TH2D("h_Canvas","",100, 20, 70, 100, vec_Plots.at(spectrum.first.second).fYRangeLow, vec_Plots.at(spectrum.first.second).fYRangeHigh);
      h_Canvas->GetXaxis()->CenterTitle();
      h_Canvas->GetYaxis()->CenterTitle();
      h_Canvas->GetYaxis()->SetTitle("Slices / 10^{9} POT");
      h_Canvas->GetXaxis()->SetTitle("Average POT / 10^{12}");
      h_Canvas->SetLineColor(0);
      h_Canvas->SetLineWidth(0);

      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_AllPeriods"]->cd();
      h_Canvas->Draw();
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
      h_Canvas->Draw();

      g->Draw("P");
    }
    else
    {
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
      g->Draw("P");
    }
    if(spectrum.first.second==0)
    {
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
      leg->AddEntry(g, (TString)datasets[spectrum.first.first].fLabel, "P");
    }
    if(spectrum.first.first==datasets.size()-1)
    {
      leg->Draw();
      lastUpdated(endTime);
      map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_AllPeriods"]->cd();
      g_AllPeriods->Draw("P");
      g_AllPeriods->Fit("pol1");
      lastUpdated(endTime);
    }
  }

  for(auto it : map_VarToCanvas)
  {
    it.second->SaveAs((TString)(outDir+it.first+"_IntensityEffect.png"));
    it.second->Write();
    it.second->Clear();
  }

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