hadEFrac_nd_data_mc_systs.C
Go to the documentation of this file.
1 // Fills spectra for systematic error band. Makes both pot norm with full syst band and area norm with shape only syst band
2 // have to manually make nd_syst_plots directory before running
3 
4 #include "CAFAna/Cuts/Cuts.h"
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Core/Var.h"
13 #include "CAFAna/Systs/Systs.h"
19 #include "CAFAna/Analysis/Plots.h"
20 #include "CAFAna/Analysis/Style.h"
21 //#include "Utilities/rootlogon.C"
22 
23 #include "TCanvas.h"
24 #include "TFile.h"
25 #include "TGraph.h"
26 #include "TH1.h"
27 #include "TH2.h"
28 #include "TMath.h"
29 #include "TGaxis.h"
30 #include "TMultiGraph.h"
31 #include "TLegend.h"
32 #include "TLegendEntry.h"
33 #include "TLatex.h"
34 #include "TStyle.h"
35 #include "THStack.h"
36 #include "TPaveText.h"
37 #include "TList.h"
38 #include "TGaxis.h"
39 #include "TAttLine.h"
40 #include "TAttMarker.h"
41 
42 #include <cmath>
43 #include <iostream>
44 #include <vector>
45 #include <list>
46 #include <sstream>
47 #include <string>
48 #include <sstream>
49 #include <fstream>
50 #include <iomanip>
51 
52 
53 
54 using namespace ana;
55 
56 
57 // weight bins by bin width for neutrino energy plots:
58 const Var kBinWidthWeight([](const caf::SRProxy *sr){
59  const double Enu = kCCE(sr); // If want bin normalised
60  //const double Enu = 6; // If don't want bin normalised.
61  if(Enu < 0.75) return 0.1 / 0.75;
62  if(Enu >= 0.75 && Enu < 1) return 0.1 / 0.25;
63  if(Enu >= 1 && Enu < 2) return 0.1 / 0.1;
64  if(Enu >= 2 && Enu < 3) return 0.1 / 0.25;
65  if(Enu >= 3 && Enu < 4) return 0.1 / 0.5;
66  if(Enu >= 4 && Enu < 5) return 0.1 / 1;
67  else return 1.;
68  });
69 
70 
71 // Put a "NOvA Preliminary" tag in the corner
73 {
74  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
75  prelim->SetTextColor(kBlue);
76  prelim->SetNDC();
77  prelim->SetTextSize(2/30.);
78  prelim->SetTextAlign(32);
79  prelim->Draw();
80 }
81 /*
82 void CenterTitles(TH1* histo)
83 {
84  histo->GetXaxis()->CenterTitle();
85  histo->GetYaxis()->CenterTitle();
86  histo->GetZaxis()->CenterTitle();
87  histo->GetYaxis()->SetDecimals();
88  histo->GetXaxis()->SetNoExponent();
89 }
90 */
92 {
93  for(const auto& obj:*(gPad->GetListOfPrimitives()))
94  {
95 
96  if(obj->InheritsFrom("TH1")) CenterTitles((TH1*)obj);
97  }
98 }
99 
100 
101 
103 {
104  // --- Loop through the objects in the pad.
105  for(const auto& obj:*(gPad->GetListOfPrimitives())) {
106  // --- Want to repaint all TH1's.
107  if(obj->InheritsFrom("TH1")) {
108  // --- Check if I have the Reco Energy canvas?
109  std::string XTit = ((TH1*)obj)->GetXaxis()->GetTitle();
110  if ( XTit.find("Reconstructed Neutrino Energ") != std::string::npos ) {
111  ((TH1*)obj)->GetYaxis()->SetTitle("Events / 0.1 GeV");
112  }
113  // Repaint the axes.
114  ((TH1*)obj)->GetXaxis()->Paint();
115  //((TH1*)obj)->GetYaxis()->SetTitleOffset(1.1);
116  ((TH1*)obj)->GetYaxis()->Paint();
117  ((TH1*)obj)->GetZaxis()->Paint();
118  }
119  }
120 }
121 
122 
123 TCanvas* LogClone(const TCanvas* can)
124 {
125 
126  TCanvas* clone = (TCanvas*)can->DrawClone();
127  clone->cd();
128 
129  clone->SetCanvasSize(can->GetWw(), can->GetWh());
130 
132  float max = 0;
133 
134  for(const auto& obj:*(clone->GetListOfPrimitives()))
135  {
136  if(obj->InheritsFrom("TH1"))
137  {
138  TH1* h = (TH1*)obj;
139  for(int i = 1; i <= h->GetNbinsX(); ++i)
140  {
141  int c = h->GetBinContent(i);
142  if(c != 0 && c < min) min = c;
143  if (c > max) max = c;
144  }
145  }
146  }
147  min *= 0.8;
148  float span = log10(max) - log10(min);
149  span *= 1.5;
150  max = pow(10, span + log10(min));
151  for(const auto& obj:*(clone->GetListOfPrimitives()))
152  {
153  if(obj->InheritsFrom("TH1"))
154  {
155  ((TH1*)obj)->SetMinimum(min);
156  ((TH1*)obj)->SetMaximum(max);
157  }
158  }
159 
160 
161 
162  TString name = TString(can->GetTitle()) + "_logy";
163  clone->SetName(name);
164  clone->SetTitle(name);
165  clone->SetLogy();
166  return clone;
167 }
168 
169 template<class T> class Tangible
170 {
171 
172 public:
173  Tangible(const T& obj, const std::string& shortName,
174  const std::string& blurb ):
175  fObj(obj),
176  fShortName(shortName),
177  fBlurb(blurb)
178  {};
179 
180  T fObj;
183 
184 };
185 
186 
189 
190 
191 
193 {
194 public:
195 
196  DataMCPair(Selection sel, TangibleAxis tanAxis,
197  SpectrumLoader& loaderData, SpectrumLoader& loaderMC,
198  std::vector<const ISyst*> systs, const Var BinWeight, const Cut& bkg=kNoCut ):
199  fSel(sel),
200  fAxis(tanAxis),
201  //fData(loaderData, fAxis.fObj, fSel.fObj, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight),
202  //fMC(loaderMC, fAxis.fObj, fSel.fObj, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight),
203  //fMCBkg(loaderMC, fAxis.fObj, fSel.fObj && bkg, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight),
204  fData (loaderData, fAxis.fObj, fSel.fObj , kNoShift, BinWeight),
205  fMC (loaderMC , fAxis.fObj, fSel.fObj , kNoShift, BinWeight),
206  fMCBkg(loaderMC , fAxis.fObj, fSel.fObj && bkg, kNoShift, BinWeight),
207  fShortName(fAxis.fShortName + "_" + fSel.fShortName)
208  {
209  fUps.reserve(systs.size());
210  fDowns.reserve(systs.size());
211  for(const auto& syst:systs){
212  //fUps.emplace_back(loaderMC, fAxis.fObj, fSel.fObj, SystShifts(syst, +1), kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight);
213  //fDowns.emplace_back(loaderMC, fAxis.fObj, fSel.fObj, SystShifts(syst, -1), kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight);
214  fUps .emplace_back(loaderMC, fAxis.fObj, fSel.fObj, SystShifts(syst, +1), BinWeight);
215  fDowns.emplace_back(loaderMC, fAxis.fObj, fSel.fObj, SystShifts(syst, -1), BinWeight);
216  }
217  };
218 
220  {return fShortName;};
221  const char* CName() const
222  {return ShortName().c_str();};
223 
224 
225 
226 
227  void DrawMCSyst(const int iSyst = -1) const
228  {
229  if(iSyst < 0)
230  PlotWithSystErrorBand(fMC, fUps, fDowns, fData.POT(),
232  else
233  {
234  assert(iSyst < fUps.size());
235  PlotWithSystErrorBand(fMC, {fUps[iSyst]}, {fDowns[iSyst]}, fData.POT(),
237  }
238 
239  }
240  void OverlayDataMCSyst(const int iSyst = -1) const
241  {
242  DrawData();
243  DrawMCSyst(iSyst);
244  DrawData();
245  DrawMCBkg();
246  TLegend* leg = DrawLegendPOT();
247  AestheticsPOT();
248 
249  };
250 
251  void DrawMCNormSyst(const int iSyst = -1) const
252  {
253  // New vectors for area normalized systs
254  std::vector<Spectrum> normUps;
255  std::vector<Spectrum> normDowns;
256 
257  Spectrum tempMC(fMC.ToTH1(fMC.POT()),fMC.POT(),fMC.Livetime());
258  tempMC.OverridePOT(fMC.POT() * fMC.ToTH1(fData.POT())->GetSumOfWeights() / fData.ToTH1(fData.POT())->GetSumOfWeights());
259 
260  if (iSyst >= 0)
261  {
262  // We'll have one entry if there's just one syst to plot.
263  normUps.reserve(1);
264  normDowns.reserve(1);
265  }
266  else
267  {
268  // Otherwise we need all systs.
269  normUps.reserve(fUps.size());
270  normDowns.reserve(fUps.size());
271  }
272  for(size_t i = 0; i < fUps.size(); ++i)
273  {
274  // bail if we want a particular syst and it's not the one we want.
275  if(iSyst >= 0 && (int)i != iSyst) continue;
276  normUps.emplace_back(fUps[i]);
277  Spectrum& norm = normUps.back();
278  norm.OverridePOT(norm.POT() * norm.Integral(1) / tempMC.Integral(1));
279  }
280  for(size_t i = 0; i < fDowns.size(); ++i)
281  {
282  // bail if we want a particular syst and it's not the one we want.
283  if(iSyst >= 0 && (int)i != iSyst) continue;
284  normDowns.emplace_back(fDowns[i]);
285  Spectrum& norm = normDowns.back();
286  norm.OverridePOT(norm.POT() * norm.Integral(1) / tempMC.Integral(1));
287  }
288  PlotWithSystErrorBand(tempMC, normUps, normDowns, fData.POT(),
290 
291  }
292 
293  //// Weird thing, normalize the systs to area of the MC
294  void OverlayDataMCSystNorm(const int iSyst = -1) const
295  {
296 
297  const int tempcolor=kBlack;
298  const int tempmarker=kFullCircle;
299  TH1D* hDummy = fData.ToTH1(fData.POT());
300  hDummy->GetYaxis()->SetTitle("Events / 0.1 GeV"); // only for cust bin neutrino plot
301  //hDummy->GetYaxis()->SetTitleOffset(1.1);
302  hDummy->SetMarkerColor(tempcolor);
303  hDummy->SetMarkerStyle(tempmarker);
304  hDummy->SetLineColor(tempcolor);
305  hDummy->Draw();
306 
307  DrawMCNormSyst(iSyst);
308  DrawData();
309  DrawMCBkg();
310  TLegend* leg = DrawLegendArea();
311  AestheticsArea();
312  };
313 
314 
315  //// Weird thing, normalize the systs to area of the MC
316  void OverlayDataMCSystExtraData(const DataMCPair& extraData,
317  const std::string dataLegendTitle,
318  const std::string extraLegendTitle,
319  const int iSyst = -1,
320  const int dataColor=kBlack,
321  const int extraColor=kGray+1,
322  const int dataMarker=kFullCircle,
323  const int extraMarker=kFullCircle) const
324  {
325 
326  DrawMCSyst(iSyst);
327  DrawData(dataColor, dataMarker);
328  extraData.DrawData(extraColor, extraMarker);
329  DrawMCBkg();
330  TLegend* leg = DrawLegendPOT(dataColor, dataLegendTitle, dataMarker);
331  TH1F* colExtraData = new TH1F();
332  colExtraData->SetMarkerColor(extraColor);
333  colExtraData->SetMarkerStyle(extraMarker);
334  colExtraData->SetLineColor(extraColor);
335  leg->AddEntry((TObject*)colExtraData, extraLegendTitle.c_str(), "ple");
336  AestheticsPOT();
337 
338  };
339 
340 
341 
342  //// Weird thing, normalize the systs to area of the MC
344  const std::string dataLegendTitle,
345  const std::string extraLegendTitle,
346  const int iSyst = -1,
347  const int dataColor=kBlack,
348  const int extraColor=kGray+1,
349  const int dataMarker=kFullCircle,
350  const int extraMarker=kFullCircle) const
351  {
352 
353 
354  DrawMCNormSyst(iSyst);
355  DrawData(dataColor, dataMarker);
356  extraData.DrawData(extraColor);
357  DrawMCBkg();
358  TLegend* leg = DrawLegendArea(dataColor, dataLegendTitle, dataMarker);
359  TH1F* colExtraData = new TH1F();
360  colExtraData->SetMarkerColor(extraColor);
361  colExtraData->SetMarkerStyle(extraMarker);
362  colExtraData->SetLineColor(extraColor);
363  leg->AddEntry((TObject*)colExtraData, extraLegendTitle.c_str(), "ple");
364  AestheticsArea();
365 
366  };
367 
368 
369  TLegend* DrawLegendArea(const int dataColor=kBlack,
370  std::string dataLegendTitle="Data",
371  const int dataMarker=kFullCircle) const
372  {
373 
374  TLegend* leg = AutoPlaceLegend(0.36, 0.19, 0.68);
375 
376  TH1F* colMC = new TH1F();
377  colMC->SetLineColor(kTotalMCColor);
378  TH1F* colMCBkg = new TH1F();
379  colMCBkg->SetLineColor(kBlue);
380  TH1F* colData = new TH1F();
381  colData->SetLineColor(dataColor);
382  colData->SetMarkerColor(dataColor);
383  colData->SetMarkerStyle(dataMarker);
384 
385  leg->AddEntry((TObject*)colMC, "Simulated Selected Events", "l");
386  leg->AddEntry((TObject*)colMCBkg, "Simulated Background", "l");
387  leg->AddEntry((TObject*)colData, dataLegendTitle.c_str(), "ple");
388  TLegendEntry *entryns1=leg->AddEntry("error","Shape-only 1-#sigma syst. range","f");
389  entryns1->SetFillStyle(1001);
390  entryns1->SetFillColor(kRed-10);
391  leg->SetTextSize(0.04);
392 
393  leg->SetFillColor(0);
394  leg->SetFillStyle(0);
395 
396  for(const auto& obj:*leg->GetListOfPrimitives())
397  {
398  if(obj->InheritsFrom("TAttFill")){
399  ((TAttFill*)obj)->SetFillStyle(0);
400  //((TAttFill*)obj)->SetFillColor(0);
401  }
402 
403  }
404  leg->Draw();
405  return leg;
406  }
407 
408  TLegend* DrawLegendPOT(const int dataColor=kBlack,
409  std::string dataLegendTitle="Data",
410  const int dataMarker=kFullCircle) const
411  {
412 
413  TLegend* leg = AutoPlaceLegend(0.36, 0.19, 0.68);
414 
415  TH1F* colMC = new TH1F();
416  colMC->SetLineColor(kTotalMCColor);
417  TH1F* colMCBkg = new TH1F();
418  colMCBkg->SetLineColor(kBlue);
419  TH1F* colData = new TH1F();
420  colData->SetLineColor(dataColor);
421  colData->SetMarkerColor(dataColor);
422  colData->SetMarkerStyle(dataMarker);
423 
424  leg->AddEntry((TObject*)colMC, "Simulated Selected Events", "l");
425  leg->AddEntry((TObject*)colMCBkg, "Simulated Background", "l");
426  leg->AddEntry((TObject*)colData, dataLegendTitle.c_str(), "ple");
427  TLegendEntry *entryns1=leg->AddEntry("error","Full 1-#sigma syst. range","f");
428  entryns1->SetFillStyle(1001);
429  entryns1->SetFillColor(kRed-10);
430  leg->SetTextSize(0.04);
431 
432  leg->SetFillColor(0);
433  leg->SetFillStyle(0);
434 
435  for(const auto& obj:*leg->GetListOfPrimitives())
436  {
437  if(obj->InheritsFrom("TAttFill")){
438  ((TAttFill*)obj)->SetFillStyle(0);
439  //((TAttFill*)obj)->SetFillColor(0);
440  }
441 
442  }
443  leg->Draw();
444  return leg;
445  }
446 
447  void AddExposureArea() const
448  {
449  std::stringstream expo;
450  std::stringstream expoData;
451  std::stringstream expoMC;
452  float pot = fData.POT();
453  float dataMean = (fData.ToTH1(pot)) -> GetMean();
454  float dataMeanError = (fData.ToTH1(pot)) -> GetMeanError();
455  float MCMean = (fMC.ToTH1(pot)) -> GetMean();
456  float MCMeanError = (fMC.ToTH1(pot)) -> GetMeanError();
457  expo.precision(2);
458  expoData.precision(2);
459  expoMC.precision(2);
460  expo << "ND area norm., "
461  << TString::Format("%.02lf #times 10^{20}", pot/1e20) << " POT " <<endl;
462  //expoData << TString::Format("Data mean: %.02lf #pm %.02lf GeV ", dataMean, dataMeanError) <<endl;
463  //expoMC << TString::Format("MC mean: %.02lf #pm %.02lf GeV ", MCMean, MCMeanError) <<endl;
464  expoData << TString::Format("Data mean: %.02lf GeV ", dataMean) <<endl;
465  expoMC << TString::Format("MC mean: %.02lf GeV ", MCMean ) <<endl;
466 
467  TLatex text;
468  TLegend* place = AutoPlaceLegend(0.0001, 0.0001, 0.52);
469  text.SetTextAlign(22);
470  text.SetTextSize(0.04);
471 
472  // text.DrawLatexNDC(place->GetX1(), place->GetY1(),
473  // (expo.str()).c_str() );
474 
475  text.DrawLatexNDC(place->GetX1(), place->GetY1(),
476  ("#splitline{" + expo.str() + "}{ #splitline{" +
477  expoData.str() + "}{" + expoMC.str() + "}}").c_str()
478  );
479 
480  }
481 
482  void AddExposurePOT() const
483  {
484  std::stringstream expo;
485  float pot = fData.POT();
486  expo.precision(2);
487  expo << "ND POT norm., "
488  << TString::Format("%.02lf #times 10^{20}", pot/1e20) << " POT";
489  TLatex text;
490  TLegend* place = AutoPlaceLegend(0.0001, 0.0001, 0.56);
491  text.SetTextAlign(22);
492  text.SetTextSize(0.04);
493 
494  text.DrawLatexNDC(place->GetX1(), place->GetY1(), expo.str().c_str());
495 
496  if (ShortName().find("numuE") != std::string::npos) {
497  std::cout << "\n" << ShortName() << "\t Data: " << fData.ToTH1(pot)->Integral()
498  << "\t MC: " << fMC.ToTH1(pot)->Integral() << "\n";
499  }
500  }
501 
502  void AestheticsArea() const
503  {
504  CenterTitles();
505  RedrawAxes();
506  AddExposureArea();
507  }
508 
509  void AestheticsPOT() const
510  {
511  CenterTitles();
512  RedrawAxes();
513  AddExposurePOT();
514  }
515 
516  void DrawData(const int color=kBlack, const int marker=kFullCircle)const
517  {
518  TH1D* hData = fData.ToTH1(fData.POT());
519  hData->GetYaxis()->SetTitle("Events / 0.1 GeV"); // only for cust bin neutrino plot
520  hData->SetMarkerColor(color);
521  hData->SetMarkerStyle(marker);
522  hData->SetLineColor(color);
523 
524  hData->Draw("SAME");
525  };
526 
527  void DrawMCBkg()const
528  {
529  TH1D* hMCBkg = fMCBkg.ToTH1(fData.POT());
530  hMCBkg->GetYaxis()->SetTitle("Events / 0.1 GeV"); // only for cust bin neutrino plot
531  hMCBkg->SetLineColor(kBlue);
532  hMCBkg->Draw("HIST SAME");
533  };
534 
535  float Purity()const
536  {
537  return 1 - fMCBkg.ToTH1(fData.POT())->GetEntries() /
538  fMC.ToTH1(fData.POT())->GetEntries();
539  }
540 
541  Selection fSel;
542  TangibleAxis fAxis;
546  std::vector<Spectrum> fUps;
547  std::vector<Spectrum> fDowns;
549 
550 
551 };
552 
553 
554 //----------------------------------------------------------------------
555  /// Get a vector of all the numu group systs
556 /*std::vector<const ISyst*> getAllNumuSystsLimited()
557 {
558  std::vector<const ISyst*> systs =
559  {&kNumuNormSyst,
560  &kNumuRelNormSyst,
561  &kNumuNCScaleSyst,
562  &kNumuTauContaminationSyst,
563  //&kNumuGEANTScaleSyst,
564  //&kNumuGEANTNormSyst,
565  &kFitHadEnergyScaleSyst,
566  &kFDFitHadEnergyScaleSyst,
567  &kMECScaleSyst,
568  &kRPACCQESyst,
569  &kFitMuEnergyScaleSyst,
570  &kFDFitMuEnergyScaleSyst
571  };
572 */
573 
574 
575  //if(includeBeam) systs.push_back(&kBeamAll);
576 
577  /*
578  systs.push_back(&kNumuSummedSmallGENIESyst);
579  std::set<rwgt::ReweightLabel_t> largeGENIEParameters;
580  largeGENIEParameters.insert(rwgt::fReweightMaNCEL);
581  // largeGENIEParameters.insert(rwgt::fReweightMaCCQE);
582  largeGENIEParameters.insert(rwgt::fReweightMaCCQEshape);
583  largeGENIEParameters.insert(rwgt::fReweightNormCCQE);
584  largeGENIEParameters.insert(rwgt::fReweightMaCCRES);
585  largeGENIEParameters.insert(rwgt::fReweightMvCCRES);
586  largeGENIEParameters.insert(rwgt::fReweightMaNCRES);
587  largeGENIEParameters.insert(rwgt::fReweightMvNCRES);
588  for (const auto & parameter:largeGENIEParameters)
589  systs.push_back(new GenieRwgtTableSyst(parameter));
590  */
591 
592 // return systs;
593 //}
594 
595 
596 // std::vector<const ISyst*> getSyst_AhtBY()
597 // {
598 // std::vector<const ISyst*> systs =
599 // {&kDIS_AhtBY_Syst};
600 // return systs;
601 // }
602 // std::vector<const ISyst*> getSyst_BhtBY()
603 // {
604 // std::vector<const ISyst*> systs =
605 // {&kDIS_BhtBY_Syst};
606 // return systs;
607 // }
608 // std::vector<const ISyst*> getSyst_CV1uBY()
609 // {
610 // std::vector<const ISyst*> systs =
611 // {&kDIS_CV1uBY_Syst};
612 // return systs;
613 // }
614 // std::vector<const ISyst*> getSyst_CV2uBY()
615 // {
616 // std::vector<const ISyst*> systs =
617 // {&kDIS_CV2uBY_Syst};
618 // return systs;
619 // }
620 
621 
622 /// Return an iterator if all keys are found in only one entry
623 //std::vector<DataMCPair>::iterator findPair(std::vector<DataMCPair>& list,
624  DataMCPair* findPair(std::vector<DataMCPair>& list,
625  std::vector<std::string> keys)
626 {
627  int nFound = 0;
628 
629  std::vector<DataMCPair>::iterator it = list.end();
630 
631  for(std::vector<DataMCPair>::iterator entry = list.begin();
632  entry != list.end();
633  ++entry)
634  {
635  /// Assume we're going to find it, innocent until guilty
636  bool found = true;
637  for(const auto& key:keys)
638  {
639  std::cout << entry->ShortName() << " " << key << " " << found << std::endl;
640  std::cout << (entry->ShortName().find(key) == std::string::npos) << " " << entry->ShortName().find(key) << " " << std::string::npos << std::endl;
641  // If any key is not present, this isn't the one.
642  if(entry->ShortName().find(key) == std::string::npos)
643  {
644 
645  found = false;
646  std::cout << entry->ShortName() << " " << key << " " << found << " " <<entry->ShortName().find(key) << std::endl;
647  break;
648 
649  }
650 
651  }
652  if(found)
653  {
654  it = entry;
655  nFound += 1;
656  }
657 
658  }
659  std::cout << "nFound: "<< nFound << std::endl;
660  if(nFound == 0) throw "Failed to find match.";
661  if(nFound > 1) throw "Found too many matches.";
662 
663  return &(*it);
664 }
665 
666 
668  std::string extra = "", std::string beginning = "")
669 {
670  std::ofstream out(dir + name);
671  out << pair.fAxis.fBlurb << pair.fSel.fBlurb << extra;
672 }
673 
674 ////////////////////////////////////////
675 // Making a text file for blessing:
676 void MakeTextFile( std::string OutName ) {
677  // Determine output name, and caption.
678  std::string FNa = "nd_syst_plots/"+OutName+".txt";
679  // Replace underscores in the canvas name
680  while (OutName.find("_") != std::string::npos) OutName.replace(OutName.find("_"),1," ");
681  // Add a string about it being a cut level.
682  if (OutName.find("Qual") != std::string::npos) OutName.insert(OutName.find(" Qual"), ", Cut level:");
683  // Add a string about it being a cut level.
684  if (OutName.find("pot") != std::string::npos) OutName.replace(OutName.find(" pot"),4,". With POT normalisation.");
685  // Add a string about it being a cut level.
686  if (OutName.find("area") != std::string::npos) OutName.replace(OutName.find(" area"),5,". With Area normalisation.");
687 
688  std::string Cap = "Plot showing the number of ND NuMi events in #nu_{#mu} Ana17 for variable: " + OutName;
689  std::cout << "\nFile name: " << FNa << "\n\tCaption: " << Cap << std::endl;
690  // Write to file.
691  std::ofstream TxtOut ( FNa.c_str(), std::ofstream::out );
692  TxtOut << Cap;
693  TxtOut.close();
694  // Done.
695  return;
696 }
697 
698 ////////////////////////////////////////
699 // Main class function:
700 void hadEFrac_nd_data_mc_systs( bool RunBinNorm = false, bool IsTest=false )
701 {
702 
703  ////////////////////////////////////////
704  // Get stuff for the hadEFrac Cut:
705  ////////////////////////////////////////
706  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used
707 
708  std::string fdspecfile = "/nova/app/users/karlwarb/Workspace/NuMuSystematics/FDQuantileHists/final_full_FD_histo_for_quantile_cuts.root";
709  TFile* inFile = TFile::Open( fdspecfile.c_str() );
710  gDirectory->cd("dir_FDSpec2D");
711  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
712  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
713 
714  ////////////////////////////////////////
715  // Define Axes and cuts.
716  ////////////////////////////////////////
717 
718  std::vector<Cut> TieredCuts;
719  std::vector<std::string> CutNames;
720  //CutNames.emplace_back("Qual"); TieredCuts.emplace_back(kNumuQuality);
721  CutNames.emplace_back("Qual_Cont"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017);
722  CutNames.emplace_back("Qual_Cont_PID"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2017);
723  //CutNames.emplace_back("Qual_Cont_PID_True"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2017 && kIsNumuCC);
724 
725  const HistAxis EnergyAxis ( "Reconstructed Muon Energy (GeV)", Binning::Simple(50, 0, 5 ), kMuE );
726  const HistAxis LengthAxis ( "Length of Primary Track (m)" , Binning::Simple(50, 0, 16 ), kTrkLength );
727  const HistAxis HadEAxis ( "Hadronic Energy (GeV)" , Binning::Simple(50, 0, 3 ), kHadE );
728  const HistAxis HadEFracAxis( "Hadronic Energy Fraction" , Binning::Simple(50, 0, 1.01), kHadEFrac );
729  const HistAxis CosNuMiAxis ( "Kalman Track Cos #theta_{NuMI}" , Binning::Simple(50, 0, 1.01), kCosNumi );
730  const HistAxis ReMIDAxis ( "ReMID score" , Binning::Simple(50, 0, 1.01), kRemID );
731  const HistAxis CVNNuMuAxis ( "CVN NuMu score" , Binning::Simple(50, 0, 1.01), kCVNm );
732 
733  ////////////////////////////////////////
734  // Now make the plots....
735  ////////////////////////////////////////
736  gStyle->SetMarkerStyle(kFullCircle);
737  TGaxis::SetMaxDigits(3);
738 
739  // --- Configure the definition to use.
740  std::string fnameNDMC = "prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
741  std::string fnameNDData = "prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns";
742  if (!IsTest) {
743  //fnameNDMC = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017"; // Takes a VERY long time to run in one macro
744  //fnameNDMC = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017_stride10";
745  fnameNDMC = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017_stride22";
746  //fnameNDMC = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017_lim1";
747  fnameNDData= "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns_numu2017";
748  }
749  // --- Set the Loaders, along with spill cuts.
750  SpectrumLoader loaderMC(fnameNDMC);
751  SpectrumLoader loaderData(fnameNDData);
753  loaderData.SetSpillCut(kStandardSpillCuts);
754 
755  // --- Load the systematics
756  std::vector<const ISyst*> systs = getAllDirectNumuSysts2017();
757  //std::vector<const ISyst*> systs = getSyst_AhtBY();
758  //std::vector<const ISyst*> systs = getSyst_BhtBY();
759  //std::vector<const ISyst*> systs = getSyst_CV1uBY();
760  //std::vector<const ISyst*> systs = getSyst_CV2uBY();
761  //std::vector<const ISyst*> systs = getAllXsecSysts_2017();
762  //std::vector<const ISyst*> systs = {&kMAQAGenieReducedSyst2017};
763 
764  // --- Make a vector of "Selections" which are a combination of Quantile && Tiered cuts.
765  std::vector<Selection> selections;
766  for (size_t cuts=0; cuts < TieredCuts.size(); ++cuts) {
767  for (size_t quant=0; quant < HadEFracQuantCuts.size(); ++quant) {
768  // --- Want no quantile cut selections too.
769  if (quant == 0) {
770  selections.emplace_back( TieredCuts[cuts], CutNames[cuts]+"_QuantAll",
771  "Selected events pass "+CutNames[cuts]+" all quantiles");
772  }
773  selections.emplace_back( TieredCuts[cuts] && HadEFracQuantCuts[quant], CutNames[cuts]+"_Quant"+std::to_string(quant+1),
774  "Selected events pass "+CutNames[cuts]+" quantile "+std::to_string(quant+1)+".");
775  }
776  }
777 
778  // --- Make a vector of "TangibleAxis" which have the axis and variables for which I want to make spectra.
779  std::vector<TangibleAxis> variables;
780 
781  if ( RunBinNorm ) {
782  variables.emplace_back( kNumuCCOptimisedAxis, "NuMuE" , "CC energy estimator.");
783  } else {
784  variables.emplace_back( EnergyAxis , "MuonE" , "Reconstructed energy of primary muon track.");
785  variables.emplace_back( HadEAxis , "HadE" , "Hadronic enery, i.e. numu energy estimate minus muon track energy. " );
786  variables.emplace_back( HadEFracAxis , "HadEFrac" , "Hadronic enery fraction, i.e. ration of hadronic energy to numu energy. " );
787  variables.emplace_back( ReMIDAxis , "ReMIDScore", "ReMId kNN score." );
788  variables.emplace_back( CVNNuMuAxis , "CVNNuMu" , "CVN muon identification score." );
789  }
790 
791  //variables.emplace_back( LengthAxis , "TrkLength" , "Primary track length." );
792  //variables.emplace_back( CosNuMiAxis , "CosNuMi" , "Beam direction of muon track. " );
793 
794 
795  // // // copied from sa script
796  // // variables.emplace_back(HistAxis("Slice N_{Hit}", Binning::Simple(50, 0, 500), kNHit),
797  // // "slcNHit", "Number of hits in slice. " );
798 
799  // // variables.emplace_back(HistAxis("Average Hadronic Energy Per Hit (GeV)", Binning::Simple(40, 0, 0.04), kHadEPerNHit),
800  // // "hadEPerNHit", "Average energy per hit in hadronic cluster, i.e. had_E/had_n_hit. " );
801 
802  // // variables.emplace_back(HistAxis("Track Energy Per Hit (GeV)", Binning::Simple(40, 0, 0.04), kTrkEPerNHit),
803  // // "trkEPerNHit", "Average energy per hit on primary track, i.e. trk_E/trk_n_hit. " );
804 
805  // // variables.emplace_back(HistAxis("Hadronic N_{Hit}", Binning::Simple(50, 0, 100), kHadNHit),
806  // // "hadNHit", "Number of hits in hadronic cluster. ");
807 
808 
809  // variables.emplace_back(HistAxis("Slice Maximum Y (m)", kXYBins, kSlcMaxY),
810  // "maxy", "Maximum Y position of slice. " );
811  // variables.emplace_back(HistAxis("Number of Tracks in Slice", Binning::Simple(14, 1, 15), kNKalman),
812  // "nkal", "Number of tracks in slice. " );
813  // // variables.emplace_back(HistAxis("Number of Hits in Slice", Binning::Simple(50, 0, 500), kNHit),
814  // // "nhit", "Number of hits in slice. " );
815  // variables.emplace_back(HistAxis("Muon Track cos(#theta_{Z})", Binning::Simple(50, 0,1), kDirZ),
816  // "dirZ", "Z-direction of muon track. " );
817  // variables.emplace_back(HistAxis("Track Start X Position (m)", kXYBins, kTrkStartX),
818  // "trkStartX", "Track start x position. " );
819  // variables.emplace_back(HistAxis("Track Start Y Position (m)", kXYBins, kTrkStartY),
820  // "trkStartY", "Track start y position. " );
821  // variables.emplace_back(
822  // HistAxis("Track Start Z Position (m)", kZBins, kTrkStartZ),
823  // "trkStartZ", "Track start z position. " );
824  // variables.emplace_back(
825  // HistAxis("Track End X Position (m)", kXYBins, kTrkEndX),
826  // "trkEndX", "Track stop x position. " );
827  // variables.emplace_back(
828  // HistAxis("Track End Y Position (m)", kXYBins, kTrkEndY),
829  // "trkEndY", "Track stop y position. " );
830  // variables.emplace_back(
831  // HistAxis("Track End Z Position (m)", kZBins, kTrkEndZ),
832  // "trkEndZ", "Track stop z position. " );
833  // // variables.emplace_back(
834  // // HistAxis("Slice Duration [ns]", Binning::Simple(50,0,3000), kSliceDuration),
835  // // "sliceDuration", "Slice duration. " );
836  // // variables.emplace_back(
837  // // HistAxis("Number of Hits in Primary Track", Binning::Simple(50,0,500), kTrkNhits),
838  // // "trkNhits", "Number of hits on primary track. ");
839  // // variables.emplace_back(
840  // // HistAxis("Length of Primary Track (m)", Binning::Simple(50,0,16), kTrkLength),
841  // // "trkLength", "Primary track length. " );
842 
843  // variables.emplace_back(
844  // HistAxis("Scattering Log-likelihood", Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
845  // "scatLL", "ReMId scattering log log-likelihood for primary track. " );
846  // variables.emplace_back(
847  // HistAxis("dE/dx Log-likelihood", Binning::Simple(50,-3,1), kReMIdDEDxLLH),
848  // "dedxLL", "ReMId dE/dx log log-likelihood for primary track. " );
849  // variables.emplace_back(
850  // HistAxis("Non-hadronic Plane Fraction",
851  // Binning::Simple(50,0,1), kReMIdMeasFrac),
852  // "nonHadPlaneFrac", "ReMId Non-hadronic plane fraction. " );
853 
854  // // finished copy from sa script
855 
856  // --- What is my weight going to be for the spectra?
857  Var MyWeight = kPPFXFluxCVWgt * kXSecCVWgt2017;
858  if (RunBinNorm) MyWeight = kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight;
859 
860  std::vector<DataMCPair> pairs;
861 
862  pairs.reserve(selections.size() * variables.size());
863  for(const auto& sel:selections){
864  for(const auto& variable:variables){
865  pairs.emplace_back(sel, variable, loaderData, loaderMC, systs, MyWeight, !kIsNumuCC );
866  }
867  }
868 
869  loaderMC.Go();
870  loaderData.Go();
871 
872  std::vector<TCanvas*> cans;
873 
874  for(const auto& pair:pairs) {
875 
876  std::cout << "\nNow looking pairs, " << pair.ShortName() << " - " << pair.CName() << ", it has purity " << pair.Purity() << std::endl;
877  TString name(TString(pair.CName()) + "_allSysts_pot");
878 
879  //cans.push_back(new TCanvas(name,name,1200,1000));
880  cans.push_back(new TCanvas(name,name));
881  pair.OverlayDataMCSyst();
882  WriteBlurb(pair, "nd_syst_plots/",
883  std::string(cans.back()->GetName())+".txt",
884  "Full systematic band shown. ",
885  "Near detector data/MC comparison for numu first analysis. " );
886  /*
887  cans.push_back(LogClone(cans.back()));
888  WriteBlurb(pair, "nd_syst_plots/",
889  std::string(cans.back()->GetName())+".txt",
890  "Full systematic band shown. Semi-log scale.",
891  "Near detector data/MC comparison for numu first analysis. " );
892  */
893 
894  // Other plot used by Vlad.
895  //*
896  TString nameNorm(TString(pair.CName()) + "_allSysts_area");
897  //cans.push_back(new TCanvas(nameNorm,nameNorm,1200,1000));
898  cans.push_back(new TCanvas(nameNorm,nameNorm));
899  pair.OverlayDataMCSystNorm();
900  WriteBlurb(pair, "nd_syst_plots/",
901  std::string(cans.back()->GetName())+".txt",
902  std::string("Each systematically shifted histogram (both up and down) ")
903  +std::string("is normalized to the area of the MC distribution, ")
904  +std::string("then summed in quadrature." ),
905  "Near detector data/MC comparison for numu first analysis. " );
906 
907  //*/
908 
909  /*
910  cans.push_back(LogClone(cans.back()));
911  WriteBlurb(pair, "nd_syst_plots/",
912  std::string(cans.back()->GetName())+".txt",
913  std::string("Each systematically shifted histogram (both up and down) ")
914  +std::string("is normalized to the area of the MC distribution, ")
915  +std::string("then summed in quadrature. Semi-log scale." ),
916  "Near detector data/MC comparison for numu first analysis. ");
917  */
918  }
919 
920  // Open a file to save all of the canvases to.
921  TFile *OutFile = new TFile("nd_syst_plots/NearDet_NuMiData_Systs_NuMuAna17.root", "RECREATE");
922  OutFile->cd();
923 
924  for(const auto& can:cans) {
925  can->cd();
926  Preliminary();
927  //can -> SetLeftMargin(0.13);
928  can->SaveAs (TString("nd_syst_plots/") + can->GetName() + ".pdf");
929  can->Write ( can->GetName() );
930  MakeTextFile( can->GetName() );
931  }
932 }
const Var kHadE
Definition: NumuVars.h:23
const XML_Char * name
Definition: expat.h:151
DataMCPair * findPair(std::vector< DataMCPair > &list, std::vector< std::string > keys)
Get a vector of all the numu group systs.
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
keys
Reco plots.
Definition: caf_analysis.py:46
enum BeamMode kRed
void DrawData(const int color=kBlack, const int marker=kFullCircle) const
void DrawMCBkg() const
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kNumuContainND2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100 &&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
Definition: NumuCuts2017.h:11
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:233
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
constexpr T pow(T x)
Definition: pow.h:75
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
std::string fShortName
Definition: DataMCPair.h:47
void Preliminary()
Tangible< Cut > Selection
Definition: DataMCPair.h:52
std::string ShortName() const
void SetSpillCut(const SpillCut &cut)
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
std::vector< Spectrum > fUps
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
cout<< t1-> GetEntries()
Definition: plottest35.C:29
void OverlayDataMCSystExtraData(const DataMCPair &extraData, const std::string dataLegendTitle, const std::string extraLegendTitle, const int iSyst=-1, const int dataColor=kBlack, const int extraColor=kGray+1, const int dataMarker=kFullCircle, const int extraMarker=kFullCircle) const
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void AddExposureArea() const
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:304
void DrawMCSyst(const int iSyst=-1) const
ifstream inFile
Definition: AnaPlotMaker.h:34
const Var kBinWidthWeight([](const caf::SRProxy *sr){const double Enu=kCCE(sr); if(Enu< 0.75) return 0.1/0.75;if(Enu >=0.75 &&Enu< 1) return 0.1/0.25;if(Enu >=1 &&Enu< 2) return 0.1/0.1;if(Enu >=2 &&Enu< 3) return 0.1/0.25;if(Enu >=3 &&Enu< 4) return 0.1/0.5;if(Enu >=4 &&Enu< 5) return 0.1/1;else return 1.;})
const Var kRemID
PID
Definition: Vars.cxx:81
void DrawData(const int color=kBlack, EBinType bintype=kBinContent) const
Draw data on plots, mostly for internal use.
Definition: DataMCPair.cxx:214
Tangible(const T &obj, const std::string &shortName, const std::string &blurb)
std::vector< Spectrum > fDowns
void RedrawAxes()
const Var kHadEFrac
Definition: NumuVars.h:24
Tangible< HistAxis > TangibleAxis
Definition: DataMCPair.h:53
const Var kCCE
Definition: NumuVars.h:21
#define pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
DataMCPair(Selection sel, TangibleAxis tanAxis, SpectrumLoader &loaderData, SpectrumLoader &loaderMC, std::vector< const ISyst * > systs, const Var BinWeight, const Cut &bkg=kNoCut)
const Cut kNumuPID2017([](const caf::SRProxy *sr){return(sr->sel.remid.pid > 0.5 &&sr->sel.cvn.numuid > 0.5);})
Definition: NumuCuts2017.h:27
const char * CName() const
TCanvas * LogClone(const TCanvas *can)
TLatex * prelim
Definition: Xsec_final.C:133
double POT() const
Definition: Spectrum.h:227
list cans
Definition: canMan.py:12
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
TLegend * DrawLegendArea(const int dataColor=kBlack, std::string dataLegendTitle="Data", const int dataMarker=kFullCircle) const
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
std::string fBlurb
Definition: DataMCPair.h:46
void hadEFrac_nd_data_mc_systs(bool RunBinNorm=false, bool IsTest=false)
std::string fShortName
Selection fSel
Definition: DataMCPair.h:182
void AestheticsPOT() const
void AestheticsArea() const
Float_t norm
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
T log10(T number)
Definition: d0nt_math.hpp:120
TDirectory * dir
Definition: macro.C:5
TLegend * DrawLegendPOT(const int dataColor=kBlack, std::string dataLegendTitle="Data", const int dataMarker=kFullCircle) const
void OverlayDataMCSystNorm(const int iSyst=-1) const
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
assert(nhit_max >=nhit_nbins)
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
float Purity() const
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
std::vector< const ISyst * > getAllDirectNumuSysts2017()
Get a vector of all the numu group systs.
Definition: NumuSysts.cxx:169
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Color_t kTotalMCColor
Definition: Style.h:16
TFile * OutFile
double T
Definition: Xdiff_gwt.C:5
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const Cut kNumuQuality
Definition: NumuCuts.h:18
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
const Var kMuE
Definition: NumuVars.h:22
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void MakeTextFile(std::string OutName)
void DrawMCNormSyst(const int iSyst=-1) const
enum BeamMode kBlue
const Var kCVNm
PID
Definition: Vars.cxx:39
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
def entry(str)
Definition: HTMLTools.py:26
short int place(double jd_tt, object *cel_object, observer *location, double delta_t, short int coord_sys, short int accuracy, sky_pos *output)
void AddExposurePOT() const
void OverlayDataMCSyst(const int iSyst=-1) const
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
const Var kCosNumi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);}}return-5.f;})
Definition: NumuVars.h:43
void OverlayDataMCSystNormExtraData(const DataMCPair &extraData, const std::string dataLegendTitle, const std::string extraLegendTitle, const int iSyst=-1, const int dataColor=kBlack, const int extraColor=kGray+1, const int dataMarker=kFullCircle, const int extraMarker=kFullCircle) const
void WriteBlurb(const DataMCPair &pair, std::string dir, std::string name, std::string extra="", std::string beginning="")
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49
enum BeamMode string