-
Notifications
You must be signed in to change notification settings - Fork 4.3k
/
SideBandSubtraction.cc
494 lines (456 loc) · 16.6 KB
/
SideBandSubtraction.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
///////////////////////////////////////////////////////////////////////
// Author: David Bjergaard
//
// Library for Generic Sideband Subtraction Methods.
//
// This library is designed to be a tool which serves two purposes.
// Its primary purpose is to provide a generic framework for doing
// sideband subtraction. It will also plug into current tag and probe
// modules to prevent code duplication and redundancy.
//
///////////////////////////////////////////////////////////////////////
#include "PhysicsTools/Utilities/interface/SideBandSubtraction.h"
// System includes
#include <iostream>
#include <ios>
#include <fstream>
#include <cstdlib>
#include <sstream>
// ROOT includes
#include <TCanvas.h>
#include <TFile.h>
#include <TF1.h>
#include <TH1F.h>
#include <TString.h>
#include <TKey.h>
#include <TClass.h>
// RooFit includes
#include <RooFitResult.h>
#include <RooRealVar.h>
#include <RooAbsPdf.h>
#include <RooDataSet.h>
#include <RooPlot.h>
using namespace RooFit;
using std::cerr;
using std::cout;
using std::endl;
using std::string;
using std::vector;
template <class T>
inline std::string stringify(const T& t) {
std::ostringstream o;
if (!(o << t))
return "err";
return o.str();
}
Double_t SideBandSubtract::getYield(const std::vector<SbsRegion>& Regions, RooAbsPdf* PDF) {
if (PDF == nullptr || SeparationVariable == nullptr)
return 0.0;
Double_t yield = 0;
RooAbsReal* intPDF;
for (unsigned int i = 0; i < Regions.size(); i++) {
if (verbose)
cout << "Integrating over Range: " << Regions[i].RegionName << " from " << Regions[i].min << " to "
<< Regions[i].max << endl;
intPDF = PDF->createIntegral(*SeparationVariable, Range(Regions[i].RegionName.c_str()));
yield += intPDF->getVal();
if (verbose)
cout << "Current yield: " << yield << endl;
}
return yield;
}
static void setHistOptions(TH1F* histo, string name, string title, string axis_label) {
histo->SetName(name.c_str());
histo->SetTitle(title.c_str());
if (axis_label == "GeV/c^2")
axis_label = "Mass (" + axis_label + ")";
if (axis_label == "GeV/c")
axis_label = "Momentum (" + axis_label + ")";
histo->GetXaxis()->SetTitle(axis_label.c_str());
}
int SideBandSubtract::doSubtraction(RooRealVar* variable,
Double_t stsratio,
Int_t index) //stsratio -> signal to sideband ratio
{
if (Data == nullptr || SeparationVariable == nullptr) {
cerr << "ERROR: Data or SeparationVariable is NULL returning now!\n";
return -1;
}
TH1F* SideBandHist = (TH1F*)BaseHistos[index]->Clone();
setHistOptions(SideBandHist,
(string)variable->GetName() + "Sideband",
(string)SideBandHist->GetTitle() + " Sideband",
(string)variable->getUnit());
TH1F* SignalHist = (TH1F*)BaseHistos[index]->Clone();
setHistOptions(SignalHist,
(string)variable->GetName() + "SignalHist",
(string)SignalHist->GetTitle() + " Raw Signal",
(string)variable->getUnit());
//Begin a loop over the data to fill our histograms. I should figure
//out how to do this in one shot to avoid a loop
//O(N_vars*N_events)...
RooRealVar* sep_var = nullptr;
for (const auto& var : *Data->get()) {
if ((string)var->GetName() == (string)SeparationVariable->GetName()) {
sep_var = (RooRealVar*)var;
break;
}
}
for (int i = 0; i < Data->numEntries(); i++) {
Data->get(i);
Double_t value = variable->getVal();
Double_t cutval = sep_var->getVal();
for (unsigned int j = 0; j < SideBandRegions.size(); j++) //UGLY! Find better solution!
{
if (cutval > SideBandRegions[j].min && cutval < SideBandRegions[j].max)
SideBandHist->Fill(value);
}
for (unsigned int j = 0; j < SignalRegions.size(); j++) {
if (cutval > SignalRegions[j].min && cutval < SignalRegions[j].max)
SignalHist->Fill(value);
}
}
//Save pre-subtracted histo
SignalHist->Sumw2();
SideBandHist->Sumw2();
//SignalHist->SetDirectory(0); SideBandHist->SetDirectory(0);
RawHistos.push_back(*SignalHist);
SignalHist->Add(SideBandHist, -stsratio);
SignalHist->SetTitle(((string)BaseHistos[index]->GetTitle() + " SBS Signal").c_str());
//Save subtracted histo
SBSHistos.push_back(*SignalHist);
//Save Sideband histo
SideBandHistos.push_back(*SideBandHist);
if (SideBandHist)
delete SideBandHist;
return 0;
}
static void print_histo(TH1F* plot, string outname) {
TCanvas genericCanvas;
plot->Draw("E1P0");
outname = outname + ".eps";
genericCanvas.SaveAs(outname.c_str());
outname.replace(outname.size() - 3, 3, "gif");
genericCanvas.SaveAs(outname.c_str());
}
void SideBandSubtract::printResults(string prefix) { //handles *all* printing
//spool over vectors of histograms and print them, then print
//separation variable plots and the results text file.
if (SeparationVariable == nullptr) {
cerr << "ERROR: printResults, SeparationVariable is NULL!\n";
return;
}
string filename; //output file name
for (unsigned int i = 0; i < RawHistos.size(); ++i) {
filename = prefix + "Raw_" + (string)RawHistos[i].GetName();
print_histo(&RawHistos[i], filename);
}
for (unsigned int i = 0; i < SBSHistos.size(); ++i) {
filename = prefix + "SBS_" + (string)RawHistos[i].GetName();
print_histo(&SBSHistos[i], filename);
}
if (verbose) {
for (unsigned int i = 0; i < SideBandHistos.size(); ++i) {
filename = prefix + "SideBand_" + (string)RawHistos[i].GetName();
print_histo(&SideBandHistos[i], filename);
}
}
string outname = prefix + (string)SeparationVariable->GetName() + "_fitted.eps";
if (Data != nullptr && ModelPDF != nullptr) {
RooPlot* SepVarFrame = SeparationVariable->frame();
Data->plotOn(SepVarFrame);
ModelPDF->plotOn(SepVarFrame);
TCanvas Canvas;
SepVarFrame->Draw();
Canvas.SaveAs(outname.c_str());
outname.replace(outname.size() - 3, 3, "gif");
Canvas.SaveAs(outname.c_str());
} else
cerr << "ERROR: printResults, Data or ModelPDF is NULL!\n";
string result_outname = prefix + "_fit_results.txt";
std::ofstream output(result_outname.c_str(), std::ios::out);
if (!output) {
cout << "ERROR: Could not open file for writing!\n";
return;
}
if (fit_result != nullptr) {
#if ROOT_VERSION_CODE <= ROOT_VERSION(5, 19, 0)
fit_result->printToStream(output, Verbose);
#else
fit_result->printMultiline(output, kTRUE);
#endif
}
}
void SideBandSubtract::saveResults(string outname) {
//saves the ephemeral stuff to a root file for future
//use/modification (ie everything printed by printResults())
TFile output(outname.c_str(), "UPDATE"); //open the output file,
//create it if it doesn't
//exist
//Since keys are only available from files on disk, we need to write
//out a new file. If the file already existed, then we opened to
//update, and are writing nothing new.
output.Write();
TString dirname;
TIter nextkey(output.GetListOfKeys());
TKey* key;
TDirectory* curDir = nullptr;
while ((key = (TKey*)nextkey.Next())) {
if (key == nullptr)
break;
TObject* obj = key->ReadObj();
if (obj->IsA()->InheritsFrom("TDirectory"))
dirname = obj->GetName();
}
if (dirname == "") {
//we didn't find any directories so, we'll make a new one
curDir = output.mkdir("run0", "Run 0");
output.cd("run0");
} else {
//manipulate last dir string, make a new one, and get ready to fill
dirname.Remove(0, 3);
Int_t run_num = dirname.Atoi();
run_num++;
dirname = "run" + stringify(run_num);
curDir = output.mkdir(dirname.Data(), ("Run " + stringify(run_num)).c_str());
output.cd(dirname.Data());
}
if (curDir == nullptr)
curDir = output.GetDirectory("", kTRUE, "");
//these should all be the same size, but to be pedantic we'll loop
//over each one individually, also, we need to associate them with
//the directory because by default they "float" in memory to avoid
//conflicts with other root files the user has open. If they want to
//write to those files, they need to close their file, pass the name
//here, and then let us work.
for (unsigned int i = 0; i < SideBandHistos.size(); ++i) {
SideBandHistos[i].SetDirectory(curDir);
SideBandHistos[i].Write();
}
for (unsigned int i = 0; i < RawHistos.size(); ++i) {
RawHistos[i].SetDirectory(curDir);
RawHistos[i].Write();
}
for (unsigned int i = 0; i < SBSHistos.size(); ++i) {
SBSHistos[i].SetDirectory(curDir);
SBSHistos[i].Write();
}
if (Data != nullptr && ModelPDF != nullptr && BackgroundPDF != nullptr && SeparationVariable != nullptr) {
RooPlot* sep_varFrame = SeparationVariable->frame();
Data->plotOn(sep_varFrame);
ModelPDF->plotOn(sep_varFrame);
BackgroundPDF->plotOn(sep_varFrame);
sep_varFrame->Write();
} else
cerr << "ERROR: saveResults, did not save RooPlot of data and fit\n";
output.Write();
}
void SideBandSubtract::setDataSet(RooDataSet* newData) {
if (newData != nullptr)
Data = newData;
}
void SideBandSubtract::print_plot(RooRealVar* printVar, string outname) {
if (Data == nullptr || ModelPDF == nullptr) {
cerr << "ERROR: print_plot, Data or ModelPDF are NULL\n";
return;
}
RooPlot* genericFrame = printVar->frame();
Data->plotOn(genericFrame);
ModelPDF->plotOn(genericFrame);
TCanvas genericCanvas;
genericFrame->Draw();
outname = outname + ".eps";
genericCanvas.SaveAs(outname.c_str());
outname.replace(outname.size() - 3, 3, "gif");
genericCanvas.SaveAs(outname.c_str());
}
SideBandSubtract::SideBandSubtract()
: BackgroundPDF(nullptr),
ModelPDF(nullptr),
Data(nullptr),
SeparationVariable(nullptr),
verbose(false),
SignalRegions(),
SideBandRegions(),
SideBandHistos(0),
RawHistos(0),
SBSHistos(0),
BaseHistos(0),
fit_result(nullptr),
SignalSidebandRatio(0) {
// Default constructor so we can do fast subtraction
}
SideBandSubtract::SideBandSubtract(RooAbsPdf* model_shape,
RooAbsPdf* bkg_shape,
RooDataSet* data,
RooRealVar* sep_var,
const vector<TH1F*>& base,
bool verb)
: BackgroundPDF(bkg_shape),
ModelPDF(model_shape),
Data(data),
SeparationVariable(sep_var),
verbose(verb),
SignalRegions(),
SideBandRegions(),
SideBandHistos(),
RawHistos(),
SBSHistos(),
BaseHistos(base),
base_histo(nullptr),
fit_result(nullptr),
SignalSidebandRatio(0) {}
/*
SideBandSubtract::SideBandSubtract(RooAbsPdf *model_shape,
RooAbsPdf *bkg_shape,
RooDataSet* data,
RooRealVar* sep_var,
bool verb
)
: BackgroundPDF(bkg_shape),
ModelPDF(model_shape),
Data(data),
SeparationVariable(sep_var),
verbose(verb),
SignalRegions(),
SideBandRegions(),
SideBandHistos(),
RawHistos(),
SBSHistos(),
BaseHistos(),
base_histo(0),
fit_result(0),
SignalSidebandRatio(0)
{
// Build BaseHistos from dataset
TIterator* iter = (TIterator*) Data->get()->createIterator();
RooAbsArg *var=NULL;
RooRealVar *variable=NULL;
cout <<"Starting while loop\n";
while((var = (RooAbsArg*)iter->Next()))
{
if((string)var->ClassName() != "RooRealVar")
continue;
variable = (RooRealVar*)var;
//we own the data this points to so we need to delete it at the end...
assert(variable!=NULL);
string title = "base_"+(string)variable->GetName();
base_histo = (TH1F*)Data->createHistogram(title.c_str(), *variable, Binning(variable->getBinning("default",verb,kTRUE)) );
//cout <<"Made histo with name: "<<base_histo->GetName()<<endl;
//base_histo->SetDirectory(0);
BaseHistos.push_back(*base_histo);
cout <<"Added histo to BaseHistos!\n Deleting local copy...";
if(base_histo) delete base_histo;
cout <<"Done!\n";
}
}
*/
SideBandSubtract::~SideBandSubtract() {
//**WARNING**
// We don't delete objects that we don't own (duh) so, all of our
// pointers just hang out and get handled by other people :)
// We DO own the BaseHistos IFF the user didn't provide us with
// them and we constructed them from the dataset... in this case
// base_histo will be non-null and we want to delete the vector
// of histograms...
resetSBSProducts();
/*
if(base_histo!=NULL)
{
BaseHistos.clear();
}
*/
}
void SideBandSubtract::addSignalRegion(Double_t min, Double_t max) {
SbsRegion signal;
signal.min = min;
signal.max = max;
signal.RegionName = "Signal" + stringify(SignalRegions.size());
if (SeparationVariable != nullptr)
SeparationVariable->setRange(signal.RegionName.c_str(), signal.min, signal.max);
SignalRegions.push_back(signal);
return;
}
void SideBandSubtract::addSideBandRegion(Double_t min, Double_t max) {
SbsRegion sideband;
sideband.min = min;
sideband.max = max;
sideband.RegionName = "SideBand" + stringify(SideBandRegions.size());
if (SeparationVariable != nullptr)
SeparationVariable->setRange(sideband.RegionName.c_str(), sideband.min, sideband.max);
SideBandRegions.push_back(sideband);
return;
}
int SideBandSubtract::doGlobalFit() {
if (verbose)
cout << "Beginning SideBand Subtraction\n";
if (ModelPDF != nullptr && Data != nullptr && SeparationVariable != nullptr) {
fit_result = ModelPDF->fitTo(*Data);
} else {
cerr << "ERROR: doGobalFit, no ModelPDF, SeparationVariable or Data specified\n";
return -1;
}
Double_t SideBandYield = getYield(SideBandRegions, BackgroundPDF);
Double_t BackgroundInSignal = getYield(SignalRegions, BackgroundPDF);
SignalSidebandRatio = BackgroundInSignal / SideBandYield;
if (verbose) {
cout << "Finished fitting background!\n";
cout << "Attained a Signal to Sideband ratio of: " << SignalSidebandRatio << endl;
}
//I can't see a way around a double loop for each variable. Maybe I
//can get a C++/RooFit guru to save me the trouble here?
//need to grab sbs objects after each global fit, because they get reset
resetSBSProducts();
for (const auto& variable : *Data->get()) {
for (unsigned int i = 0; i < BaseHistos.size(); i++) {
if ((string)variable->GetName() != (string)SeparationVariable->GetName() &&
(string)variable->GetName() == (string)BaseHistos[i]->GetName())
doSubtraction((RooRealVar*)variable, SignalSidebandRatio, i);
}
}
return 0;
}
void SideBandSubtract::doFastSubtraction(TH1F& Total, TH1F& Result, SbsRegion& leftRegion, SbsRegion& rightRegion) {
Int_t binMin = Total.FindBin(leftRegion.max, 0.0, 0.0);
Int_t binMax = Total.FindBin(leftRegion.min, 0.0, 0.0);
double numLeft = Total.Integral(binMin, binMax);
binMin = Total.FindBin(rightRegion.max, 0.0, 0.0);
binMax = Total.FindBin(rightRegion.min, 0.0, 0.0);
double numRight = Total.Integral(binMin, binMax);
const unsigned int nbinsx = Total.GetNbinsX();
const double x1 = (leftRegion.max + leftRegion.min) / 2.0;
const double x2 = (rightRegion.max + rightRegion.min) / 2.0;
const double y1 = numLeft / (leftRegion.max - leftRegion.min);
const double y2 = numRight / (rightRegion.max - rightRegion.min);
const double Slope = (y2 - y1) / (x2 - x1);
const double Intercept = y1 - Slope * x1;
// Evantually we want to handle more complicated functions, but for
// now... just use a linear one
TF1 function("sbs_function_line", "[0]*x + [1]", Total.GetMinimum(), Total.GetMaximum());
for (unsigned int binx = 1; binx <= nbinsx; ++binx) {
double binWidth = Total.GetBinWidth(binx);
function.SetParameter(0, Slope * binWidth);
function.SetParameter(1, Intercept * binWidth);
double xx = Total.GetBinCenter(binx);
double cu = Total.GetBinContent(binx) - function.Eval(xx);
// TODO: Propogate the error on the parameters in function.
double error1 = Total.GetBinError(binx);
Result.SetBinContent(binx, cu);
Result.SetBinError(binx, error1);
}
Result.SetEntries(Result.Integral());
}
RooFitResult* SideBandSubtract::getFitResult() { return fit_result; }
vector<TH1F*> SideBandSubtract::getBaseHistos() { return BaseHistos; }
vector<TH1F> SideBandSubtract::getRawHistos() { return RawHistos; }
vector<TH1F> SideBandSubtract::getSBSHistos() { return SBSHistos; }
Double_t SideBandSubtract::getSTSRatio() { return SignalSidebandRatio; }
void SideBandSubtract::resetSBSProducts() {
//cout <<"Starting to reset products \n";
if (!SideBandHistos.empty())
SideBandHistos.clear();
if (!RawHistos.empty())
RawHistos.clear();
if (!SBSHistos.empty())
SBSHistos.clear();
}