DataMCPair.cxx
Go to the documentation of this file.
2 
3 
4 #include "CAFAna/Core/Binning.h"
5 #include "CAFAna/Core/Var.h"
7 
10 
11 #include "TDirectory.h"
12 #include "TH1.h"
13 #include "TH2.h"
14 #include "TKey.h"
15 #include "TLatex.h"
16 #include "TLegend.h"
17 #include "TList.h"
18 #include "TObjString.h"
19 #include "TParameter.h"
20 
21 #include <iostream>
22 
23 namespace ana
24 {
25 
26  // Replace all instances of from with to in str
28  const std::string& from,
29  const std::string& to) {
30  std::string ret(str);
31  if(from.empty())
32  return ret;
33  size_t start_pos = 0;
34  while((start_pos = ret.find(from, start_pos)) != std::string::npos) {
35  ret.replace(start_pos, from.length(), to);
36  start_pos += to.length(); // In case 'to' contains 'from', like replacing 'x' with 'yx'
37  }
38  return ret;
39  }
40 
41 
42  TString TableNum(float num)
43  {
44  return TString::Format("%.02lf", num);
45  }
46 
47  //-------------------------------------------------------------------------
49  SpectrumLoaderBase& loaderData,
50  SpectrumLoaderBase& loaderMC,
51  const Var & wei,
52  const std::vector<const ISyst*>& systs,
53  const std::map<std::string, SpectrumComponents::Component>& mcCategories)
54  :
55  fSel(sel),
56  fData(loaderData, tanAxis.fObj, fSel.fObj, kNoShift),
57  fMC(fSel.fObj, tanAxis.fObj, loaderMC, wei, mcCategories),
58  fShortName(tanAxis.fShortName + "_" + fSel.fShortName),
59  fSysts(systs)
60  {
61  fUps.reserve(systs.size());
62  fDowns.reserve(systs.size());
63  for(const auto& syst:systs){
64  fUps.emplace_back(loaderMC, tanAxis.fObj, fSel.fObj, SystShifts(syst,+1), wei);
65  fDowns.emplace_back(loaderMC, tanAxis.fObj, fSel.fObj, SystShifts(syst,-1), wei);
66  }
67  }
68 
69  //-------------------------------------------------------------------------
70 
72  SpectrumLoaderBase& loaderData,
73  SpectrumLoaderBase& loaderMC,
74  const Var & wei,
75  const std::vector<const ISyst*>& systs,
76  const std::map<std::string, SpectrumComponents::Component>& mcCategories)
77  : DataMCPair(Selection(cut, "", ""), TangibleAxis(axis, "", ""),
78  loaderData, loaderMC, wei, systs, mcCategories)
79  {
80  }
81 
82 
83  //-------------------------------------------------------------------------
84  void DataMCPair::DrawMCSyst(const int iSyst, EBinType bintype) const
85  {
86 
87  if(iSyst < 0)
89  kFullPredColor, kPredErrColor, 1.5, true, bintype);
90  else
91  {
92  assert(iSyst < (int)fUps.size());
93  PlotWithSystErrorBand(fMC.Tot(), {fUps[iSyst]}, {fDowns[iSyst]}, fData.POT(),
94  kFullPredColor, kPredErrColor, 1.5, true, bintype);
95  }
96  }
97 
98  //-------------------------------------------------------------------------
99  void DataMCPair::OverlayDataMCSyst(const int iSyst, bool drawMCBreakdown, EBinType bintype, bool drawLegend) const
100  {
101  DrawMCSyst(iSyst, bintype);
102  if (drawMCBreakdown)
103  DrawMCComponents(true, bintype);
104  DrawData(kBlack, bintype);
105  if (drawLegend)
106  DrawLegend();
107  }
108 
109  //-------------------------------------------------------------------------
110  void DataMCPair::DrawMCNormSyst(const int iSyst, EBinType bintype) const
111  {
112  // New vectors for area normalized systs
113  std::vector<Spectrum> normUps;
114  std::vector<Spectrum> normDowns;
115 
116  if (iSyst >= 0)
117  {
118  // We'll have one entry if there's just one syst to plot.
119  normUps.reserve(1);
120  normDowns.reserve(1);
121  }
122  else
123  {
124  // Otherwise we need all systs.
125  normUps.reserve(fUps.size());
126  normDowns.reserve(fUps.size());
127  }
128  for(size_t i = 0; i < fUps.size(); ++i)
129  {
130  // bail if we want a particular syst and it's not the one we want.
131  if(iSyst >= 0 && (int)i != iSyst) continue;
132  normUps.emplace_back(fUps[i]);
133  Spectrum& norm = normUps.back();
134  norm.OverridePOT(norm.POT() * norm.Integral(1) / fMC.Tot().Integral(1));
135  }
136  for(size_t i = 0; i < fDowns.size(); ++i)
137  {
138  // bail if we want a particular syst and it's not the one we want.
139  if(iSyst >= 0 && (int)i != iSyst) continue;
140  normDowns.emplace_back(fDowns[i]);
141  Spectrum& norm = normDowns.back();
142  norm.OverridePOT(norm.POT() * norm.Integral(1) / fMC.Tot().Integral(1));
143  }
144  PlotWithSystErrorBand(fMC.Tot(), normUps, normDowns, fData.POT(),
145  kFullPredColor, kPredErrColor, 1.5,true, bintype);
146  }
147 
148  //-------------------------------------------------------------------------
149  void DataMCPair::OverlayDataMCSystNorm(const int iSyst, bool drawMCBreakdown, EBinType bintype, bool drawLegend) const
150  {
151  DrawMCNormSyst(iSyst, bintype);
152  if (drawMCBreakdown)
153  DrawMCComponents(true, bintype);
154  DrawData(kBlack, bintype);
155  if (drawLegend)
156  DrawLegend();
157  }
158 
159  //-------------------------------------------------------------------------
160  TLegend* DataMCPair::DrawLegend(const int dataColor, bool includeComponents, double drawThreshold) const
161  {
162  TLegend* leg;
163  if (includeComponents)
164  // needs to be a bit smaller if components included
165  leg = AutoPlaceLegend(0.25, 0.18 + 0.04*fMC.ComponentDescrs().size());
166  else
167  leg = AutoPlaceLegend(0.36, 0.18, 0.67);
168  //TLegend* leg = new TLegend(0.6, 0.7, 0.85, 0.85);
169  TH1F* colMC = new TH1F();
170  colMC->SetLineColor(kFullPredColor);
171  TH1F* colMCBkg = new TH1F();
172  colMCBkg->SetLineColor(kTotalBkgColor);
173  TH1F* colData = new TH1F();
174  colData->SetLineColor(dataColor);
175  colData->SetMarkerColor(dataColor);
176 
177  leg->AddEntry((TObject*)colData, "Data", "ple");
178  leg->AddEntry((TObject*)colMC, (includeComponents ? "Total Simulation" : "Simulation"), "l");
179 
180  if (includeComponents)
181  fMC.DrawLegend(drawThreshold, leg);
182 
183 
184  leg->SetFillColor(0);
185  leg->SetFillStyle(0);
186 
187  for(const auto& obj:*leg->GetListOfPrimitives())
188  {
189  if(obj->InheritsFrom("TAttFill")){
190  ((TAttFill*)obj)->SetFillStyle(0);
191  //((TAttFill*)obj)->SetFillColor(0);
192  }
193 
194  }
195  leg->Draw();
196  return leg;
197  }
198 
199  //-------------------------------------------------------------------------
201  {
202  std::stringstream expo;
203  float pot = fData.POT();
204  expo.precision(2);
205  expo << "ND, "
206  << TString::Format("%.02lf#times10^{20}", pot/1e20) << " POT";
207  TLatex text;
208  text.SetTextAlign(12);
209 
210  text.DrawLatexNDC(0.2, 0.95, expo.str().c_str());
211  }
212 
213  //-------------------------------------------------------------------------
214  void DataMCPair::DrawData(const int color, EBinType bintype) const
215  {
216  TH1D* hData = fData.ToTH1(fData.POT(), color, kFullCircle, kPOT, bintype);
217 
218  hData->Draw("SAME");
219  }
220 
221  //-------------------------------------------------------------------------
222  void DataMCPair::DrawMCComponents(bool stacked, EBinType bintype) const
223  {
224  fMC.DrawComponents(stacked, bintype, fData.POT());
225  }
226 
227  //-------------------------------------------------------------------------
228  float DataMCPair::Purity(const std::set<std::string>& signalCatNames) const
229  {
230  return fMC.Purity(signalCatNames);
231  }
232 
233  //-------------------------------------------------------------------------
235  {
236 
237  std::ofstream table(fname);
238  table <<
239  "\\documentclass[12pt]{article}\n"
240  "\\usepackage{longtable}\n"
241  "\\begin{document}\n"
242  << fShortName <<
243  "\\begin{longtable}\n"
244  "{|p{0.4\\textwidth}|p{0.3\\textwidth}|p{0.3\\textwidth}|}\n"
245  "\\hline\n"
246  "Parameter & $+1 \\sigma$ Mean & $-1 \\sigma$ Mean \\\\\n"
247  "\\hline \\hline\n"
248  "\\endfirsthead\n"
249  "\\multicolumn{3}{l}\n"
250  "{{\\bfseries\\tablename\\ \\thetable{} -- "
251  "continued from previous page}}\\\\\n"
252  "\\hline\n"
253  "Parameter & $+1 \\sigma$ Mean & $-1 \\sigma$ Mean \\\\\n"
254  "\\hline \\hline\n"
255  "\\endhead\n"
256  "\\hline \\hline \\multicolumn{3}{|r|}{{Continued on next page}}"
257  " \\\\ \\hline\n"
258  "\\endfoot\n"
259  "\\hline\n"
260  "\\endlastfoot\n"
261  "\\hline\n"
262  << std::endl;
263 
264  float meanNom = fMC.Tot().Mean();
265  for(size_t iSyst = 0; iSyst < fSysts.size(); ++ iSyst)
266  {
267 
268  float meanUp = fUps[iSyst].Mean();
269  float meanDown = fDowns[iSyst].Mean();
270  float percUp = 100 * (meanUp - meanNom)/meanNom;
271  float percDown = 100 * (meanDown - meanNom)/meanNom;
272 
273  table << replaceAll(fSysts[iSyst]->LatexName(), "_", "\\_") << " & "
274  << TableNum(meanUp) << " (" << TableNum(percUp) << "\\%) &"
275  << TableNum(meanDown) << " (" << TableNum(percDown) << "\\%) "
276  << " \\\\ \\hline" << std::endl;
277  }
278  table << "\\end{longtable}\n\\end{document}" << std::endl;
279  }
280 
281  //-------------------------------------------------------------------------
282  void DataMCPair::SaveTo(TDirectory* dir, const std::string& name) const
283  {
284  TDirectory* tmp = gDirectory;
285 
286  dir = dir->mkdir(name.c_str()); // switch to subdir
287  dir->cd();
288 
289  fData.SaveTo(dir, "data");
290 
291  fMC.SaveTo(dir, "MC");
292 
293  TDirectory * upsDir = dir->mkdir("upShifts");
294  for (std::size_t idx = 0; idx < fUps.size(); idx++)
295  {
296  fUps[idx].SaveTo(upsDir, std::to_string(idx));
297  }
298  TDirectory * downsDir = dir->mkdir("downShifts");
299  for (std::size_t idx = 0; idx < fDowns.size(); idx++)
300  {
301  fDowns[idx].SaveTo(downsDir, std::to_string(idx));
302  }
303 
304  dir->cd();
305  TObjString shortName(fShortName.c_str());
306  shortName.Write("shortName");
307 
308  dir->Write();
309  delete dir;
310 
311  tmp->cd();
312  }
313 
314  //-------------------------------------------------------------------------
315  std::unique_ptr<DataMCPair> DataMCPair::LoadFrom(TDirectory* dir, const std::string& name)
316  {
317  dir = dir->GetDirectory(name.c_str()); // switch to subdir
318  assert(dir);
319 
320  Spectrum data = *Spectrum::LoadFrom(dir, "data");
322 
323  std::map<std::string, std::vector<Spectrum>> specs;
324  for (const auto & shift : {"up", "down"})
325  {
326  TDirectory * shiftsDir = dynamic_cast<TDirectory*>(dir->Get(Form("%sShifts", shift)));
327  // maybe there weren't any shifts, and a hadd didn't create the empty dir?
328  if (!shiftsDir)
329  continue;
330 
331  for ( TObject * obj : *(shiftsDir->GetListOfKeys()) )
332  {
333  TKey * key = dynamic_cast<TKey*>(obj);
334 
335  // pray that they're read out in the same order they went in
336  specs[shift].emplace_back( *(Spectrum::LoadFrom(shiftsDir, key->GetName()).release()) );
337  }
338  }
339  auto ret = std::make_unique<DataMCPair>(
340  std::move(data),
341  std::move(mc),
342  std::move(specs["up"]),
343  std::move(specs["down"]),
344  dynamic_cast<TObjString*>(dir->Get("shortName"))->String().Data()
345  );
346 
347  delete dir;
348 
349  return std::move(ret);
350  }
351 
352  //-------------------------------------------------------------------------
353 
355  {
356  assert (fMC.ComponentDescrs().find(name) != fMC.ComponentDescrs().end());
357 
358  fMC.ComponentDescrs().at(name).color = color;
359  }
360 
361  //-------------------------------------------------------------------------
362 
364  {
365  assert (fMC.ComponentDescrs().find(name) != fMC.ComponentDescrs().end());
366 
367  fMC.ComponentDescrs().at(name).blurb = blurb;
368 
369  }
370 
371 
372 } // namespace
TString TableNum(float num)
Format numbers for tables, with two nice digits after decimal.
Definition: DataMCPair.cxx:42
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string fShortName
Definition: DataMCPair.h:187
static std::unique_ptr< SpectrumComponents > LoadFrom(TDirectory *dir, const std::string &name)
Spectrum fData
Definition: DataMCPair.h:183
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
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
void DrawMCSyst(const int iSyst=-1, EBinType bintype=kBinContent) const
Draw MC with error band.
Definition: DataMCPair.cxx:84
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:225
void DrawMCNormSyst(const int iSyst=-1, EBinType bintype=kBinContent) const
Draw MC with error band, but with each syst normalized to total area of MC to highlight shape differe...
Definition: DataMCPair.cxx:110
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: DataMCPair.cxx:282
void DrawMCComponents(bool stacked=true, EBinType bintype=kBinContent) const
Draw MC components distribution.
Definition: DataMCPair.cxx:222
void OverlayDataMCSystNorm(const int iSyst=-1, bool drawMCBreakdown=true, EBinType bintype=kBinContent, bool drawLegend=true) const
Draw data and MC with error band, but with each syst normalized to total area of MC to highlight shap...
Definition: DataMCPair.cxx:149
std::vector< Spectrum > fDowns
Definition: DataMCPair.h:186
Float_t tmp
Definition: plot.C:36
float Purity(const std::set< std::string > &signalCatNames={}) const
Purity of the MC selection based on the MC subcategories.
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:272
void SetComponentColor(const std::string &name, int color)
Definition: DataMCPair.cxx:354
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const std::vector< const ISyst * > fSysts
Definition: DataMCPair.h:188
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
const XML_Char const XML_Char * data
Definition: expat.h:268
EBinType
Definition: UtilsExt.h:28
const Color_t kFullPredColor
Definition: Style.h:36
DataMCPair(Selection sel, TangibleAxis tanAxis, SpectrumLoaderBase &loaderData, SpectrumLoaderBase &loaderMC, const Var &wei=kUnweighted, const std::vector< const ISyst * > &systs={}, const std::map< std::string, SpectrumComponents::Component > &mcCategories={})
DataMCPair wraps up spectrum creation for data/MC comparison
Definition: DataMCPair.cxx:48
void DrawData(const int color=kBlack, EBinType bintype=kBinContent) const
Draw data on plots, mostly for internal use.
Definition: DataMCPair.cxx:214
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
const Color_t kPredErrColor
Definition: Style.h:37
void AddExposure() const
Write exposure on plots, used internally.
Definition: DataMCPair.cxx:200
TLegend * DrawLegend(double drawThreshold=0, TLegend *leg=nullptr) const
Draw legend on plots.
#define pot
const Color_t kTotalBkgColor
Definition: Style.h:39
void DrawComponents(bool stacked=true, EBinType bintype=kBinContent, double POT=0, bool sameAll=true) const
Draw MC components distribution.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
float Purity(const std::set< std::string > &signalCatNames={}) const
Purity of the MC selection based on the MC subcategories.
Definition: DataMCPair.cxx:228
double POT() const
Definition: Spectrum.h:219
void OverlayDataMCSyst(const int iSyst=-1, bool drawMCBreakdown=true, EBinType bintype=kBinContent, bool drawLegend=true) const
Draw data and MC with error band.
Definition: DataMCPair.cxx:99
const SystShifts kNoShift
Definition: SystShifts.cxx:21
void CreateSystTable(const std::string &fname) const
Write LaTeX table of syst means: |par|mean up(%)|mean down(%)|.
Definition: DataMCPair.cxx:234
Base class for the various types of spectrum loader.
static std::unique_ptr< DataMCPair > LoadFrom(TDirectory *dir, const std::string &name)
Definition: DataMCPair.cxx:315
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
std::vector< Spectrum > fUps
Definition: DataMCPair.h:185
Selection fSel
Definition: DataMCPair.h:182
Float_t norm
TDirectory * dir
Definition: macro.C:5
int num
Definition: f2_nu.C:119
TLegend * DrawLegend(const int dataColor=kBlack, bool includeComponents=true, double drawThreshold=-1) const
Draw legend on plots, mostly for internal use.
Definition: DataMCPair.cxx:160
std::string replaceAll(const std::string &str, const std::string &from, const std::string &to)
Definition: DataMCPair.cxx:27
assert(nhit_max >=nhit_nbins)
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void SetComponentBlurb(const std::string &name, const std::string &blurb)
Definition: DataMCPair.cxx:363
SpectrumComponents fMC
Definition: DataMCPair.h:184
const std::map< std::string, Component > & ComponentDescrs() const
void SaveTo(TDirectory *dir, const std::string &name) const
double Mean() const
Return mean of 1D histogram.
Definition: Spectrum.cxx:290