drawLongTerm_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/NumuVars.h"
10 #include "CAFAna/Vars/Vars.h"
11 #include "CAFAna/Analysis/Calcs.h"
12 #include "CAFAna/Analysis/Plots.h"
13 
14 #include "TStyle.h"
15 #include "TCanvas.h"
16 #include "TGraph.h"
17 #include "TGraphErrors.h"
18 #include "TGraphAsymmErrors.h"
19 #include "TFile.h"
20 #include "TH1.h"
21 #include "TH2.h"
22 #include "TLegend.h"
23 #include "TLatex.h"
24 #include "TProfile.h"
25 #include "TSpectrum.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 drawLongTerm_2019(std::string samDef, unsigned int endTime, unsigned int nDays, std::string outDir)
34 {
35  //SET SIG FIGS.
36  gStyle->SetStatFormat("6.3g");
37  gStyle->SetPaintTextFormat("6.3g");
38  gStyle->SetFitFormat("6.3g");
39 
40  //DEFINE A TFILE TO WRITE ALL CANVASES TO.
41  TFile *f_Out = new TFile((TString)(outDir+"BeamMonitoring_LongTerm.root"),"RECREATE");
42  TDirectory *dir_Profile = f_Out->mkdir("Profiles");
43  TDirectory *dir_1D = f_Out->mkdir("1Ds");
44  TDirectory *dir_Ratio = f_Out->mkdir("Ratios");
45 
46  //APPLY TIME CUTS TO DATASET, DEFINE THE LOADER AND SET SPILL CUTS.
47  unsigned int startTime = endTime - 86400*nDays;
48  const std::string fFiles = "dataset_def_name_newest_snapshot " + samDef + Form(" and online.subrunendtime >= '%u' and online.subrunendtime <= '%u'", startTime, endTime);
49  SpectrumLoader loader(fFiles);
51 
52  const Binning kTimeBinning = Binning::Simple(nDays, startTime, endTime);
53  const Cut kNumuNDCut = kNumuCutND2018;
54 
55  //PLOT SPILL TIME FOR A PARTICULAR PERIOD SO THAT WE HAVE A POT TO NORMALISE BY.
56  //CUT USING getYearAndMonth.
57  std::vector<TH1D*> vec_hPOT;
58 
59  for(int i = 0; i < kNumPeriods; i++)
60  {
61  TH1D *hPOT = new TH1D(Form("hPOT_Period%i", i),";Central Time;POT (10^{15})", nDays, startTime, endTime);
62  vec_hPOT.push_back(hPOT);
63 
64  loader.AddSpillHistogram(vec_hPOT.at(i), kSpillTime, getYearAndMonth(datasets[i].fYear, datasets[i].fMonth), kPOTe15);
65  }
66 
67  //DEFINE PLOTS OF THE VARIABLES WE WANT TO SHOW.
68  std::vector<Plot> vec_Plots = {
69  {"Central Time", "Reconstructed Neutrino Energy, (GeV)", "ccE", Binning::Simple(10,0,5), kCCE, kNumuNDCut, 1.5, 2.2 },
70  {"Central Time", "Reconstructed Hadronic Energy, (GeV)", "hadE", Binning::Simple(12,0,3), kHadE, kNumuNDCut, 0.2, 0.6 },
71  {"Central Time", "Reconstructed Muon Energy, (GeV)", "muonE", Binning::Simple(15,0,4.5), kMuE, kNumuNDCut, 0.9, 1.7 },
72  {"Central Time", "Hadronic Energy Fraction", "hadEfrac", Binning::Simple(10,0,1.0), kHadEFrac, kNumuNDCut, 0.05, 0.35},
73  {"Central Time", "No. of Hits", "nhit", Binning::Simple(10,0,250), kNHit, kNumuNDCut, 85, 105 },
74  {"Central Time", "Vertex X, (m)", "startX", Binning::Simple(10,-2,2), kTrkStartX, kNumuNDCut, -0.40, 0.0 },
75  {"Central Time", "Vertex Y, (m)", "startY", Binning::Simple(10,-2,2), kTrkStartY, kNumuNDCut, -0.25, 0.15},
76  {"Central Time", "Vertex Z, (m)", "startZ", Binning::Simple(12,0,12), kTrkStartZ, kNumuNDCut, 4.5, 5.8 }
77  };
78 
79  std::map<std::pair<unsigned int,unsigned int>,Spectrum*> map_PeriodAndVarToSpectrum;
80  std::map<std::string,TCanvas*> map_VarToCanvas;
81  std::map<std::string,TH1D*> map_VarToMeanHist;
82 
83  for(unsigned int i = 0; i < vec_hPOT.size(); i++)
84  {
85  for(unsigned int j = 0; j < vec_Plots.size(); j++)
86  {
87  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,
88  kYear == datasets[i].fYear && kMonth == datasets[i].fMonth && vec_Plots.at(j).fCut);
89  map_PeriodAndVarToSpectrum[{i,j}] = h_Spectrum;
90  if(i==0)
91  {
92  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"));
93  map_VarToCanvas[vec_Plots.at(j).fName]->SetBottomMargin(0.15);
94  map_VarToCanvas[vec_Plots.at(j).fName]->SetLeftMargin (0.15);
95  map_VarToCanvas[vec_Plots.at(j).fName]->SetRightMargin (0.15);
96  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(),
97  map_VarToCanvas[vec_Plots.at(j).fName]->GetWh());
98  map_VarToCanvas[vec_Plots.at(j).fName+"_1D"]->SetBottomMargin(0.15);
99  map_VarToCanvas[vec_Plots.at(j).fName+"_1D"]->SetLeftMargin (0.15);
100  map_VarToCanvas[vec_Plots.at(j).fName+"_1D"]->SetRightMargin (0.15);
101  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(),
102  map_VarToCanvas[vec_Plots.at(j).fName]->GetWh());
103  map_VarToCanvas[vec_Plots.at(j).fName+"_Ratio"]->SetBottomMargin(0.15);
104  map_VarToCanvas[vec_Plots.at(j).fName+"_Ratio"]->SetLeftMargin (0.15);
105  map_VarToCanvas[vec_Plots.at(j).fName+"_Ratio"]->SetRightMargin (0.15);
106 
107  TH1D *h_Mean = new TH1D();
108  h_Mean->SetName((TString)vec_Plots.at(j).fName);
109  h_Mean->SetLineWidth(1);
110  h_Mean->SetLineColor(kBlack);
111  map_VarToMeanHist[vec_Plots.at(j).fName] = h_Mean;
112  }
113  }
114  }
115 
116  loader.Go();
117 
118  //PROFILES.
119  TLegend* leg = new TLegend(0.17,0.17,0.83,0.4);
120  leg->SetFillColor(0);
121  leg->SetNColumns(3);
122  TLegend* leg_1D_1 = new TLegend(0.6,0.6,0.85,0.85);
123  leg_1D_1->SetFillColor(0);
124  leg_1D_1->SetNColumns(3);
125  TLegend* leg_1D_2 = new TLegend(0.35,0.20,0.62,0.48);
126  leg_1D_2->SetFillColor(0);
127  leg_1D_2->SetNColumns(3);
128 
129  for(auto spectrum : map_PeriodAndVarToSpectrum)
130  {
131  TH2 *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
132  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)));
133  h_Prof->SetMarkerStyle(datasets[spectrum.first.first].fStyle);
134  h_Prof->SetMarkerColor(datasets[spectrum.first.first].fColor);
135  h_Prof->SetLineColor (datasets[spectrum.first.first].fColor);
136  h_Prof->SetMarkerSize (1.4);
137  if(spectrum.first.first==0)
138  {
139  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
140  h_Prof->SetTitle("");
141  h_Prof->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(0));
142  h_Prof->GetYaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
143  h_Prof->GetXaxis()->CenterTitle();
144  h_Prof->GetYaxis()->CenterTitle();
145  h_Prof->GetXaxis()->SetTimeDisplay(1);
146  h_Prof->GetXaxis()->SetTimeFormat("%y-%m-%d%F1970-01-01 00:00:00");
147  h_Prof->Draw();
148 
149  unsigned int minBin = h_Prof->GetXaxis()->FindBin(startTime);
150  unsigned int maxBin = h_Prof->GetXaxis()->FindBin(endTime);
151  double minZoom = h_Prof->GetXaxis()->GetBinLowEdge(minBin-1);
152  double maxZoom = h_Prof->GetXaxis()->GetBinLowEdge(maxBin+2);
153  h_Prof->GetXaxis()->SetRangeUser(minZoom, maxZoom);
154  h_Prof->GetYaxis()->SetRangeUser(vec_Plots.at(spectrum.first.second).fYRangeLow, vec_Plots.at(spectrum.first.second).fYRangeHigh);
155  }
156  else
157  {
158  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
159  h_Prof->Draw("SAME");
160  }
161 
162  if(spectrum.first.second==0)
163  {
164  leg->AddEntry(h_Prof, (TString)datasets[spectrum.first.first].fLabel, "P");
165  leg_1D_1->AddEntry(h_Prof, (TString)datasets[spectrum.first.first].fLabel, "P");
166  leg_1D_2->AddEntry(h_Prof, (TString)datasets[spectrum.first.first].fLabel, "P");
167  }
168  if(spectrum.first.first==datasets.size()-1)
169  {
170  leg->Draw();
171  lastUpdated(endTime);
172  }
173  }
174  TH1D *h_Dummy = new TH1D("h_Dummy","h_Dummy",1,0,1);
175  h_Dummy->SetLineColor(kBlack);
176  leg_1D_1->AddEntry(h_Dummy, "Average (Mean)", "L");
177  leg_1D_2->AddEntry(h_Dummy, "Average (Mean)", "L");
178 
179  //1D.
180  for(auto spectrum : map_PeriodAndVarToSpectrum)
181  {
182  TH2 *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
183  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)));
184  if(spectrum.first.first==0)
185  {
186  map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName] = (TH1D*)h_1D->Clone((TString)("h_Mean"+vec_Plots.at(spectrum.first.second).fName));
187  }
188  else
189  {
190  map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Add(h_1D);
191  }
192 
193  if(vec_hPOT.at(spectrum.first.first)->Integral()!=0.){h_1D->Scale(1.0/vec_hPOT.at(spectrum.first.first)->Integral());}
194  TGraphErrors *gr = new TGraphErrors(h_1D);
195  gr->SetMarkerStyle(datasets[spectrum.first.first].fStyle);
196  gr->SetMarkerColor(datasets[spectrum.first.first].fColor);
197  gr->SetLineColor (datasets[spectrum.first.first].fColor);
198 
199  for(int j = 0; j < gr->GetN(); j++)
200  {
201  double x, y, dy;
202  double shift = (h_1D->GetBinWidth(1)/(double)datasets.size())*((spectrum.first.first+1)*1. - datasets.size()*0.5);
203  gr->GetPoint(j, x, y);
204  dy = gr->GetErrorY(j);
205  gr->SetPoint(j, x + shift, y);
206  gr->SetPointError(j, 0, dy);
207  }
208  if(spectrum.first.first==0)
209  {
210  gr->SetTitle("");
211  gr->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
212  gr->GetYaxis()->SetTitle("Slices / 10^{15} POT");
213  gr->GetXaxis()->CenterTitle();
214  gr->GetYaxis()->CenterTitle();
215  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_1D"]->cd();
216  gr->Draw("AP");
217  }
218  else
219  {
220  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_1D"]->cd();
221  gr->Draw("P");
222  }
223 
224  if(spectrum.first.first==datasets.size()-1)
225  {
226  map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Scale(1.0e15/map_PeriodAndVarToSpectrum[{0, spectrum.first.second}]->POT());
227  map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->SetTitle("");
228  map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Draw("HIST SAME");
229  if(vec_Plots.at(spectrum.first.second).fName=="ccE" ||
230  vec_Plots.at(spectrum.first.second).fName=="muonE" ||
231  vec_Plots.at(spectrum.first.second).fName=="hadE" ||
232  vec_Plots.at(spectrum.first.second).fName=="hadEfrac" ||
233  vec_Plots.at(spectrum.first.second).fName=="nhit") {
234  leg_1D_1->Draw();
235  }
236  else{
237  leg_1D_2->Draw();
238  }
239 
240  lastUpdated(endTime);
241  }
242  }
243  //RATIOS.
244  for(auto spectrum : map_PeriodAndVarToSpectrum)
245  {
246  TH2 *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
247  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)));
248  if(vec_hPOT.at(spectrum.first.first)->Integral()!=0.){h_1D->Scale(1.0/vec_hPOT.at(spectrum.first.first)->Integral());}
249  TGraphErrors *gr = new TGraphErrors(h_1D);
250  gr->SetMarkerStyle(datasets[spectrum.first.first].fStyle);
251  gr->SetMarkerColor(datasets[spectrum.first.first].fColor);
252  gr->SetLineColor (datasets[spectrum.first.first].fColor);
253 
254  for(int j = 0; j < gr->GetN(); j++)
255  {
256  double x, y, my, dy;
257  double shift = (h_1D->GetBinWidth(1)/(double)datasets.size())*((spectrum.first.first+1)*1. - datasets.size()*0.5);
258  my = map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->GetBinContent(j+1);
259  gr->GetPoint(j, x, y);
260  if(my>0)
261  {
262  y /= my;
263  dy = gr->GetErrorY(j)/my;
264  gr->SetPoint(j, x + shift, y);
265  gr->SetPointError(j, 0, dy);
266  }
267  }
268  double xMin, xMax;
269  if(spectrum.first.first==0)
270  {
271  map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName] = (TH1D*)h_1D->Clone((TString)("h_Mean"+vec_Plots.at(spectrum.first.second).fName));
272  gr->SetTitle("");
273  gr->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
274  gr->GetYaxis()->SetTitle("Slices per 10^{15} POT / Average");
275  gr->GetXaxis()->CenterTitle();
276  gr->GetYaxis()->CenterTitle();
277  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_Ratio"]->cd();
278  gr->Draw("AP");
279  gr->GetYaxis()->SetRangeUser(0.8,1.2);
280 
281  double binsArrary[h_1D->GetNbinsX()];
282  h_1D->GetLowEdge(binsArrary);
283  for(int j = 0; j < h_1D->GetNbinsX(); j++)
284  {
285  TLine *l = new TLine(binsArrary[j], 0.8, binsArrary[j], 1.2);
286  l->SetLineStyle(kDashed);
287  l->Draw();
288  }
289  }
290  else
291  {
292  map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+"_Ratio"]->cd();
293  gr->Draw("P");
294  }
295 
296  if(spectrum.first.first==datasets.size()-1)
297  {
298  lastUpdated(endTime);
299  TGraph *xAxis = new TGraph;
300  xAxis->SetPoint(0, -1.e99, 1);
301  xAxis->SetPoint(1, 1.e99, 1);
302  xAxis->Draw("L");
303  }
304  }
305 
306  for(auto it : map_VarToCanvas)
307  {
308  it.second->SaveAs((TString)(outDir+it.first+"_LongTerm.png"));
309  if (((std::string)(it.second->GetName())).find("Profile")!=std::string::npos){dir_Profile->cd();}
310  else if(((std::string)(it.second->GetName())).find("1D") !=std::string::npos){dir_1D ->cd();}
311  else if(((std::string)(it.second->GetName())).find("Ratio") !=std::string::npos){dir_Ratio ->cd();}
312  it.second->Write();
313  it.second->Clear();
314  }
315 
316  std::cout << "ALL PLOTS MADE" << std::endl;
317 
318 }
const Var kHadE
Definition: NumuVars.h:23
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 SpillVar kSpillTime([](const caf::SRSpillProxy *spill){return spill->spilltimesec;})
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
const Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
Definition: NumuVars.h:52
void SetSpillCut(const SpillCut &cut)
std::string outDir
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kTime
Definition: VarsAndCuts.h:33
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.
const int kNumPeriods
Definition: Periods.h:8
const Var kHadEFrac
Definition: NumuVars.h:24
const Var kNHit
Definition: Vars.cxx:71
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
Definition: NumuVars.h:53
const Var kCCE
Definition: NumuVars.h:21
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
const Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
Definition: NumuVars.h:51
OStream cout
Definition: OStream.cxx:6
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const SpillVar kPOTe15([](const caf::SRSpillProxy *spill){return spill->spillpot/1e15;})
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const Var kYear
Definition: VarsAndCuts.h:34
const Var kMuE
Definition: NumuVars.h:22
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
void drawLongTerm_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