drawIntensityEffect_2019.C
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
4 #include "CAFAna/Cuts/NumuCuts2018.h"
5 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Vars/Vars.h"
10 #include "CAFAna/Analysis/Calcs.h"
11 #include "CAFAna/Analysis/Plots.h"
12 
13 #include "TStyle.h"
14 #include "TCanvas.h"
15 #include "TGraph.h"
16 #include "TGraphErrors.h"
17 #include "TGraphAsymmErrors.h"
18 #include "TFile.h"
19 #include "TH1.h"
20 #include "TH2.h"
21 #include "TLegend.h"
22 #include "TLatex.h"
23 #include "TProfile.h"
24 #include "TTimeStamp.h"
25 #include "TPaveLabel.h"
26 
27 #include "VarsAndCuts_2019.h"
28 #include "Periods_2019.h"
29 #include "CustomFunctions_2019.h"
30 
31 using namespace ana;
32 
33 void drawIntensityEffect_2019(std::string samdef, unsigned int endTime, unsigned int nDays, std::string outDir)
34 {
35  gStyle->SetOptFit(1111);
36 
37  //DEFINE A FILE TO WRITE ALL OF THE CANVASES TO.
38  TFile *f_Out = new TFile((TString)(outDir+"BeamMonitoring_IntensityEffect.root"),"RECREATE");
39 
40  //APPLY TIME CUTS TO DATAASET, DEFINE THE LOADER AND SET SPILL CUTS.
41  unsigned int startTime = endTime - 86400*nDays;
42  const std::string fFiles = "dataset_def_name_newest_snapshot " + samdef + Form(" and online.subrunendtime >= '%u' and online.subrunendtime <= '%u'", startTime, endTime);
43  SpectrumLoader loader(fFiles);
45 
46  const Binning kSingleBin = Binning::Simple(1,0.5,1.5);
47  const Cut kNumuNDCut = kNumuCutND2018;
48 
49 
50  //PLOT SPILL TIME FOR A PARTICULAR PERIOD SO THAT WE HAVE A POT TO NORMALISE BY.
51  //CUT USING getYearAndMonth.
52  std::vector<TH1D*> vec_hPOT;
53  std::vector<TH1D*> vec_hNumPOT;
54 
55  for(unsigned int i = 0; i < datasets.size(); i++)
56  {
57  TH1D *hPOT = new TH1D(Form("hPOT_Period%i", i),";POT (10^{15});POT (10^{15})", 1000, 0, 100);
58  TH1D *hNumPOT = new TH1D(Form("hNumPOT_Period%i", i),";POT (10^{15});POT (10^{15})", 1, 0.5, 1.5);
59  vec_hPOT .push_back(hPOT);
60  vec_hNumPOT.push_back(hNumPOT);
61 
62  loader.AddSpillHistogram(vec_hPOT.at(i), kPOTe12, getYearAndMonth(datasets[i].fYear, datasets[i].fMonth) );
63  loader.AddSpillHistogram(vec_hNumPOT.at(i), kSpillUnweighted, getYearAndMonth(datasets[i].fYear, datasets[i].fMonth), kPOTe12);
64  }
65 
66 
67  //DEFINE PLOTS OF THE VARIABLES WE WANT TO SHOW.
68  std::vector<Plot> vec_Plots = {
69  {"Central Time", "Number of slices / 10^{15} POT", "slices", kSingleBin, kUnweighted, kNumuNDCut, 0, 2.0}
70  };
71 
72  std::map<std::pair<unsigned int,unsigned int>,Spectrum*> map_PeriodAndVarToSpectrum;
73  std::map<std::string,TCanvas*> map_VarToCanvas;
74 
75  for(unsigned int i = 0; i < vec_hPOT.size(); i++)
76  {
77  for(unsigned int j = 0; j < vec_Plots.size(); j++)
78  {
79  Spectrum *h_Spectrum = new Spectrum(vec_Plots.at(j).fLabelX, vec_Plots.at(j).fLabelX, loader, kSingleBin, kUnweighted,
80  vec_Plots.at(j).fBins, vec_Plots.at(j).fVar, kYear == datasets[i].fYear && kMonth == datasets[i].fMonth && vec_Plots.at(j).fCut);
81  map_PeriodAndVarToSpectrum[{i,j}] = h_Spectrum;
82  if(i==0)
83  {
84  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"));
85  map_VarToCanvas[vec_Plots.at(j).fName]->SetBottomMargin(0.15);
86  map_VarToCanvas[vec_Plots.at(j).fName]->SetLeftMargin (0.15);
87  map_VarToCanvas[vec_Plots.at(j).fName]->SetRightMargin (0.15);
88  map_VarToCanvas[vec_Plots.at(j).fName+"_AllPeriods"] = new TCanvas((TString)("c_"+vec_Plots.at(j).fName+"_AllPeriods_IntensityEffect"),
89  (TString)("c_"+vec_Plots.at(j).fName+"_AllPeriods_IntensityEffect"));
90  map_VarToCanvas[vec_Plots.at(j).fName+"_AllPeriods"]->SetBottomMargin(0.15);
91  map_VarToCanvas[vec_Plots.at(j).fName+"_AllPeriods"]->SetLeftMargin (0.15);
92  map_VarToCanvas[vec_Plots.at(j).fName+"_AllPeriods"]->SetRightMargin (0.15);
93  }
94  }
95  }
96 
97  loader.Go();
98 
99  TLegend* leg = new TLegend(0.17,0.17,0.83,0.4);
100  leg->SetFillColor(0);
101  leg->SetNColumns(5);
102 
103  TGraphErrors *g_AllPeriods = new TGraphErrors;
104  g_AllPeriods->SetTitle("");
105  g_AllPeriods->SetMarkerStyle(20);
106  g_AllPeriods->SetMarkerColor(kBlack);
107  g_AllPeriods->SetLineColor(kBlack);
108  g_AllPeriods->GetXaxis()->CenterTitle();
109  g_AllPeriods->GetYaxis()->CenterTitle();
110  g_AllPeriods->GetXaxis()->SetTitle("Slices / 10^{9} POT");
111  g_AllPeriods->GetYaxis()->SetTitle("Average POT / 10^{12}");
112 
113  f_Out->cd();
114  for(auto spectrum : map_PeriodAndVarToSpectrum)
115  {
116  TH2 *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
117  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)));
118  h_Prof->SetMarkerStyle(datasets[spectrum.first.first].fStyle);
119  h_Prof->SetMarkerColor(datasets[spectrum.first.first].fColor);
120  h_Prof->SetLineColor (datasets[spectrum.first.first].fColor);
121  h_Prof->SetMarkerSize (1.4);
122  if(spectrum.first.first==0)
123  {
124  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
125  h_Prof->SetTitle("");
126  h_Prof->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(0));
127  h_Prof->GetYaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
128  h_Prof->GetXaxis()->CenterTitle();
129  h_Prof->GetYaxis()->CenterTitle();
130  h_Prof->Draw();
131  }
132 
133  TGraphErrors *g = new TGraphErrors;
134  g->SetName((TString)("1D_"+vec_Plots.at(spectrum.first.second).fName+"_"+Form("%i_%i", datasets[spectrum.first.first].fYear, datasets[spectrum.first.first].fMonth)));
135  double y = h_2D->Integral();
136  double dy = std::sqrt(y);
137  double x = vec_hPOT.at(spectrum.first.first) ->GetMean();
138  double npot = vec_hNumPOT.at(spectrum.first.first)->Integral()/1000.;
139  g->SetTitle("");
140  g->SetMarkerStyle(datasets[spectrum.first.first].fStyle);
141  g->SetMarkerColor(datasets[spectrum.first.first].fColor);
142  g->SetLineColor (datasets[spectrum.first.first].fColor);
143  if(x > 0 && npot > 0)
144  {
145  g->SetPoint (0, x, y/npot);
146  g->SetPointError(0, 0, dy/npot);
147 
148  g_AllPeriods->SetPoint (g_AllPeriods->GetN(), x, y/npot );
149  g_AllPeriods->SetPointError(g_AllPeriods->GetN()-1, 0, dy/npot);
150  }
151 
152  if(spectrum.first.first==0)
153  {
154  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);
155  h_Canvas->GetXaxis()->CenterTitle();
156  h_Canvas->GetYaxis()->CenterTitle();
157  h_Canvas->GetYaxis()->SetTitle("Slices / 10^{9} POT");
158  h_Canvas->GetXaxis()->SetTitle("Average POT / 10^{12}");
159  h_Canvas->SetLineColor(0);
160  h_Canvas->SetLineWidth(0);
161 
162  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_AllPeriods"]->cd();
163  h_Canvas->Draw();
164  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
165  h_Canvas->Draw();
166 
167  g->Draw("P");
168  }
169  else
170  {
171  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
172  g->Draw("P");
173  }
174  if(spectrum.first.second==0)
175  {
176  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
177  leg->AddEntry(g, (TString)datasets[spectrum.first.first].fLabel, "P");
178  }
179  if(spectrum.first.first==datasets.size()-1)
180  {
181  leg->Draw();
182  lastUpdated(endTime);
183  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_AllPeriods"]->cd();
184  g_AllPeriods->Draw("P");
185  g_AllPeriods->Fit("pol1");
186  lastUpdated(endTime);
187  }
188  }
189 
190  for(auto it : map_VarToCanvas)
191  {
192  it.second->SaveAs((TString)(outDir+it.first+"_IntensityEffect.png"));
193  it.second->Write();
194  it.second->Clear();
195  }
196 
197  std::cout << "ALL PLOTS MADE" << std::endl;
198 }
const SpillVar kPOTe12([](const caf::SRSpillProxy *spill){return spill->spillpot/1e12;})
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetSpillCut(const SpillCut &cut)
std::string outDir
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const SpillVar kSpillUnweighted
Definition: Var.h:98
void lastUpdated(unsigned int timestamp)
double dy[NP][NC]
virtual void AddSpillHistogram(TH1 *h, const SpillVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
Uses include counting the total POT or spills in a run.
SpillCut getYearAndMonth(unsigned int iYear, unsigned int iMonth)
Definition: VarsAndCuts.h:71
virtual void Go() override
Load all the registered spectra.
const double j
Definition: BetheBloch.cxx:29
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
OStream cout
Definition: OStream.cxx:6
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const Var kYear
Definition: VarsAndCuts.h:34
std::vector< Dataset > datasets
Definition: Periods.h:3
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
void drawIntensityEffect_2019(std::string samdef, unsigned int endTime, unsigned int nDays, std::string outDir)
const Var kMonth
Definition: VarsAndCuts.h:35
std::vector< Plot > vec_Plots
enum BeamMode string