SpectrumComponents.cxx
Go to the documentation of this file.
2 
3 #include <map>
4 #include <string>
5 
6 #include "TDirectory.h"
7 #include "TH1F.h"
8 #include "TKey.h"
9 #include "TLegend.h"
10 #include "TObjString.h"
11 #include "TParameter.h"
12 
13 #include "CAFAna/Analysis/Plots.h"
14 #include "CAFAna/Core/Spectrum.h"
15 
16 namespace ana
17 {
18  //-------------------------------------------------------------------------
20  const ana::HistAxis & axis,
22  const Var & wei,
23  const std::map<std::string, Component>& categories)
24  :
25  fSel(baseCut),
26  fTot(loader, axis, fSel, kNoShift, wei)
27  {
28  if (!categories.empty())
29  {
30  for (const auto & componentPair : categories)
31  {
32  fComponentDescr.emplace(componentPair.first, componentPair.second);
33  fComponentSpecs.emplace(std::piecewise_construct,
34  std::forward_as_tuple(componentPair.first),
35  std::forward_as_tuple(loader, axis, fSel && componentPair.second.cut, kNoShift, wei));
36  }
37  }
38  }
39 
40  //-------------------------------------------------------------------------
44  const Var & wei,
45  const std::map<std::string, Component>& categories)
46  :
47  fSel(baseCut),
48  fTot(loader, axis, fSel, kNoShift, wei)
49  {
50  if (!categories.empty())
51  {
52  for (const auto & componentPair : categories)
53  {
54  fComponentDescr.emplace(componentPair.first, componentPair.second);
55  fComponentSpecs.emplace(std::piecewise_construct,
56  std::forward_as_tuple(componentPair.first),
57  std::forward_as_tuple(loader, axis, fSel && componentPair.second.cut, kNoShift, wei));
58  }
59  }
60  }
61 
62  //-------------------------------------------------------------------------
63  TLegend* SpectrumComponents::DrawLegend(double drawThreshold, TLegend * leg) const
64  {
65  if (!leg)
66  leg = AutoPlaceLegend(0.25, 0.18 + 0.04*fComponentDescr.size());
67 
68  double thresh = 0;
69  if (drawThreshold >= 0 && drawThreshold <= 1.0)
70  thresh = fTot.Integral(fTot.POT()) * drawThreshold;
71  else
72  // this is supposed to be a fraction
73  assert(drawThreshold <= 1 && "SpectrumComponents::DrawLegend(): legend draw thresholds must be between 0.0 and 1.0");
74  for (const auto & mcSpecPair : fComponentDescr)
75  {
76  auto spec = fComponentSpecs.at(mcSpecPair.first);
77  if (spec.Integral(spec.POT()) < thresh)
78  continue;
79 
80  TH1F* comp = new TH1F();
81  comp->SetLineColor(kBlack);
82  comp->SetFillColor(mcSpecPair.second.color);
83  leg->AddEntry((TObject*)comp, mcSpecPair.second.blurb.c_str(), "f");
84  }
85 
86 
87  leg->SetFillColor(0);
88  leg->SetFillStyle(0);
89 
90  for(const auto& obj:*leg->GetListOfPrimitives())
91  {
92  if(obj->InheritsFrom("TAttFill")){
93  ((TAttFill*)obj)->SetFillStyle(0);
94  //((TAttFill*)obj)->SetFillColor(0);
95  }
96 
97  }
98  leg->Draw();
99  return leg;
100  }
101 
102  //-------------------------------------------------------------------------
103  void SpectrumComponents::DrawComponents(bool stacked, EBinType bintype, double POT, bool sameAll) const
104  {
105  if (POT == 0)
106  POT = fComponentSpecs.begin()->second.POT();
107 
108  TH1 * last = nullptr;
109  std::vector<TH1*> hists;
110 
111  // iterate backwards so that they're drawn in the specified order from top to bottom
112  for (auto it_comp = fComponentSpecs.rbegin(); it_comp != fComponentSpecs.rend(); it_comp++)
113  {
114  auto h = it_comp->second.ToTH1(POT, kBlack, kSolid, kPOT, bintype);
115  if (stacked)
116  {
117  if (last)
118  h->Add(last);
119  h->SetFillColor(fComponentDescr.at(it_comp->first).color);
120  h->SetLineWidth(1);
121  last = h;
122  }
123 
124  hists.push_back(h);
125  }
126  bool drawn = false;
127  for (auto it_hist = hists.rbegin(); it_hist != hists.rend(); it_hist++)
128  {
129  std::string opt = "hist";
130  if (drawn || sameAll)
131  opt += " same";
132  (*it_hist)->Draw(opt.c_str());
133  if (!drawn)
134  drawn = true;
135  }
136  }
137 
138  //-------------------------------------------------------------------------
139  float SpectrumComponents::Purity(const std::set<std::string>& signalCatNames) const
140  {
141  double bkgCounts = 0;
142  double POT = fComponentSpecs.begin()->second.POT();
143  for (const auto & mcCatPair : fComponentSpecs)
144  {
145  if (signalCatNames.find(mcCatPair.first) != signalCatNames.end())
146  continue;
147  bkgCounts += mcCatPair.second.ToTH1(POT)->GetEntries();
148  }
149  return 1 - bkgCounts / fTot.ToTH1(POT)->GetEntries();
150  }
151 
152  //-------------------------------------------------------------------------
153  void SpectrumComponents::SaveTo(TDirectory* dir, const std::string& name) const
154  {
155  TDirectory* tmp = gDirectory;
156 
157  dir = dir->mkdir(name.c_str()); // switch to subdir
158  dir->cd();
159 
160  fTot.SaveTo(dir, "Total");
161 
162  TDirectory * compsDir = dir->mkdir("Components");
163  for (const auto & specPair : fComponentSpecs)
164  {
165  specPair.second.SaveTo(compsDir, specPair.first);
166 
167  TDirectory* compDir = compsDir->GetDirectory(specPair.first.c_str());
168  compDir->cd();
169  const auto & comp = fComponentDescr.at(specPair.first);
170  TObjString blurb(comp.blurb.c_str());
171  blurb.Write("blurb");
172 
173  TParameter<int> color;
174  color.SetVal(comp.color);
175  color.Write("color");
176  }
177 
178  dir->Write();
179  delete dir;
180 
181  tmp->cd();
182  }
183 
184  //-------------------------------------------------------------------------
185  std::unique_ptr<SpectrumComponents> SpectrumComponents::LoadFrom(TDirectory* dir, const std::string& name)
186  {
187  dir = dir->GetDirectory(name.c_str()); // switch to subdir
188  assert(dir);
189 
190  Spectrum tot = *Spectrum::LoadFrom(dir, "Total");
191 
192  std::map<std::string, Component> compDescrs;
193  std::map<std::string, Spectrum> compSpecs;
194  TDirectory * compDir = dynamic_cast<TDirectory*>(dir->Get("Components"));
195  for ( TObject * obj : *(compDir->GetListOfKeys()) )
196  {
197  TKey * key = dynamic_cast<TKey*>(obj);
198  TDirectory * subdir = dynamic_cast<TDirectory*>(key->ReadObj());
199  if (!subdir)
200  continue;
201 
202  compSpecs.emplace(key->GetName(), *(Spectrum::LoadFrom(compDir, key->GetName()).release()));
203  compDescrs.emplace(std::piecewise_construct,
204  std::forward_as_tuple(key->GetName()),
205  std::forward_as_tuple(kNoCut, // Cut object won't be re-used anyhow since Spectra are already made?
206  dynamic_cast<TObjString*>(subdir->Get("blurb"))->String().Data(),
207  dynamic_cast<TParameter<int>*>(subdir->Get("color"))->GetVal() )
208  );
209  }
210 
211  auto ret = std::make_unique<SpectrumComponents>(
212  std::move(tot),
213  std::move(compDescrs),
214  std::move(compSpecs)
215  );
216 
217  delete dir;
218 
219  return std::move(ret);
220 
221  }
222 
223 }
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
static std::unique_ptr< SpectrumComponents > LoadFrom(TDirectory *dir, const std::string &name)
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
subdir
Definition: cvnie.py:7
std::map< std::string, Component > fComponentDescr
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.
TString hists[nhists]
Definition: bdt_com.C:3
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:272
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
EBinType
Definition: UtilsExt.h:28
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
TLegend * DrawLegend(double drawThreshold=0, TLegend *leg=nullptr) const
Draw legend on plots.
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
loader
Definition: demo0.py:10
dictionary categories
Naming schema histograms should follow a logical naming convention: <category>-<observable>-<reconstr...
Definition: naming_schema.py:9
double POT() const
Definition: Spectrum.h:219
const SystShifts kNoShift
Definition: SystShifts.cxx:21
std::vector< double > POT
Base class for the various types of spectrum loader.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
SpectrumComponents(const ana::Cut &baseCut, const ana::HistAxis &axis, ana::SpectrumLoaderBase &loader, const Var &wei=kUnweighted, const std::map< std::string, Component > &categories={})
SpectrumComponents makes a collection of Spectra differing by a single (often truth) Cut...
assert(nhit_max >=nhit_nbins)
std::map< std::string, Spectrum > fComponentSpecs
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
void SaveTo(TDirectory *dir, const std::string &name) const