| 1 |
#include "binneddata.hh"
|
| 2 |
#include "fit.hh"
|
| 3 |
|
| 4 |
#include <cassert>
|
| 5 |
#include <iostream>
|
| 6 |
#include <fstream>
|
| 7 |
|
| 8 |
#include <TFile.h>
|
| 9 |
#include <TH1D.h>
|
| 10 |
|
| 11 |
TH1D* getData(const std::vector<std::string>& filenames, const char* histname, int nbins, double* bins)
|
| 12 |
{
|
| 13 |
TH1D* hist=new TH1D("normalized_data", "normalized data", nbins, bins);
|
| 14 |
|
| 15 |
unsigned nFiles = filenames.size();
|
| 16 |
|
| 17 |
TFile *files[nFiles];
|
| 18 |
TH1D *hists[nFiles];
|
| 19 |
|
| 20 |
// open input files and get input histograms
|
| 21 |
for(unsigned i=0; i<nFiles; ++i)
|
| 22 |
{
|
| 23 |
files[i] = new TFile(filenames[i].c_str());
|
| 24 |
hists[i] = (TH1D*)files[i]->Get(histname);
|
| 25 |
}
|
| 26 |
|
| 27 |
// normalize the data in the histogram
|
| 28 |
for(int i=1; i<=hist->GetNbinsX(); ++i)
|
| 29 |
{
|
| 30 |
|
| 31 |
int bin_min = hists[0]->GetXaxis()->FindBin(hist->GetXaxis()->GetBinLowEdge(i)+0.5);
|
| 32 |
int bin_max = hists[0]->GetXaxis()->FindBin(hist->GetXaxis()->GetBinUpEdge(i)-0.5);
|
| 33 |
double val=hists[0]->Integral(bin_min,bin_max);
|
| 34 |
double err=Fitter::histError(val);
|
| 35 |
double width=hist->GetBinWidth(i);
|
| 36 |
hist->SetBinContent(i, val/width);
|
| 37 |
hist->SetBinError(i, err/width);
|
| 38 |
|
| 39 |
// for simultaneous fit in 3 distribution (used in di-b-jet analysis, currently commented out)
|
| 40 |
//if(hist->GetXaxis()->GetBinUpEdge(i)<=5890.)
|
| 41 |
//{
|
| 42 |
// double offset = 0.;
|
| 43 |
// int bin_min = hists[0]->GetXaxis()->FindBin(hist->GetXaxis()->GetBinLowEdge(i)-offset+0.5);
|
| 44 |
// int bin_max = hists[0]->GetXaxis()->FindBin(hist->GetXaxis()->GetBinUpEdge(i)-offset-0.5);
|
| 45 |
// double val=hists[0]->Integral(bin_min,bin_max);
|
| 46 |
// double err=Fitter::histError(val);
|
| 47 |
// double width=hist->GetBinWidth(i);
|
| 48 |
// hist->SetBinContent(i, val/width);
|
| 49 |
// hist->SetBinError(i, err/width);
|
| 50 |
//}
|
| 51 |
//else if(hist->GetXaxis()->GetBinUpEdge(i)>5890. && hist->GetXaxis()->GetBinUpEdge(i)<=10890. && nFiles>1)
|
| 52 |
//{
|
| 53 |
// double offset = 5000.;
|
| 54 |
// int bin_min = hists[1]->GetXaxis()->FindBin(hist->GetXaxis()->GetBinLowEdge(i)-offset+0.5);
|
| 55 |
// int bin_max = hists[1]->GetXaxis()->FindBin(hist->GetXaxis()->GetBinUpEdge(i)-offset-0.5);
|
| 56 |
// double val=hists[1]->Integral(bin_min,bin_max);
|
| 57 |
// double err=Fitter::histError(val);
|
| 58 |
// double width=hist->GetBinWidth(i);
|
| 59 |
// hist->SetBinContent(i, val/width);
|
| 60 |
// hist->SetBinError(i, err/width);
|
| 61 |
//}
|
| 62 |
//else if(hist->GetXaxis()->GetBinUpEdge(i)>10890. && nFiles>2)
|
| 63 |
//{
|
| 64 |
// double offset = 10000.;
|
| 65 |
// int bin_min = hists[2]->GetXaxis()->FindBin(hist->GetXaxis()->GetBinLowEdge(i)-offset+0.5);
|
| 66 |
// int bin_max = hists[2]->GetXaxis()->FindBin(hist->GetXaxis()->GetBinUpEdge(i)-offset-0.5);
|
| 67 |
// double val=hists[2]->Integral(bin_min,bin_max);
|
| 68 |
// double err=Fitter::histError(val);
|
| 69 |
// double width=hist->GetBinWidth(i);
|
| 70 |
// hist->SetBinContent(i, val/width);
|
| 71 |
// hist->SetBinError(i, err/width);
|
| 72 |
//}
|
| 73 |
}
|
| 74 |
|
| 75 |
// do some checks
|
| 76 |
if(hist->GetBinContent(0)>0.0 || hist->GetBinContent(hist->GetNbinsX()+1)>0.0)
|
| 77 |
{
|
| 78 |
std::cerr << "There is an underflow/overflow. This should be fixed. Please rebin!" << std::endl;
|
| 79 |
assert(0);
|
| 80 |
}
|
| 81 |
if(hist->GetBinContent(1)<=0.0)
|
| 82 |
{
|
| 83 |
std::cerr << "The first data bin is empty. This will give incorrect results. Please rebin!" << std::endl;
|
| 84 |
assert(0);
|
| 85 |
}
|
| 86 |
|
| 87 |
for(unsigned i=0; i<nFiles; ++i)
|
| 88 |
{
|
| 89 |
files[i]->Close();
|
| 90 |
}
|
| 91 |
|
| 92 |
return hist;
|
| 93 |
}
|
| 94 |
|
| 95 |
|
| 96 |
TH1D* getSignalCDF(const char* filename1, const char* histname1, const char* filename2, const char* histname2, const double BR, const double eff_h, const double eff_l, const std::string& postfix)
|
| 97 |
{
|
| 98 |
TH1D* h_pdf=new TH1D("h_pdf", "Resonance Shape", 14000, 0, 14000);
|
| 99 |
TH1D* h_cdf=new TH1D(("h_cdf"+postfix).c_str(), "Resonance Shape CDF", 14000, 0, 14000);
|
| 100 |
|
| 101 |
TFile* histfile1=new TFile(filename1);
|
| 102 |
TFile* histfile2=new TFile(filename2);
|
| 103 |
|
| 104 |
TH1D* hist1=(TH1D*)histfile1->Get(histname1);
|
| 105 |
TH1D* hist2=(TH1D*)histfile2->Get(histname2);
|
| 106 |
|
| 107 |
hist1->Scale(eff_h*BR);
|
| 108 |
hist2->Scale(eff_l*(1-BR));
|
| 109 |
hist1->Add(hist2); // this adds together two resonance shapes with appropriate branching fractions and efficiencies
|
| 110 |
|
| 111 |
for(int i=1; i<=hist1->GetNbinsX(); i++)
|
| 112 |
{
|
| 113 |
|
| 114 |
int bin_min = h_pdf->GetXaxis()->FindBin(hist1->GetXaxis()->GetBinLowEdge(i)+0.5);
|
| 115 |
int bin_max = h_pdf->GetXaxis()->FindBin(hist1->GetXaxis()->GetBinUpEdge(i)-0.5);
|
| 116 |
double bin_content = hist1->GetBinContent(i)/double(bin_max-bin_min+1);
|
| 117 |
if(bin_content<0.) bin_content=0.; // some extrapolated resonance shapes can have a small negative bin content in some of the bins. this protects against such cases
|
| 118 |
for(int b=bin_min; b<=bin_max; b++){
|
| 119 |
h_pdf->SetBinContent(b, bin_content);
|
| 120 |
}
|
| 121 |
}
|
| 122 |
|
| 123 |
h_pdf->Scale(1./h_pdf->Integral()); // normalize final PDF
|
| 124 |
|
| 125 |
// produce CDF from the normalized PDF
|
| 126 |
for(int i=1; i<=h_cdf->GetNbinsX(); i++)
|
| 127 |
{
|
| 128 |
int bin_min = h_pdf->GetXaxis()->FindBin(h_cdf->GetXaxis()->GetBinLowEdge(i)+0.5);
|
| 129 |
int bin_max = h_pdf->GetXaxis()->FindBin(h_cdf->GetXaxis()->GetBinUpEdge(i)-0.5);
|
| 130 |
|
| 131 |
double curr = 0;
|
| 132 |
for(int b=bin_min; b<=bin_max; b++){
|
| 133 |
curr+=h_pdf->GetBinContent(b);
|
| 134 |
}
|
| 135 |
|
| 136 |
double prev=h_cdf->GetBinContent(i-1);
|
| 137 |
|
| 138 |
h_cdf->SetBinContent(i, prev+curr);
|
| 139 |
}
|
| 140 |
|
| 141 |
histfile1->Close();
|
| 142 |
histfile2->Close();
|
| 143 |
delete h_pdf;
|
| 144 |
|
| 145 |
return h_cdf;
|
| 146 |
}
|
