drawLongTerm.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void drawLongTerm(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming){
3  std::cout << "Sorry, you must run in compiled mode" << std::endl;
4 }
5 #else
6 
7 #include "CAFAna/Cuts/Cuts.h"
8 #include "CAFAna/Cuts/SpillCuts.h"
9 #include "CAFAna/Cuts/TimingCuts.h"
10 #include "CAFAna/Cuts/NumuCuts2018.h"
11 #include "CAFAna/Core/Loaders.h"
12 #include "CAFAna/Core/SpectrumLoader.h"
13 #include "CAFAna/Core/Spectrum.h"
14 #include "CAFAna/Core/ReweightableSpectrum.h"
15 #include "CAFAna/Vars/NumuVars.h"
16 #include "CAFAna/Vars/Vars.h"
17 #include "CAFAna/Analysis/Calcs.h"
18 #include "CAFAna/Analysis/Plots.h"
19 
20 #include "TStyle.h"
21 #include "TCanvas.h"
22 #include "TGraph.h"
23 #include "TGraphErrors.h"
24 #include "TGraphAsymmErrors.h"
25 #include "TFile.h"
26 #include "TH1.h"
27 #include "TH2.h"
28 #include "TLegend.h"
29 #include "TLatex.h"
30 #include "TProfile.h"
31 #include "TSpectrum.h"
32 
33 #include "VarsAndCuts.h"
34 #include "Periods.h"
35 #include "CustomFunctions.h"
36 
37 using namespace ana;
38 
39 void drawLongTerm(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming)
40 {
41 
42  std::cout << "Drawing plots..." << std::endl;
43 
44  gStyle->SetStatFormat("6.3g"); // Significant figures. Fix if this is not satisfactory
45  gStyle->SetPaintTextFormat("6.3g");
46  gStyle->SetFitFormat("6.3g");
47 
48  // Output directory
49  const std::string outDir = "/nova/ana/users/novadq/BeamMonitoring/Plots/stability/";
50  // Loader for data
51  const std::string fFiles = "dataset_def_name_newest_snapshot " + samdef + Form(" and online.subrunendtime >= '%u' and online.subrunendtime <= '%u'", timestamp - 86400*nDays, timestamp);
52 
53  // int nBins = nDays/7 + 1; // weekly looks better than daily
54  int nBins = nDays;
55 
56  const Binning kTimeBinning = Binning::Simple(nBins, timestamp - 86400*nDays, timestamp);
57 
58  const Cut kNumuNDCut = kNumuCutND2018;
59 
60  SpectrumLoader loader(fFiles);
61  loader.SetSpillCut(kStandardSpillCuts);
62 
63  TH1D* hPOT[kNumPeriods];
64 
65  for(int j = 0; j < kNumPeriods; ++j)
66  {
67  hPOT[j] = new TH1D(Form("pot_%i",j+1),";Central Time;POT (10^{15})", nBins, timestamp - 86400*nDays, timestamp);
68 
69  loader.AddSpillHistogram(hPOT[j], kSpillTime, getYearAndMonth(datasets[j].year, datasets[j].month), kPOTe15);
70  }
71 
72  const int kNumPlots = 8;
73 
74  Plot plots[kNumPlots] = {
75  {";Central Time;Reconstructed neutrino energy (GeV)","ccE", Binning::Simple(10,0,5), kCCE, kNumuNDCut, 1.5, 2.0},
76  {";Central Time;Reconstructed hadronic energy (GeV)","hadE", Binning::Simple(12,0,3), kHadE, kNumuNDCut, 0.15, 0.50},
77  {";Central Time;Reconstructed muon energy (GeV)","muonE", Binning::Simple(15,0,4.5), kMuE, kNumuNDCut, 1.0, 2.0},
78  {";Central Time;Hadronic energy fraction","hadEfrac", Binning::Simple(10,0,1.0), kHadEFrac, kNumuNDCut, 0, 0.25},
79  {";Central Time;No. of hits","nhit", Binning::Simple(10,0,250), kNHit, kNumuNDCut, 90, 110},
80  {";Central Time;Vertex X (m)","startX", Binning::Simple(10,-2,2), kTrkStartX, kNumuNDCut, -0.40, 0.0 },
81  {";Central Time;Vertex Y (m)","startY", Binning::Simple(10,-2,2), kTrkStartY, kNumuNDCut, -0.25, 0.15},
82  {";Central Time;Vertex Z (m)","startZ", Binning::Simple(12,0,12), kTrkStartZ, kNumuNDCut, 5.0, 6.25}
83  };
84 
85  /// ====================================================================== ///
86 
87  Spectrum* hData[kNumPlots][kNumPeriods];
88 
89  for(int i = 0; i < kNumPlots; ++i){
90  for(int j = 0; j < kNumPeriods; ++j){
91  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);
92  } // end loop over j
93  } // end loop over i
94 
95  // GO!
96  loader.Go();
97 
98 for(int i = 0; i < kNumPlots; ++i){
99 
100  std::string fullName = plots[i].fname+"_"+naming;
101  TCanvas* c = new TCanvas("c",(fullName).c_str());
102  c->SetBottomMargin(0.15);
103  c->SetLeftMargin(0.15);
104  c->SetRightMargin(0.15);
105 
106  TLegend* leg = new TLegend(0.17,0.17,0.83,0.4);
107  leg->SetFillColor(0);
108  leg->SetNColumns(3);
109 
110  // First the profile
111  for(int j = 0; j < kNumPeriods; ++j){
112 
113  std::string monthName = plots[i].fname+Form("_period%i_",j+1)+naming;
114 
115  TH2* Histogram = hData[i][j]->ToTH2(hData[i][j]->POT());
116 
117  TProfile* Profile = Histogram->ProfileX((monthName+"Profile").c_str());
118  Profile->SetMarkerStyle(datasets[j].style);
119  Profile->SetMarkerColor(datasets[j].color);
120  Profile->SetLineColor(datasets[j].color);
121  Profile->SetMarkerSize(1.4);
122 
123  Profile->GetXaxis()->CenterTitle();
124  Profile->GetYaxis()->CenterTitle();
125  Profile->GetYaxis()->SetTitle(Histogram->GetYaxis()->GetTitle());
126  Profile->GetXaxis()->SetTimeDisplay(1);
127  Profile->GetXaxis()->SetTimeFormat("%y-%m-%d%F1970-01-01 00:00:00");
128 
129  unsigned int minBin = Profile->GetXaxis()->FindBin(timestamp - 86400*nDays);
130  unsigned int maxBin = Profile->GetXaxis()->FindBin(timestamp);
131 
132  double minZoom = Profile->GetXaxis()->GetBinLowEdge(minBin-1);
133  double maxZoom = Profile->GetXaxis()->GetBinLowEdge(maxBin+2);
134 
135  Profile->GetXaxis()->SetRangeUser(minZoom, maxZoom);
136 
137  if(j == 0)
138  {
139  Profile->GetYaxis()->SetRangeUser(plots[i].yRangeLow, plots[i].yRangeHigh);
140  Profile->Draw();
141  }
142 
143  else Profile->Draw("same");
144 
145  leg->AddEntry(Profile, datasets[j].label.c_str(), "p");
146 
147  } // End loop over j
148 
149  leg->Draw();
150  lastUpdated(timestamp);
151 
152  //c->SaveAs(("stability/"+fullName+".png").c_str());
153  //c->SaveAs(("stability/"+fullName+".pdf").c_str());
154  //c->SaveAs(("stability/"+fullName+".root").c_str());
155  c->SaveAs((outDir+fullName+".png").c_str());
156  c->SaveAs((outDir+fullName+".pdf").c_str());
157  c->SaveAs((outDir+fullName+".root").c_str());
158 
159  // And finally the 1D plots
160  double ww = c->GetWw();
161  double wh = c->GetWh();
162  c->Close();
163  TCanvas* c2 = new TCanvas("c2",(fullName).c_str(), 2*ww, wh);
164  c2->SetBottomMargin(0.15);
165  c2->SetLeftMargin(0.15);
166  c2->SetRightMargin(0.15);
167 
168  c2->SetWindowSize(2*ww, wh); // wide for spectra
169 
170  TLegend* leg2;
171  if(i>4) leg2 = new TLegend(0.35,0.20,0.62,0.48);
172  else leg2 = new TLegend(0.55,0.60,0.82,0.88);
173  leg2->SetFillColor(0);
174  leg2->SetNColumns(3);
175 
176  std::string dim1Name = plots[i].fname+"_1D_"+naming;
177 
178  TH1* hMean1D = new TH1();
179 
180  for(int j = 0; j < kNumPeriods; ++j){
181  std::string dim1NameMonth = plots[i].fname+Form("_period%i_1D_",j+1)+naming;
182 
183  TH2* Histogram = hData[i][j]->ToTH2(hData[i][j]->POT());
184 
185  TH1* Hist1D = Histogram->ProjectionY((dim1NameMonth).c_str());
186 
187  if(j == 0) hMean1D = (TH1*) Hist1D->Clone("Mean_1D");
188  else hMean1D->Add(Hist1D);
189 
190  Hist1D->Scale(1.0/hPOT[j]->Integral());
191 
192  TGraphErrors* Graph = new TGraphErrors(Hist1D);
193  Graph->SetName((dim1NameMonth+"_gr").c_str());
194  Graph->SetMarkerStyle(datasets[j].style);
195  Graph->SetMarkerColor(datasets[j].color);
196  Graph->SetLineColor(datasets[j].color);
197  Graph->GetXaxis()->SetTitle(Hist1D->GetXaxis()->GetTitle());
198  Graph->GetYaxis()->SetTitle("Slices / 10^{15} POT");
199  Graph->GetXaxis()->CenterTitle();
200  Graph->GetYaxis()->CenterTitle();
201 
202  for(int fPoint = 0; fPoint<Graph->GetN(); ++fPoint)
203  {
204  double x, y, dy;
205  double shift = (Hist1D->GetBinWidth(1) / kNumPeriods) * ((j+1)*1. - kNumPeriods*0.5);
206  Graph->GetPoint(fPoint, x, y);
207  dy = Graph->GetErrorY(fPoint);
208  Graph->SetPoint(fPoint, x + shift, y);
209  Graph->SetPointError(fPoint,0,dy);
210  }
211 
212  if(j == 0)
213  {
214  Hist1D->SetLineWidth(0);
215  Hist1D->SetLineColor(10);
216  Hist1D->SetMarkerStyle(0);
217  Hist1D->SetMarkerColor(0);
218  Hist1D->Draw("hist");
219  }
220 
221  Graph->Draw("P");
222 
223  leg2->AddEntry(Graph, datasets[j].label.c_str(), "p");
224 
225  } // End loop over j
226 
227  hMean1D->Scale(1.0e15/hData[i][0]->POT());
228  hMean1D->SetLineWidth(1);
229  hMean1D->SetLineColor(kBlack);
230  hMean1D->Draw("hist same");
231 
232  leg2->AddEntry(hMean1D,"Average","l");
233  leg2->Draw();
234  lastUpdated(timestamp);
235 
236  //c2->SaveAs(("stability/"+dim1Name+".png").c_str());
237  //c2->SaveAs(("stability/"+dim1Name+".pdf").c_str());
238  //c2->SaveAs(("stability/"+dim1Name+".root").c_str());
239  c2->SaveAs((outDir+dim1Name+".png").c_str());
240  c2->SaveAs((outDir+dim1Name+".pdf").c_str());
241  c2->SaveAs((outDir+dim1Name+".root").c_str());
242 
243  // Now the ratio
244  c2->Clear();
245 
246  std::string ratioName = plots[i].fname+"_ratio_"+naming;
247 
248  for(int j = 0; j < kNumPeriods; ++j){
249  std::string ratioNameMonth = plots[i].fname+Form("_period%i_ratio_",j+1)+naming;
250 
251  TH2* Histogram = hData[i][j]->ToTH2(hData[i][j]->POT());
252 
253  TH1* Hist1D = Histogram->ProjectionY((ratioNameMonth).c_str());
254 
255  Hist1D->Scale(1.0/hPOT[j]->Integral());
256 
257  TGraphErrors* Graph = new TGraphErrors(Hist1D);
258  Hist1D->Reset();
259  Graph->SetName((ratioNameMonth+"_gr").c_str());
260  Graph->SetMarkerStyle(datasets[j].style);
261  Graph->SetMarkerColor(datasets[j].color);
262  Graph->SetLineColor(datasets[j].color);
263  Graph->GetXaxis()->SetTitle(Hist1D->GetXaxis()->GetTitle());
264  Graph->GetYaxis()->SetTitle("Slices per 10^{15} POT / Average");
265  Graph->GetXaxis()->CenterTitle();
266  Graph->GetYaxis()->CenterTitle();
267 
268  for(int fPoint = 0; fPoint<Graph->GetN(); ++fPoint)
269  {
270  double x, y, my, dy;
271  double shift = (Hist1D->GetBinWidth(1) / kNumPeriods) * ((j+1)*1. - kNumPeriods*0.5);
272  my = hMean1D->GetBinContent(fPoint+1);
273  Graph->GetPoint(fPoint, x, y);
274  if(my > 0)
275  {
276  y = y/my;
277  dy = Graph->GetErrorY(fPoint)/my;
278  Graph->SetPoint(fPoint, x + shift, y);
279  Graph->SetPointError(fPoint,0,dy);
280  }
281  }
282 
283  if(j == 0)
284  {
285  Hist1D->GetYaxis()->SetRangeUser(0.8,1.2);
286  Hist1D->SetLineWidth(0);
287  Hist1D->SetLineColor(10);
288  Hist1D->SetMarkerStyle(0);
289  Hist1D->SetMarkerColor(0);
290  Hist1D->Draw("hist");
291  }
292 
293  Graph->Draw("P");
294  } // End loop over j
295 
296  TGraph* gr1 = new TGraph;
297  gr1->SetPoint(0,-1.e99,1.);
298  gr1->SetPoint(1, 1.e99,1.);
299  gr1->SetLineColor(kBlack);
300  gr1->SetLineWidth(1);
301  gr1->SetFillColor(10);
302  gr1->Draw("L");
303 
304  lastUpdated(timestamp);
305 
306  //c2->SaveAs(("stability/"+ratioName+".png").c_str());
307  //c2->SaveAs(("stability/"+ratioName+".pdf").c_str());
308  //c2->SaveAs(("stability/"+ratioName+".root").c_str());
309  c2->SaveAs((outDir+ratioName+".png").c_str());
310  c2->SaveAs((outDir+ratioName+".pdf").c_str());
311  c2->SaveAs((outDir+ratioName+".root").c_str());
312 
313 } // End loop over i
314 
315  std::cout << "All plots finished!" << std::endl;
316 
317 } // End of function
318 
319 #endif
ana::kHadE
const Var kHadE
Definition: NumuVars.h:23
plotROC.nBins
int nBins
Definition: plotROC.py:16
ana::Spectrum::ToTH2
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:226
PandAna.Demos.tute_pid_validation.var
list var
Definition: tute_pid_validation.py:17
ana::Binning
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
kSpillTime
const SpillVar kSpillTime([](const caf::SRSpillProxy *spill){return spill->spilltimesec;})
plot_validation_datamc_2018.color
color
Definition: plot_validation_datamc_2018.py:516
Integral
double Integral(const Spectrum &s, const double pot, int cvnregion)
Definition: syst_table_from_preds.C:46
ana::kNumuCutND2018
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
ana::kTrkStartY
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
ana::Profile
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
ana::SpectrumLoaderBase::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: SpectrumLoaderBase.cxx:64
gr1
TGraph * gr1
Definition: compare.C:42
outDir
std::string outDir
Definition: plot_ND_DataMC.C:35
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
submit_syst.y
y
Definition: submit_syst.py:249
kTime
const Var kTime
Definition: VarsAndCuts.h:33
label
const char * label
Definition: cana_time_peak.C:34
demo5.c2
c2
Definition: demo5.py:33
lastUpdated
void lastUpdated(unsigned int timestamp)
Definition: CustomFunctions.h:4
submit_syst.x
x
Definition: submit_syst.py:249
dy
double dy[NP][NC]
Definition: checkMeasured_Target.C:19
Plot
====================================================================== ///
Definition: CutFlow_Data.C:28
ana::SpectrumLoaderBase::AddSpillHistogram
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.
Definition: SpectrumLoaderBase.cxx:220
kNumPeriods
const int kNumPeriods
Definition: Periods.h:8
ana::kHadEFrac
const Var kHadEFrac
Definition: NumuVars.h:24
kNumPlots
const int kNumPlots
Definition: GetSpectra.h:1
ana::kNHit
const Var kNHit
Definition: Vars.cxx:71
plots
const std::vector< Plot > plots
Definition: make_DiF_plots.C:85
ana::kTrkStartZ
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
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
ana::kCCE
const Var kCCE
Definition: NumuVars.h:21
getYearAndMonth
SpillCut getYearAndMonth(unsigned int iYear, unsigned int iMonth)
Definition: VarsAndCuts.h:71
ana::SpectrumLoader::Go
virtual void Go() override
Load all the registered spectra.
Definition: SpectrumLoader.cxx:66
calib::j
const double j
Definition: BetheBloch.cxx:29
demo0.loader
loader
Definition: demo0.py:10
ana::kTrkStartX
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
om::cout
OStream cout
Definition: OStream.cxx:6
POT
std::vector< double > POT
Definition: neutKEsyst_plot.C:45
ana::bins
const Binning bins
Definition: NumuCC_CPiBin.h:8
Plot::fname
std::string fname
Definition: CutFlow_Data.C:30
xml_schema::string
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
cut
const Cut cut
Definition: exporter_fd.C:30
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
submit_nova_art.timestamp
timestamp
Definition: submit_nova_art.py:688
kPOTe15
const SpillVar kPOTe15([](const caf::SRSpillProxy *spill){return spill->spillpot/1e15;})
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
drawLongTerm
void drawLongTerm(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming)
Definition: drawLongTerm.C:39
kYear
const Var kYear
Definition: VarsAndCuts.h:34
ana::kMuE
const Var kMuE
Definition: NumuVars.h:22
datasets
std::vector< Dataset > datasets
Definition: Periods.h:3
ana::Binning::Simple
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
kMonth
const Var kMonth
Definition: VarsAndCuts.h:35
Munits::year
static constexpr Double_t year
Definition: Munits.h:185
Generated on Fri Sep 18 2020 05:13:49 for NOvA by   doxygen 1.8.11