FOMUtilities.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "TH1.h"
4 #include "TH2.h"
5 #include "TCanvas.h"
6 #include "TLegend.h"
8 
9 enum FomType {
10  kSOverSB, // signal over sqrt background
11  kSOverSSB, // signal over sqrt (signal plus background)
12  kBinnedSOverSSB, // signal over sqrt (signal plus background)
13  kBinnedSOverSSBSq, // above, but squared
14  kEfficiency, // efficiency (total true selected / total true)
15  // -- n.b. "total true" is after the Concat cuts, so not quite a true purity
16  kPurity, // purity (total true selected / total selected)
17  kEffTimesPurity // efficiency times purity
18 };
19 
20 // Print FOM values
21 //------------------------------------------------
22 void printValues(std::vector<double> theseValues, std::vector<double> theseFOMValues){
23 
24  for (size_t i = 0; i < theseValues.size(); ++i){
25  std::cout << " " << theseValues.at(i);
26  }
28 
29  std::cout << "FOM value: ";
30  for (size_t i = 0; i < theseFOMValues.size(); ++i){
31  std::cout << " " << theseFOMValues.at(i);
32  }
33  std::cout << "\n" << std::endl;
34 
35 }
36 
37 // find cut values where FOMs are maximised
38 //------------------------------------------------
39 void getMaximumFOM(TH1D* h){
40 
41  std::vector<double> cutValues;
42  std::vector<double> FOMValues;
43 
44  double maxFOM = h->GetMaximum();
45  for (int i = 1; i < h->GetNbinsX()+1; ++i)
46  {
47  if (h->GetBinContent(i) == maxFOM)
48  {
49  cutValues.push_back(h->GetBinLowEdge(i));
50  FOMValues.push_back(h->GetBinContent(i));
51  }
52  }
53 
54  printValues(cutValues, FOMValues);
55 
56 }
57 
58 // saves some uncerainty histograms
59 //------------------------------------------------
60 void SaveUncertaintyHistogram(std::string name, std::vector<double> stat, std::vector<float> syst, std::vector<double> total){
61 
62  if (syst.size() == 0)
63  return;
64 
65  TH1D* hstat = new TH1D((name+"hstat").c_str(), ";Energy Bin;", stat.size(), 0, stat.size());
66  TH1D* hsyst = new TH1D((name+"hsyst").c_str(), ";Energy Bin;", syst.size(), 0, syst.size());
67  TH1D* htot = new TH1D((name+"htot").c_str() , ";Energy Bin;", syst.size() , 0, syst.size() );
68  TH1D* htotal = new TH1D((name+"htotal").c_str(), ";Energy Bin;", syst.size(), 0, syst.size());
69 
70  for (size_t i = 0; i < syst.size(); ++i){
71  hstat->SetBinContent(i+1, stat.at(i));
72  hsyst->SetBinContent(i+1, syst.at(i));
73  htot->SetBinContent(i+1, std::sqrt(std::pow(stat.at(i),2)+std::pow(syst.at(i),2)));
74  htotal->SetBinContent(i+1, total.at(i));
75  }
76 
77  hstat->SetLineColor(kGreen+1);
78  hstat->SetLineWidth(2);
79  hsyst->SetLineColor(kAzure+1);
80  hsyst->SetLineWidth(2);
81  htot->SetLineColor(kBlack);
82  htot->SetLineWidth(2);
83 
84  TCanvas* c3 = new TCanvas((name+"c").c_str(), "", 600, 400);
85 
86  htot ->Draw("same");
87  hsyst->Draw("same hist");
88  hstat->Draw("same hist");
89  htot ->Draw("same hist");
90 
91  c3->SaveAs(("uncertainties_"+name+".pdf").c_str());
92  c3->SaveAs(("uncertainties_"+name+".png").c_str());
93 
94  htotal->GetYaxis()->SetRangeUser(0, 60000);
95  htotal->Draw("hist");
96  c3->SaveAs(("totaluncertainties_"+name+".pdf").c_str());
97  c3->SaveAs(("totaluncertainties_"+name+".png").c_str());
98 }
99 
100 // Styles the FOM histograms and scales them
101 //------------------------------------------------
102 void StyleFOMs(TH1D* EffTimesPur = nullptr, TH1D* SOverSB = nullptr,
103  TH1D* SOverSSBSyst = nullptr, TH1D* SOverSSB = nullptr, double scalef = 1){
104 
105  if (EffTimesPur != nullptr){
106  EffTimesPur->SetLineColor(kOrange+1);
107  EffTimesPur->SetMarkerColor(kOrange+1);
108  EffTimesPur->SetLineWidth(2);
109  //std::cout << "Printing maximum FOM value for EffTimesPur" << std::endl;
110  //getMaximumFOM(EffTimesPur);
111  if (EffTimesPur->GetMaximum() != 0)
112  EffTimesPur->Scale(scalef/EffTimesPur->GetMaximum());
113  }
114 
115  if (SOverSB != nullptr){
116  SOverSB->SetLineColor(kViolet+1);
117  SOverSB->SetMarkerColor(kViolet+1);
118  SOverSB->SetLineWidth(2);
119  //std::cout << "Printing maximum FOM value for SOverSB" << std::endl;
120  //getMaximumFOM(SOverSB);
121  if (SOverSB->GetMaximum()!=0)
122  SOverSB->Scale(scalef * 0.9 / SOverSB->GetMaximum());
123  }
124 
125  if (SOverSSBSyst != nullptr){
126  SOverSSBSyst->SetLineColor(kRed+1);
127  SOverSSBSyst->SetMarkerColor(kRed+1);
128  SOverSSBSyst->SetLineWidth(2);
129  //std::cout << "Printing maximum FOM value for SOverSSBSyst" << std::endl;
130  //getMaximumFOM(SOverSSBSyst);
131  if (SOverSSBSyst->GetMaximum() != 0)
132  SOverSSBSyst->Scale(scalef * 0.8 / SOverSSBSyst->GetMaximum());
133  }
134 
135  if (SOverSSB != nullptr){
136  SOverSSB->SetLineColor(kPink+1);
137  SOverSSB->SetMarkerColor(kPink+1);
138  SOverSSB->SetLineWidth(2);
139  //std::cout << "Printing maximum FOM value for SOverSSB" << std::endl;
140  //getMaximumFOM(SOverSSB);
141  if (SOverSSB->GetMaximum() != 0){
142  SOverSSB->Scale(scalef * 0.7 / SOverSSB->GetMaximum());
143  }
144  }
145 }
146 
147 // simple function to calculate SOverSB
148 //------------------------------------------------
149 double calcSOverSB(double sig, double bg){
150  if (bg == 0) return 0;
151  else return sig/std::sqrt(bg);
152 }
153 
154 double calcSOverSSB(double sig, double bg){
155  if (sig == 0) return 10e-5;
156  return sig/sqrt(sig+bg);
157 }
158 
159 // simple function to calculate binned SOverSSB
160 //------------------------------------------------
161 double calcBinnedSOverSSB(TH1D* projSig, std::vector<double> totUnc){
162 
163  if (projSig->GetNbinsX() != (int)totUnc.size())
164  return 0;
165 
166  std::vector<double> metric = {};
167  for (int i = 0; i < projSig->GetNbinsX(); ++i){
168  metric.push_back( projSig->GetBinContent(i+1)/totUnc.at(i));
169  }
170 
171  double med = TMath::Mean(metric.size(), &metric[0]);
172 
173  return med;
174 }
175 
176 double calcBinnedSOverSSBSqr(TH1D* projSig, std::vector<double> totUnc){
177 
178  if (projSig->GetNbinsX() != (int)totUnc.size())
179  return 0;
180 
181  std::vector<double> metric = {};
182  for (int i = 0; i < projSig->GetNbinsX(); ++i){
183  metric.push_back( std::pow(projSig->GetBinContent(i+1)/totUnc.at(i),2));
184  }
185 
186  double med = TMath::Mean(metric.size(), &metric[0]);
187 
188  return med;
189 }
190 
191 // calculate efficiency
192 //------------------------------------------------
193 double calcEfficiency(double selSig, double totSig){
194  if (totSig == 0) return 0;
195  else return selSig/totSig;
196 }
197 
198 // calculate purity
199 //------------------------------------------------
200 double calcPurity(double selSig, double selTot){
201  if (selTot == 0) return 0;
202  else return selSig/selTot;
203 }
204 
205 // handler function to get correct FOM given a FOMType
206 //------------------------------------------------
207 double getFOM(FomType thisFomType, double selSig, double selBg, double totSig, double selTot, TH1D* projSig, std::vector<double> totUnc, int i = 0){
208 
209  switch (thisFomType){
210  case kSOverSB:
211  return calcSOverSB(selSig, selBg);
212  case kSOverSSB:
213  return calcSOverSSB(selSig, selBg);
214  case kBinnedSOverSSB:
215  return calcBinnedSOverSSB(projSig, totUnc);
216  case kBinnedSOverSSBSq:
217  return std::pow(calcBinnedSOverSSB(projSig, totUnc),2);
218  case kEfficiency:
219  return calcEfficiency(selSig, totSig);
220  case kPurity:
221  return calcPurity(selSig, selTot);
222  case kEffTimesPurity:
223  return calcEfficiency(selSig, totSig) * calcPurity(selSig, selTot);
224  default:
225  throw std::logic_error("thisFomType not accepted");
226  }
227 }
228 
229 // get fractional "statistical" uncertainty
230 // which is just sqrt(N)/N
231 //------------------------------------------------
232 std::vector<double> getStatUncertainty(TH1D* h){
233 
234  std::vector<double> statUnc;
235  for (int i = 0; i < h->GetNbinsX(); ++i){
236  statUnc.push_back(std::sqrt(h->GetBinContent(i+1))/h->GetBinContent(i+1));
237  }
238 
239  return statUnc;
240 
241 }
242 
243 // take fractional statistical and systematic uncertainty,
244 // add in quaderature, convert to total uncertainty
245 //------------------------------------------------
246 std::vector<double> getTotalUncertainty(TH1D* sel, std::vector<double> stat, std::vector<float> syst){
247  std::vector<double> totUncertainty;
248  if (sel == nullptr)
249  return totUncertainty;
250  for (size_t i = 0; i < stat.size(); ++i){
251  double statsq = std::pow(stat.at(i),2);
252  double systsq = std::pow(syst.at(i),2);
253  totUncertainty.push_back(std::sqrt(statsq + systsq)*sel->GetBinContent(i+1));
254  }
255  return totUncertainty;
256 }
257 
258 // function to return the histogram of the Figure of Merit
259 //------------------------------------------------
260 TH1D* getFOMHist(std::string fomName,
261  CutSide thisCutSide,
262  FomType thisFomType,
263  TH1D* all,
264  TH1D* sig,
265  TH1D* trueSig,
266  TH2D* energyCorrTotal = nullptr,
267  TH2D* energyCorrSignal = nullptr,
268  std::vector<float> systUncertainty = {}){
269 
270  TH1D* fomReturnHist = new TH1D(fomName.c_str(),
271  ";;FOM",
272  all->GetNbinsX(),
273  all->GetBinLowEdge(1),
274  all->GetBinLowEdge(all->GetNbinsX()+1));
275 
276  if ((thisFomType == FomType::kBinnedSOverSSB || thisFomType == FomType::kBinnedSOverSSBSq) &&
277  energyCorrTotal == nullptr)
278  return fomReturnHist;
279  if (thisCutSide == CutSide::kBoth ||
280  thisCutSide == CutSide::kNeither)
281  return fomReturnHist;
282 
283 
284  double selSig, totSig, selTot, selBg;
285  if (thisCutSide == CutSide::kLow){
286  for (int i = 1; i < fomReturnHist->GetNbinsX()+1; i++){
287  selTot = all->Integral(i, fomReturnHist->GetNbinsX());
288  selSig = sig->Integral(i, fomReturnHist->GetNbinsX());
289  selBg = selTot - selSig;
290  totSig = trueSig->Integral();
291 
292  // project out TH2 for this cut
293  TH1D* projTot = nullptr;
294  if (energyCorrTotal != nullptr){
295  projTot = energyCorrTotal->ProjectionX("projTot", i, fomReturnHist->GetNbinsX()+1);
296  }
297  TH1D* projSig = nullptr;
298  if (energyCorrSignal != nullptr){
299  projSig = energyCorrSignal->ProjectionX("projSig", i, fomReturnHist->GetNbinsX()+1);
300  }
301 
302  // get statistical uncertainty distribution
303  // get systematic uncertainty distribution
304  // pass those to the getFOM
305  // calcaulate average FOM across distribution
306  std::vector<double> statUncertainty;
307  if (projTot != nullptr)
308  statUncertainty = getStatUncertainty(projTot);
309 
310  std::vector<double> totalUncertainty;
311  if (systUncertainty.size() != 0)
312  totalUncertainty = getTotalUncertainty(projTot, statUncertainty, systUncertainty);
313  else totalUncertainty = statUncertainty;
314 
315  fomReturnHist->SetBinContent(i, getFOM(thisFomType, selSig, selBg, totSig, selTot, projSig, totalUncertainty, i));
316  }
317  }
318  else if (thisCutSide == CutSide::kHigh){
319  for (int i = fomReturnHist->GetNbinsX()+1; i > 0; i--){
320  selTot = all->Integral(1,i);
321  selSig = sig->Integral(1,i);
322  selBg = selTot - selSig;
323  totSig = trueSig->Integral();
324 
325  TH1D* projTot = nullptr;
326  if (energyCorrTotal != nullptr){
327  projTot = energyCorrTotal->ProjectionX("projTot", 1, i);
328  }
329  TH1D* projSig = nullptr;
330  if (energyCorrSignal != nullptr){
331  projSig = energyCorrSignal->ProjectionX("projSig", 1, i);
332  }
333 
334  // get statistical uncertainty distribution
335  // get systematic uncertainty distribution
336  // pass those to the getFOM
337  // calcaulate average FOM across distribution
338  std::vector<double> statUncertainty;
339  if (projTot != nullptr)
340  statUncertainty = getStatUncertainty(projTot);
341  std::vector<double> totalUncertainty = getTotalUncertainty(projTot, statUncertainty, systUncertainty);
342 
343  //SaveUncertaintyHistogram(std::string(all->GetName()), statUncertainty, systUncertainty);
344 
345  fomReturnHist->SetBinContent(i, getFOM(thisFomType, selSig, selBg, totSig, selTot, projSig, totalUncertainty));
346  }
347  }
348 
349  return fomReturnHist;
350 }
351 
352 
void SaveUncertaintyHistogram(std::string name, std::vector< double > stat, std::vector< float > syst, std::vector< double > total)
Definition: FOMUtilities.h:60
void getMaximumFOM(TH1D *h)
Definition: FOMUtilities.h:39
double calcPurity(double selSig, double selTot)
Definition: FOMUtilities.h:200
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
enum BeamMode kRed
void StyleFOMs(TH1D *EffTimesPur=nullptr, TH1D *SOverSB=nullptr, TH1D *SOverSSBSyst=nullptr, TH1D *SOverSSB=nullptr, double scalef=1)
Definition: FOMUtilities.h:102
T sqrt(T number)
Definition: d0nt_math.hpp:156
double getFOM(FomType thisFomType, double selSig, double selBg, double totSig, double selTot, TH1D *projSig, std::vector< double > totUnc, int i=0)
Definition: FOMUtilities.h:207
constexpr T pow(T x)
Definition: pow.h:75
std::vector< double > getStatUncertainty(TH1D *h)
Definition: FOMUtilities.h:232
void printValues(std::vector< double > theseValues, std::vector< double > theseFOMValues)
Definition: FOMUtilities.h:22
TH1D * getFOMHist(std::string fomName, CutSide thisCutSide, FomType thisFomType, TH1D *all, TH1D *sig, TH1D *trueSig, TH2D *energyCorrTotal=nullptr, TH2D *energyCorrSignal=nullptr, std::vector< float > systUncertainty={})
Definition: FOMUtilities.h:260
double calcBinnedSOverSSB(TH1D *projSig, std::vector< double > totUnc)
Definition: FOMUtilities.h:161
double calcBinnedSOverSSBSqr(TH1D *projSig, std::vector< double > totUnc)
Definition: FOMUtilities.h:176
double calcSOverSB(double sig, double bg)
Definition: FOMUtilities.h:149
OStream cout
Definition: OStream.cxx:6
double calcSOverSSB(double sig, double bg)
Definition: FOMUtilities.h:154
std::vector< double > getTotalUncertainty(TH1D *sel, std::vector< double > stat, std::vector< float > syst)
Definition: FOMUtilities.h:246
enum BeamMode kViolet
double calcEfficiency(double selSig, double totSig)
Definition: FOMUtilities.h:193
FomType
Definition: FOMUtilities.h:9
Float_t e
Definition: plot.C:35
enum BeamMode kGreen
enum BeamMode string