DataMCNDAna_nus17.C
Go to the documentation of this file.
1 // A macro for Near Detector N-1 analysis
2 // cafe -bq DataMCNDAna_nus17.C <folder>
3 // where folder = the directory the input file lives
4 // and the output file will be saved too.
5 //
6 // Gavin S. Davies
7 
10 #include "CAFAna/Core/Spectrum.h"
13 
14 #include "TCanvas.h"
15 #include "TFile.h"
16 #include "TH1.h"
17 #include "TKey.h"
18 #include "TLegend.h"
19 
20 #include <string>
21 
22 using namespace ana;
23 
24 void make_dataMC(TFile* SavedFile,
25  TFile* OutFile);
26 
27 void make_plots(const ana::Spectrum& data,
28  const ana::CheatDecomp& MC,
29  std::string plotlabel,
30  TFile* plotsOut);
31 
33 
35 {
36  TH1::AddDirectory(0);
37  // load in loaded specs
38  std::string filein = folder + "/dataMCNDLoadOut.root";
39  std::string fileout = folder + "/dataMCNDAnaOut.root";
40 
41  TFile* spSaved = new TFile(filein.c_str());
42  TFile* plotsOut = new TFile(fileout.c_str(), "RECREATE");
43 
44  // make plots
45  make_dataMC(spSaved, plotsOut);
46 
47  spSaved->Close();
48  plotsOut->Close();
49  std::cout << "plotted all" << std::endl;
50 }
51 
52 //..............................................................................
53 void make_dataMC(TFile* SavedFile,
54  TFile* OutFile)
55 {
56  TIter next(SavedFile->GetListOfKeys());
57  TKey *key;
58  while ((key=(TKey*)next()))
59  {
60  const char* fDatKey = key->GetName();
61  std::string strDatKey(fDatKey);
62 
63  if(strDatKey.find("MC") != std::string::npos) { continue; }
64 
65  std::string strMCKey(fDatKey);
66  strMCKey.replace(strMCKey.size()-4, 4, "MC");
67 
68  std::string fPlotLabel(fDatKey);
69  fPlotLabel.erase(fPlotLabel.end()-4, fPlotLabel.end());
70 
71  std::cout << strDatKey << ", " << strMCKey << std::endl;
72  std::cout << fPlotLabel << std::endl;
73 
74  CheatDecomp sMC = *CheatDecomp::LoadFrom(SavedFile, strMCKey.c_str());
75  Spectrum sData = *Spectrum::LoadFrom(SavedFile, strDatKey.c_str());
76 
77  make_plots(sData, sMC, fPlotLabel.c_str(), OutFile);
78  }
79 };
80 
81 //..............................................................................
83  const ana::CheatDecomp& MC,
84  std::string plotlabel,
85  TFile* plotsOut)
86 {
87  Spectrum empty = MC.NCTotalComponent();
88  empty.Clear();
89 
90  Spectrum mcspecs[MAXSPEC] = {
91  MC.NCTotalComponent() +
92  MC.NueComponent() + MC.AntiNueComponent() +
93  MC.NumuComponent() + MC.AntiNumuComponent(),
94  MC.NCTotalComponent(),
95  MC.NumuComponent() + MC.AntiNumuComponent(),
96  MC.NueComponent() + MC.AntiNueComponent(),
97  empty,
98  empty
99  };
100 
101  double pot = data.POT();
102  std::cout << "data POT: " << pot << std::endl;
103  TH1* hDat = data .ToTH1(pot);
104  TH1* hMCTot = mcspecs[0].ToTH1(pot);
105  TH1* hMCNC = mcspecs[1].ToTH1(pot);
106  TH1* hMCNuMu = mcspecs[2].ToTH1(pot);
107  TH1* hMCNuE = mcspecs[3].ToTH1(pot);
108 
109  std::string xLabel(hMCNC->GetXaxis()->GetTitle());
110  if(xLabel.find("PT/P") != std::string::npos) {
111  xLabel.replace(xLabel.find("PT/P"), 4, "p_{T}/p");
112  }
113  if(xLabel.find("RecoE") != std::string::npos) {
114  xLabel = "Visible Energy (GeV)";
115  }
116  if(xLabel.find("RecoELong") != std::string::npos) {
117  xLabel = "Visible Energy (GeV)";
118  }
119 
120  std::string yLabel = "";
121  if(plotlabel.find("RecoE") != std::string::npos) { yLabel = "10^{3} Events / 0.1 GeV / 8.04 #times10^{20} POT"; }
122  if(plotlabel.find("RecoELong") != std::string::npos) { yLabel = "10^{3} Events / 0.1 GeV / 8.04 #times10^{20} POT"; }
123  else yLabel = "10^{3} Events / 8.04 #times10^{20} POT";
124  std::string histTitle = ";" + xLabel + ";" + yLabel;
125 
126  // Numu Stack = Numu
127  TH1* hNumuStack = (TH1*)hMCNuMu->Clone();
128 
129  // Nue Stack = Numu Stack + Nue = Numu + Nue
130  TH1* hNueStack = (TH1*)hNumuStack->Clone();
131  hNueStack->Add(hMCNuE);
132 
133  // NC Stack = Nue Stack + NC = Numu + Nue + NC
134  TH1* hNCStack = (TH1*)hNueStack->Clone();
135  hNCStack->Add(hMCNC);
136 
137  double maxi = 0;
138  double maxDat = hDat ->GetBinContent(hDat ->GetMaximumBin());
139  double maxMCTot = hMCTot ->GetBinContent(hMCTot ->GetMaximumBin());
140  double maxMCNC = hMCNC ->GetBinContent(hMCNC ->GetMaximumBin());
141  double maxMCNuE = hMCNuE ->GetBinContent(hMCNuE ->GetMaximumBin());
142  double maxMCNuMu = hMCNuMu->GetBinContent(hMCNuMu->GetMaximumBin());
143 
144  maxi = std::max(maxDat, maxMCTot);
145  maxi = std::max(maxi, maxMCNC);
146  maxi = std::max(maxi, maxMCNuE);
147  maxi = std::max(maxi, maxMCNuMu);
148 
149  int nbins = hDat->GetNbinsX();
150  for(int i = 1; i <= nbins; ++i) {
151  hDat->SetBinError(i, sqrt(hDat->GetBinContent(i)));
152  }
153 
154  // Scale down the ND by 10e3 (don't forget the max values!)
155  double ndScale = 1./1000.;
156  hDat ->Scale(ndScale);
157  hMCTot ->Scale(ndScale);
158  hMCNC ->Scale(ndScale);
159  hMCNuE ->Scale(ndScale);
160  hMCNuMu->Scale(ndScale);
161  hNCStack ->Scale(ndScale);
162  hNueStack ->Scale(ndScale);
163  hNumuStack->Scale(ndScale);
164  maxi *= ndScale;
165 
166  SetHistOptions(hDat, maxi, histTitle, 0, kBlack, false);
167  SetHistOptions(hMCTot, maxi, histTitle, 0, kTotalMCColor, false);
168  SetHistOptions(hMCNC, maxi, histTitle, 0, kNCBackgroundColor, false);
169  SetHistOptions(hMCNuE, maxi, histTitle, 0, kBeamNueBackgroundColor, false);
170  SetHistOptions(hMCNuMu, maxi, histTitle, 0, kNumuBackgroundColor, false);
171 
172  hDat->SetMarkerStyle(kFullCircle);
173 
174  SetHistOptions(hNueStack, maxi, histTitle, 0, kBeamNueBackgroundColor, true);
175  SetHistOptions(hNumuStack, maxi, histTitle, 0, kNumuBackgroundColor, true);
176  SetHistOptions(hNCStack, maxi, histTitle, 0, kNCBackgroundColor, true);
177 
178  // Create a legend for the histograms
179  double xL = 0.5, xR = 0.7, yB = 0.64, yT = 0.84;
180  TLegend* leg = new TLegend(xL, yB, xR, yT);
181  SetLegendOptions(leg);
182  leg->AddEntry(hDat, "ND Data", "lep");
183  leg->AddEntry(hMCTot, "Total Prediction", "l");
184  leg->AddEntry(hMCNC, "NC Prediction", "l");
185  leg->AddEntry(hMCNuE, "#nu_{e} CC Background", "l");
186  leg->AddEntry(hMCNuMu, "#nu_{#mu} CC Background", "l");
187  leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);
188 
189  TH1* hRat = (TH1*)hDat->Clone("hRat");
190  hRat->Divide(hMCTot);
191  TH1* hLine = (TH1*)hMCTot->Clone("");
192 
193  double maxvalRat = 1., minvalRat = 1.;
194  for(int i = 1; i <= nbins; ++i) {
195  if(hDat->GetBinContent(i) > 0) {
196  double errpct = hDat->GetBinError(i)/hDat->GetBinContent(i);
197  hRat->SetBinError(i, errpct);
198  }
199 
200  double binval = hRat->GetBinContent(i);
201  double errval = hRat->GetBinError(i);
202  if(binval > 0) {
203  maxvalRat = std::max(maxvalRat, binval + errval);
204  minvalRat = std::min(minvalRat, binval - errval);
205  }
206 
207  hLine->SetBinContent(i, 1.);
208  }
209 
210  hRat->GetYaxis()->SetTitle("Data / Total MC");
211  hRat->SetLineColor(kBlack);
212  hLine->GetYaxis()->SetTitle("Data / Total MC");
213  hLine->SetLineColor(kTotalMCColor);
214 
215  hRat ->SetMaximum(1.05*maxvalRat);
216  hLine->SetMaximum(1.05*maxvalRat);
217  hRat ->SetMinimum(0.95*minvalRat);
218  hLine->SetMinimum(0.95*minvalRat);
219 
220  CenterTitles(hRat);
221 
222  if(plotlabel.find("NM1") == std::string::npos &&
223  plotlabel.find("ECalLongAllE") == std::string::npos) {
224  TCanvas* cnorat = new TCanvas(plotlabel.c_str(), plotlabel.c_str(), 800, 500);
225  gPad->SetFillStyle(0);
226  hDat->SetMinimum(0.001);
227  hMCTot->SetMinimum(0.001);
228  hMCNC->SetMinimum(0.001);
229  hMCNuE->SetMinimum(0.001);
230  hMCNuMu->SetMinimum(0.001);
231  //gPad->SetLogy();
232  hDat->Draw("E1");
233  hMCTot->Draw("hist same");
234  hMCNC->Draw("hist same");
235  hMCNuE->Draw("hist same");
236  hMCNuMu->Draw("hist same");
237  leg->Draw();
238  gPad->RedrawAxis();
239 
240  Preliminary();
241  plotsOut->WriteTObject(cnorat);
242 
243  // Create a canvas for and draw the stacked histograms
244  TCanvas* cStack = new TCanvas((plotlabel + "Stack").c_str(),
245  (plotlabel + "Stack").c_str(), 800, 500);
246  gPad->SetFillStyle(0);
247  //gPad->SetLogy();
248  TLegend* legStack = new TLegend(xL, yB, xR, yT);
249  SetLegendOptions(legStack);
250  legStack->AddEntry(hDat, "ND Data", "lep");
251  legStack->AddEntry(hNCStack, "NC 3 Flavor Prediction", "f");
252  legStack->AddEntry(hNueStack, "#nu_{e} CC Background", "f");
253  legStack->AddEntry(hNumuStack, "#nu_{#mu} CC Background", "f");
254  legStack->SetY1(legStack->GetY2() - legStack->GetNRows()*0.05);
255 
256  hNCStack ->Draw("hist");
257  hNueStack ->Draw("hist same");
258  hNumuStack->Draw("hist same");
259  hDat->Draw("same");
260  legStack->Draw();
261  gPad->RedrawAxis();
262 
263  if(plotlabel.find("RecoELong") != std::string::npos) {
264  FILE* textFO;
265  textFO = fopen("out2.txt", "a+");
266  // Spectrum sData = hDat.FakeData(pot);
267  fprintf(textFO, "Event numbers at the %2.2s for %s %.2fe20 POT:\n", "ND", plotlabel.c_str(), 8.04);
268  if(hDat) { fprintf(textFO, "%12.12s, ", "Data"); }
269  else { fprintf(textFO, "%14.14s", ""); }
270  fprintf(textFO, "%12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
271  "All", "NC (3 Flav)", "Numu", "Nue", "FOM");
272 
273  double nAll = hMCTot ->Integral();
274  double nNC3F = hMCNC ->Integral();
275  double nNumu = hMCNuMu->Integral();
276  double nNue = hMCNuE ->Integral();
277  double fom = nNC3F/sqrt(nAll);
278  if(hDat) {
279  double nData = hDat->Integral(1, hDat->GetNbinsX());
280  fprintf(textFO, "%10E, ", nData);
281  }
282  else {
283  fprintf(textFO, "%14.14s", "");
284  }
285  fprintf(textFO, "%10E, %10E, %10E, %10E, %10E\n\n",
286  nAll, nNC3F, nNumu, nNue, fom);
287  fclose(textFO);
288  }
289  Preliminary();
290  plotsOut->WriteTObject(cStack);
291  }
292  else {
293  FILE* textF;
294  textF = fopen("out.txt", "a+");
295  Spectrum sData = data.FakeData(pot);
296  std::cout << "POT: " << pot << std::endl;
297  /// TODO: The pot values need to not be hard-coded here
298  PlotStack(mcspecs, plotsOut, textF, plotlabel, "ND",
299  8.04e20, 8.04, GetPlotOptions(plotlabel), &sData);
300  fclose(textF);
301  }
302 
303  TCanvas* canvas = new TCanvas((plotlabel + "Rat").c_str(),
304  (plotlabel + "Rat").c_str(), 800, 800);
305  gPad->SetFillStyle(0);
306  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
307  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
308  pSpecs->Draw();
309  pRatio->Draw();
310 
311  pSpecs->cd();
312  hDat->Draw("E1");
313  hDat->GetXaxis()->SetRangeUser(0,5);
314  hDat->GetXaxis()->SetNdivisions(-405);
315  hMCTot->Draw("hist same");
316  hMCNC->Draw("hist same");
317  hMCNuE->Draw("hist same");
318  hMCNuMu->Draw("hist same");
319  leg->Draw();
320  gPad->RedrawAxis();
321 
322  Preliminary();
323  pRatio->cd();
324  hLine->Draw("hist ][");
325  hRat->Draw("E1 same");
326 
327  plotsOut->WriteTObject(canvas);
328  std::cout << "done\n";
329 };
330 
331 //------------------------------------------------------------------------------
333 {
335 
336  int div = -404;
337  double xL = 0.5, xR = 0.7, yB = 0.62, yT = 0.87;
338  double texX = xL + .01;
339  double texY = yB - 0.045;
340  double cutL = 1.;
341  double cutR = 0.;
342 
343  if(name.find("RecoELong") != std::string::npos) {
344  div = -406;
345  xL = 0.35;
346  xR = 0.55;
347  texX = xL + .01;
348  cutL = 0.5;
349  cutR = 4.0;
350  }
351  if(name.find("VtxX") != std::string::npos) {
352  div = -404;
353  xL = 0.4;
354  xR = 0.6;
355  texY = yT - 0.045;
356  cutL = -100.;
357  cutR = 100.;
358  }
359  if(name.find("VtxY") != std::string::npos) {
360  div = -404;
361  xL = 0.4;
362  xR = 0.6;
363  texX = 0.2;
364  texY = yT - 0.045;
365  cutL = -100.;
366  cutR = 100.;
367  }
368  if(name.find("VtxZ") != std::string::npos) {
369  div = -503;
370  xL = 0.28;
371  xR = 0.48;
372  texX = xL + 0.01;
373  cutL = 200.;
374  cutR = 1000.;
375  }
376  if(name.find("CVN") != std::string::npos) {
377  div = -505;
378  xL = 0.35;
379  xR = 0.55;
380  texX = xL + 0.01;
381  cutL = 0.2;
382  }
383  if(name.find("NHit") != std::string::npos) {
384  div = -505;
385  cutL = 20.;
386  cutR = 200.;
387  }
388  if(name.find("PTP") != std::string::npos) {
389  div = -505;
390  xL = 0.125;
391  xR = 0.325;
392  texX = 0.5;
393  texY = yT - 0.045;
394  cutR = 0.8;
395  }
396  if(name.find("EpH") != std::string::npos) {
397  div = -505;
398  cutL = 0.009;
399  }
400  if(name.find("DistTop") != std::string::npos) {
401  div = -404;
402  xL = 0.4;
403  xR = 0.6;
404  texX = 0.675;
405  texY = yT - 0.045;
406  }
407 
408  ret.ndiv = div;
409  ret.xL = xL;
410  ret.xR = xR;
411  ret.yT = yT;
412  ret.yB = yB;
413  ret.texX = texX;
414  ret.texY = texY;
415  ret.cutL = cutL;
416  ret.cutR = cutR;
417 
418  return ret;
419 }
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string det, double POT, double potEquiv, PlotOptions opt, Spectrum *dataspec)
void SetHistOptions(TH1 *h, double max, std::string title, int ndiv, Color_t col, bool fill)
Set common options for a TLegend.
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Spectrum NCTotalComponent() const override
Definition: CheatDecomp.h:27
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
#define MAXSPEC
Definition: NusLoadProd3.h:18
T sqrt(T number)
Definition: d0nt_math.hpp:156
void Clear()
Definition: Spectrum.cxx:361
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
const Color_t kNumuBackgroundColor
Definition: Style.h:30
void make_plots(const ana::Spectrum &data, const ana::CheatDecomp &MC, std::string plotlabel, TFile *plotsOut)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char const XML_Char * data
Definition: expat.h:268
Double_t maxi
Definition: plot.C:28
static std::unique_ptr< CheatDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Definition: CheatDecomp.cxx:73
Spectrum NueComponent() const override
Definition: CheatDecomp.h:32
fclose(fg1)
const int nbins
Definition: cellShifts.C:15
Spectrum AntiNumuComponent() const override
Definition: CheatDecomp.h:33
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
#define pot
Spectrum NumuComponent() const override
Definition: CheatDecomp.h:31
double POT() const
Definition: Spectrum.h:227
void Preliminary()
OStream cout
Definition: OStream.cxx:6
void make_dataMC(TFile *SavedFile, TFile *OutFile)
void SetLegendOptions(TLegend *leg)
Set common options for a TLegend.
void DataMCNDAna_nus17(std::string folder)
const Color_t kTotalMCColor
Definition: Style.h:16
TFile * OutFile
T min(const caf::Proxy< T > &a, T b)
const Color_t kNCBackgroundColor
Definition: Style.h:22
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
PlotOptions GetPlotOptions(std::string name)
void next()
Definition: show_event.C:84
Spectrum AntiNueComponent() const override
Definition: CheatDecomp.h:34
enum BeamMode string