plotStabilitySpectra.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 #include <sstream>
14 
15 #include "TStyle.h"
16 #include "TCanvas.h"
17 #include "TGraph.h"
18 #include "TGraphErrors.h"
19 #include "TGraphAsymmErrors.h"
20 #include "TFile.h"
21 #include "TH1.h"
22 #include "TH2.h"
23 #include "TKey.h"
24 #include "TLegend.h"
25 #include "TLatex.h"
26 #include "TProfile.h"
27 #include "TSpectrum.h"
28 
29 #include "VarsAndCuts_2019.h"
30 
31 using namespace ana;
32 
33 void lastUpdated(unsigned int timestamp){
34  TTimeStamp refreshTime(timestamp + 604800);
35  TTimeStamp includedTime(timestamp);
36  unsigned int year, month, day, hour, minute, second;
37 
38  refreshTime.GetDate(kFALSE, 0, &year, &month, &day);
39  refreshTime.GetTime(kFALSE, 0, &hour, &minute, &second);
40 
41  TPaveLabel* pl = new TPaveLabel;
42  pl->SetX1NDC(0.50);
43  pl->SetX2NDC(0.98);
44  pl->SetY1NDC(0.95);
45  pl->SetY2NDC(1.00);
46  pl->SetFillColor(0);
47  pl->SetShadowColor(0);
48  pl->SetLineColor(0);
49  pl->SetLabel(Form("Last refreshed: %u-%02u-%02u-%02u:%02u:%02u",year,month,day,hour,minute,second));
50  pl->SetTextSize(0);
51  pl->Draw();
52 
53  includedTime.GetDate(kFALSE, 0, &year, &month, &day);
54  includedTime.GetTime(kFALSE, 0, &hour, &minute, &second);
55 
56  TPaveLabel* pl2 = new TPaveLabel;
57  pl2->SetX1NDC(0.10);
58  pl2->SetX2NDC(0.35);
59  pl2->SetY1NDC(0.95);
60  pl2->SetY2NDC(1.00);
61  pl2->SetFillColor(0);
62  pl2->SetShadowColor(0);
63  pl2->SetLineColor(0);
64  pl2->SetLabel(Form("Data up to %u-%02u-%02u-%02u:%02u:%02u",year,month,day,hour,minute,second));
65  pl2->SetTextSize(0);
66  pl2->Draw();
67 }
68 
69 std::pair<double,double> getRange(std::string const &s_VarName)
70 {
71  int thePlot = 0;
72  for(unsigned int i = 0; i < vec_Plots.size(); i++)
73  {
74  if(vec_Plots.at(i).fName==s_VarName)
75  {
76  thePlot = i;
77  break;
78  }
79  }
80  return {vec_Plots.at(thePlot).fYRangeLow, vec_Plots.at(thePlot).fYRangeHigh};
81 }
82 
83 TLegend* getLegend(std::string const &s_Var, std::vector<TLegend*> const &vec_Leg)
84 {
85  if(s_Var=="ccE" || s_Var=="muonE" || s_Var=="hadE" || s_Var=="hadEfrac" || s_Var=="nhit")
86  {
87  return vec_Leg.at(1);
88  }
89  else
90  {
91  return vec_Leg.at(2);
92  }
93 }
94 
95 std::vector<TString> parseName(TString const &s_Name)
96 {
97  std::stringstream ss_Name((std::string)s_Name);
98  std::string s_Comp;
99  std::vector<TString> vec_Comp;
100  while(std::getline(ss_Name, s_Comp, '_'))
101  {
102  vec_Comp.push_back((TString)s_Comp);
103  }
104 
105  return vec_Comp;
106 }
107 
108 double getPOT(std::string const &s_TH2Name, std::vector<TH1D*> const &vec_hPOT)
109 {
110  double POT = 0.;
111 
112  for(unsigned int i = 0; i < vec_hPOT.size(); i++)
113  {
114  std::string s_POTName = ((std::string)vec_hPOT.at(i)->GetName()).substr(5);
115  if(s_TH2Name.find(s_POTName)!=std::string::npos)
116  {
117  std::string s_TheChar = s_TH2Name.substr(s_POTName.length(), 1);
118  if(s_TheChar=="_")
119  {
120  return vec_hPOT.at(i)->Integral();
121  }
122  }
123  }
124 
125  return POT;
126 }
127 
128 void removeZeroPOT(std::map<std::string,std::vector<TH2*>> &map_VarToTH2s, std::vector<TH1D*> const &vec_hPOT)
129 {
130  for(auto it : map_VarToTH2s)
131  {
132  std::vector<TH2*> vec_Temp;
133  for(unsigned int i = 0; i < it.second.size(); i++)
134  {
135  //CUT OUT LOW STATS MONTHS FROM THE CALCULATION OF THE MEANS.
136  if(getPOT((std::string)it.second.at(i)->GetName(), vec_hPOT)>1000.)
137  {
138  vec_Temp.push_back(it.second.at(i));
139  }
140  }
141  map_VarToTH2s[it.first] = vec_Temp;
142  }
143 
144  return;
145 }
146 
147 std::vector<TLegend*> makeProfiles(TFile *f_Out, std::string const &s_OutDir, std::map<std::string,std::vector<TH2*>> const &map_VarToTH2s, double const &endTime)
148 {
149  f_Out->cd();
150  TLegend *leg_Pr = new TLegend(0.17, 0.17, 0.83, 0.4); leg_Pr->SetFillStyle(0); leg_Pr->SetNColumns(3);
151  TLegend *leg_UR = new TLegend(0.6, 0.6, 0.85, 0.85); leg_UR->SetFillStyle(0); leg_UR->SetNColumns(3);
152  TLegend *leg_CM = new TLegend(0.35, 0.2, 0.62, 0.48); leg_CM->SetFillStyle(0); leg_CM->SetNColumns(3);
153 
154  bool makeLegends = true;
155  for(auto it : map_VarToTH2s)
156  {
157  TCanvas *c = new TCanvas(("c_prof_"+it.first).c_str(), ("c_prof_"+it.first).c_str());
158  for(unsigned int i = 0; i < it.second.size(); i++)
159  {
160  std::vector<TString> vec_Comp = parseName(it.second.at(i)->GetName());
161  TProfile *h_Prof = it.second.at(i)->ProfileX((TString)("prof_+"+it.first+"_")+vec_Comp.at(0)+"_"+vec_Comp.at(1));
162  h_Prof->SetMarkerStyle(std::stoi((std::string)vec_Comp.at(2)));
163  h_Prof->SetMarkerColor(std::stoi((std::string)vec_Comp.at(3)));
164  h_Prof->SetLineColor (std::stoi((std::string)vec_Comp.at(3)));
165  h_Prof->SetMarkerSize(1.4);
166  if(i == 0)
167  {
168  h_Prof->SetTitle("");
169  h_Prof->GetXaxis()->SetTitle(it.second.at(i)->GetXaxis()->GetTitle());
170  h_Prof->GetYaxis()->SetTitle(it.second.at(i)->GetYaxis()->GetTitle());
171  h_Prof->GetXaxis()->CenterTitle();
172  h_Prof->GetYaxis()->CenterTitle();
173  h_Prof->GetXaxis()->SetTimeDisplay(1);
174  h_Prof->GetXaxis()->SetTimeFormat("%y-%m-%d%F1970-01-01 00:00:00");
175  h_Prof->Draw();
176  h_Prof->GetYaxis()->SetRangeUser(getRange(it.first).first, getRange(it.first).second);
177  }
178  else
179  {
180  h_Prof->Draw("SAME");
181  }
182  if(makeLegends)
183  {
184  leg_Pr->AddEntry(h_Prof, it.second.at(i)->GetTitle(), "P");
185  leg_UR->AddEntry(h_Prof, it.second.at(i)->GetTitle(), "P");
186  leg_CM->AddEntry(h_Prof, it.second.at(i)->GetTitle(), "P");
187  }
188  if(i==it.second.size()-1)
189  {
190  leg_Pr->Draw();
191  }
192  }
193  makeLegends = false;
194 
195  lastUpdated(endTime);
196 
197  c->SaveAs((s_OutDir+"c_prof_"+it.first+".pdf").c_str());
198  c->SaveAs((s_OutDir+"c_prof_"+it.first+".png").c_str());
199  c->Write();
200  }
201 
202  return{leg_Pr, leg_UR, leg_CM};
203 }
204 
205 void make1DsAndRatio(TFile *f_Out, std::string const &s_OutDir, std::map<std::string,std::vector<TH2*>> const &map_VarToTH2s, std::vector<TH1D*> const &vec_hPOT, std::vector<TLegend*> const & vec_Leg, double const &endTime)
206 {
207  f_Out->cd();
208 
209  for(auto it : map_VarToTH2s)
210  {
211  TCanvas *c = new TCanvas(("c_1D_"+it.first).c_str(), ("c_1D_"+it.first).c_str(), 1000, 500);
212  TH1D *h_Mean;
213  unsigned int nNonZero = 0;
214  std::vector<TGraph*> vec_Graph;
215  for(unsigned int i = 0; i < it.second.size(); i++)
216  {
217  std::vector<TString> vec_Comp = parseName(it.second.at(i)->GetName());
218  TH1D *h_1D = it.second.at(i)->ProjectionY((TString)("1D_+"+it.first+"_")+vec_Comp.at(0)+"_"+vec_Comp.at(1));
219 
220  h_1D->Scale(1./getPOT((std::string)it.second.at(i)->GetName(), vec_hPOT));
221  TGraphErrors *g = new TGraphErrors(h_1D);
222  g->SetMarkerStyle(std::stoi((std::string)vec_Comp.at(2)));
223  g->SetMarkerColor(std::stoi((std::string)vec_Comp.at(3)));
224  g->SetLineColor (std::stoi((std::string)vec_Comp.at(3)));
225 
226  for(int j = 0; j < g->GetN(); j++)
227  {
228  double x, y, dy;
229  double fracShift = h_1D->GetBinWidth(1)/((double)it.second.size()+1.);
230  g->GetPoint(j, x, y);
231  dy = g->GetErrorY(j);
232  g->SetPoint(j, x - 0.5*h_1D->GetBinWidth(1) + (i+1)*fracShift, y);
233  g->SetPointError(j, 0, dy);
234  }
235 
236  if(i == 0)
237  {
238  h_Mean = (TH1D*)h_1D->Clone(("h_Mean_"+it.first).c_str());
239  nNonZero++;
240 
241  g->SetTitle("");
242  g->GetXaxis()->SetTitle(it.second.at(i)->GetYaxis()->GetTitle());
243  g->GetYaxis()->SetTitle("Slices / 10^{15} POT");
244  g->GetXaxis()->CenterTitle();
245  g->GetYaxis()->CenterTitle();
246  g->Draw("AP");
247 
248  vec_Graph.push_back(g);
249  }
250  else
251  {
252  if(getPOT((std::string)it.second.at(i)->GetName(), vec_hPOT)!=0)
253  {
254  h_Mean->Add(h_1D);
255  nNonZero++;
256  vec_Graph.push_back(g);
257  }
258  g->Draw("P");
259  }
260  if(i==it.second.size()-1)
261  {
262  h_Mean->Scale(1./(double)nNonZero);
263  h_Mean->Draw("HIST SAME");
264  getLegend(it.first, vec_Leg)->Draw();
265  }
266  }
267 
268  lastUpdated(endTime);
269 
270  c->SaveAs((s_OutDir+"c_1D_"+it.first+".pdf").c_str());
271  c->SaveAs((s_OutDir+"c_1D_"+it.first+".png").c_str());
272  c->Write();
273 
274  TCanvas *c_Rat = new TCanvas(("c_Rat_"+it.first).c_str(), ("c_Rat_"+it.first).c_str(), 1000, 500);
275  std::vector<TGraphErrors*> vec_GraphCopy;
276  for(unsigned int i = 0; i < vec_Graph.size(); i++)
277  {
278  vec_GraphCopy.push_back((TGraphErrors*)vec_Graph.at(i)->Clone());
279  for(int j = 0; j < vec_Graph.at(i)->GetN(); j++)
280  {
281  double x, y, ey;
282  vec_Graph.at(i)->GetPoint(j, x, y);
283  ey = vec_Graph.at(i)->GetErrorY(j);
284 
285  vec_GraphCopy.at(i)->SetPoint (j, x, y/h_Mean->GetBinContent(h_Mean->FindBin(x)));
286  vec_GraphCopy.at(i)->SetPointError(j, 0, ey/h_Mean->GetBinContent(h_Mean->FindBin(x)));
287  }
288  }
289 
290  for(unsigned int i = 0; i < vec_Graph.size(); i++)
291  {
292  if(i==0)
293  {
294  vec_GraphCopy.at(i)->Draw("AP");
295  vec_GraphCopy.at(i)->GetYaxis()->SetTitleOffset(0.9);
296  }
297  else
298  {
299  vec_GraphCopy.at(i)->Draw("P");
300  }
301  }
302  vec_GraphCopy.at(0)->GetYaxis()->SetRangeUser(0.8,1.2);
303  double binsArrary[h_Mean->GetNbinsX()];
304  h_Mean->GetLowEdge(binsArrary);
305  for(int j = 0; j < h_Mean->GetNbinsX(); j++)
306  {
307  TLine *l = new TLine(binsArrary[j], 0.8, binsArrary[j], 1.2);
308  l->SetLineStyle(kDashed);
309  l->Draw();
310  }
311  TLine *l_xAxis = new TLine(h_Mean->GetBinLowEdge(1),1.,h_Mean->GetBinLowEdge(h_Mean->GetNbinsX()+1),1.);
312  l_xAxis->Draw();
313  lastUpdated(endTime);
314 
315  c_Rat->SaveAs((s_OutDir+"c_Rat_"+it.first+".pdf").c_str());
316  c_Rat->SaveAs((s_OutDir+"c_Rat_"+it.first+".png").c_str());
317  c_Rat->Write();
318  }
319 
320  return;
321 }
322 
323 void plotStabilitySpectra(std::string s_InFile, double endTime, std::string s_OutDir)
324 {
325  TFile *f_In = new TFile(s_InFile.c_str(), "READ");
326  TFile *f_Out = new TFile((s_OutDir+"BeamMonitoring_StabilityCanvases.root").c_str(), "RECREATE");
327 
328  //CONTAINER TO HOLD THE TH2s.
329  std::map<std::string,std::vector<TH2*>> map_VarToTH2s;
330  std::vector<TH1D*> vec_hPOT;
331 
332  TIter iter(f_In->GetListOfKeys());
333  TKey *key;
334  while((key=(TKey*)iter()))
335  {
336  TDirectory *dir = (TDirectory*)key->ReadObj();
337 
338  TIter iter_Dir(dir->GetListOfKeys());
339  TKey *key_Dir;
340  std::string s_Dir = (std::string)dir->GetName();
341  if(s_Dir.find("Hists_")!=std::string::npos)
342  {
343  std::string s_VarName = s_Dir.substr(6, s_Dir.length()-6);
344  std::cout << "\nEXTRACTING: " << s_VarName << std::endl;
345  std::cout << std::endl;
346  while((key_Dir=(TKey*)iter_Dir()))
347  {
348  TH2 *h = (TH2*)key_Dir->ReadObj();
349  map_VarToTH2s[s_VarName].push_back(h);
350  }
351  }
352  else if(s_Dir.find("POT")!=std::string::npos)
353  {
354  std::cout << "EXTRACTING POT: " << std::endl;
355  std::cout << std::endl;
356  while((key_Dir=(TKey*)iter_Dir()))
357  {
358  TH1D *h = (TH1D*)key_Dir->ReadObj();
359  vec_hPOT.push_back(h);
360  }
361  }
362  }
363 
364  std::vector<TLegend*> vec_Leg = makeProfiles(f_Out, s_OutDir, map_VarToTH2s, endTime);
365 
366  removeZeroPOT(map_VarToTH2s, vec_hPOT);
367  make1DsAndRatio(f_Out, s_OutDir, map_VarToTH2s, vec_hPOT, vec_Leg, endTime);
368 
369  std::cout << "ALL PLOTS MADE!" << std::endl;
370 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
TLegend * getLegend(std::string const &s_Var, std::vector< TLegend * > const &vec_Leg)
void removeZeroPOT(std::map< std::string, std::vector< TH2 * >> &map_VarToTH2s, std::vector< TH1D * > const &vec_hPOT)
void lastUpdated(unsigned int timestamp)
double dy[NP][NC]
std::vector< TLegend * > makeProfiles(TFile *f_Out, std::string const &s_OutDir, std::map< std::string, std::vector< TH2 * >> const &map_VarToTH2s, double const &endTime)
const double j
Definition: BetheBloch.cxx:29
void make1DsAndRatio(TFile *f_Out, std::string const &s_OutDir, std::map< std::string, std::vector< TH2 * >> const &map_VarToTH2s, std::vector< TH1D * > const &vec_hPOT, std::vector< TLegend * > const &vec_Leg, double const &endTime)
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
std::pair< double, double > getRange(std::string const &s_VarName)
TDirectory * dir
Definition: macro.C:5
void plotStabilitySpectra(std::string s_InFile, double endTime, std::string s_OutDir)
std::vector< TString > parseName(TString const &s_Name)
double getPOT(std::string const &s_TH2Name, std::vector< TH1D * > const &vec_hPOT)
static constexpr Double_t year
Definition: Munits.h:185
std::vector< Plot > vec_Plots
enum BeamMode string