| 1 |
#include <iostream>
|
| 2 |
#include <cmath>
|
| 3 |
#include <cassert>
|
| 4 |
#include <sstream>
|
| 5 |
|
| 6 |
#include <TGraph.h>
|
| 7 |
#include <TFile.h>
|
| 8 |
#include <TH1D.h>
|
| 9 |
#include <TCanvas.h>
|
| 10 |
#include <TF1.h>
|
| 11 |
#include <TMath.h>
|
| 12 |
#include <TROOT.h>
|
| 13 |
|
| 14 |
#include "binneddata.hh"
|
| 15 |
#include "fit.hh"
|
| 16 |
#include "statistics.hh"
|
| 17 |
|
| 18 |
////////////////////////////////////////////////////////////////////////////////
|
| 19 |
// magic numbers
|
| 20 |
////////////////////////////////////////////////////////////////////////////////
|
| 21 |
|
| 22 |
// number of pseudoexperiments
|
| 23 |
const int NPES=0;
|
| 24 |
|
| 25 |
// number of samples of nuisance parameters for Bayesian MC integration
|
| 26 |
const int NSAMPLES=2E5;
|
| 27 |
|
| 28 |
// alpha (1-alpha=confidence interval)
|
| 29 |
const double ALPHA=0.05;
|
| 30 |
|
| 31 |
// left side tail
|
| 32 |
const double LEFTSIDETAIL=0.0;
|
| 33 |
|
| 34 |
// output file name
|
| 35 |
const char* OUTPUTFILE="stats.root";
|
| 36 |
|
| 37 |
// input file name
|
| 38 |
const char* INPUTFILE="dijet_mass_HT_fat_1p010fbm1.txt";
|
| 39 |
|
| 40 |
// histogram binning
|
| 41 |
const int NBINS=34;
|
| 42 |
double BOUNDARIES[NBINS] = { 838,
|
| 43 |
890, 944, 1000, 1058, 1118, 1181, 1246, 1313, 1383,
|
| 44 |
1455, 1530, 1607, 1687, 1770, 1856, 1945, 2037, 2132,
|
| 45 |
2231, 2332, 2438, 2546, 2659, 2775, 2895, 3019, 3147,
|
| 46 |
3275, 3403, 3531, 3659, 3787, 3915 };
|
| 47 |
|
| 48 |
// parameters
|
| 49 |
double SIGMASS=0;
|
| 50 |
const int NPARS=8;
|
| 51 |
const int POIINDEX=0; // which parameter is "of interest"
|
| 52 |
const char* PAR_NAMES[8] = { "xs", "lumi", "jes", "jer", "bkg norm", "p1", "p2", "p3" };
|
| 53 |
const double PAR_GUESSES[8] = { 0.1, 1010., 1.0, 1.0, 5.6E-2, 7.4, 6.3, 0.2 };
|
| 54 |
const double PAR_MIN[8] = { 0.0, 0.0, 0.0, 0.0, 0., 0.0, 0.0, 0.0 };
|
| 55 |
const double PAR_MAX[8] = { 1.E6, 5000., 2.0, 2.0, 10.0, 100.0, 100.0, 10.0 };
|
| 56 |
const double PAR_ERR[8] = { 0.01, 40.4, 0.04, 0.10, 1.E-3, 1.0, 1.0, 0.1 };
|
| 57 |
const int PAR_TYPE[8] = { 1, 1, 1, 1, 0, 0, 0, 0 }; // 1 = signal, 0 = background
|
| 58 |
const int PAR_NUIS[8] = { 0, 0, 0, 0, 0, 0, 0, 0 }; // 1 = nuisance parameter, 0 = not varied (the POI is not a nuisance parameter)
|
| 59 |
|
| 60 |
TH1D* HIST=0; // signal histogram
|
| 61 |
TH1D* HISTCDF=0; // signal CDF
|
| 62 |
|
| 63 |
////////////////////////////////////////////////////////////////////////////////
|
| 64 |
// function definition (not really used; only the integral is used)
|
| 65 |
////////////////////////////////////////////////////////////////////////////////
|
| 66 |
double FCN(double *x, double *par)
|
| 67 |
{
|
| 68 |
double invmass=x[0];
|
| 69 |
double xs=par[0];
|
| 70 |
double lumi=par[1];
|
| 71 |
double jes=par[2];
|
| 72 |
double jer=par[3];
|
| 73 |
double norm=par[4];
|
| 74 |
double p1=par[5];
|
| 75 |
double p2=par[6];
|
| 76 |
double p3=par[7];
|
| 77 |
|
| 78 |
double bkg = norm*pow(1.0-invmass/7000.0,p1)/pow(invmass/7000.0,p2+p3*log(invmass/7000.0));
|
| 79 |
double mass=jes*(jer*(invmass-SIGMASS)+SIGMASS);
|
| 80 |
int bin=HIST->GetXaxis()->FindBin(mass);
|
| 81 |
double sig = xs*lumi*HIST->GetBinContent(bin);
|
| 82 |
return bkg+sig;
|
| 83 |
}
|
| 84 |
|
| 85 |
////////////////////////////////////////////////////////////////////////////////
|
| 86 |
// function integral
|
| 87 |
////////////////////////////////////////////////////////////////////////////////
|
| 88 |
double INTEGRAL(double *x0, double *xf, double *par)
|
| 89 |
{
|
| 90 |
double xs=par[0];
|
| 91 |
double lumi=par[1];
|
| 92 |
double jes=par[2];
|
| 93 |
double jer=par[3];
|
| 94 |
double norm=par[4];
|
| 95 |
double p1=par[5];
|
| 96 |
double p2=par[6];
|
| 97 |
double p3=par[7];
|
| 98 |
|
| 99 |
// uses Simpson's 3/8th rule to compute the background integral over a short interval
|
| 100 |
// also use a power series expansion to determine the intermediate intervals since the pow() call is expensive
|
| 101 |
|
| 102 |
double dx=(xf[0]-x0[0])/3./7000.;
|
| 103 |
double x=x0[0]/7000.0;
|
| 104 |
double logx=log(x);
|
| 105 |
|
| 106 |
double a=pow(1-x,p1)/pow(x,p2+p3*logx);
|
| 107 |
double b=dx*a/x/(x-1)*(p2+p1*x-p2*x-2*p3*(x-1)*logx);
|
| 108 |
double c=0.5*dx*dx*a*( (p1-1)*p1/(x-1)/(x-1) - 2*p1*(p2+2*p3*logx)/(x-1)/x + (p2+p2*p2-2*p3+2*p3*logx*(1+2*p2+2*p3*logx))/x/x );
|
| 109 |
double d=0.166666667*dx*dx*dx*a*( (p1-2)*(p1-1)*p1/(x-1)/(x-1)/(x-1) - 3*(p1-1)*p1*(p2+2*p3*logx)/(x-1)/(x-1)/x - (1+p2+2*p3*logx)*(p2*(2+p2) - 6*p3 + 4*p3*logx*(1+p2*p3*logx))/x/x/x + 3*p1*(p2+p2*p2-2*p3+2*p3*logx*(1+2*p2+2*p3*logx))/(x-1)/x/x );
|
| 110 |
|
| 111 |
double bkg=(xf[0]-x0[0])*norm*(a+0.375*(b+c+d)+0.375*(2*b+4*c+8*d)+0.125*(3*b+9*c+27*d));
|
| 112 |
|
| 113 |
if(xs==0.0) return bkg;
|
| 114 |
|
| 115 |
double xprimef=jes*(jer*(xf[0]-SIGMASS)+SIGMASS);
|
| 116 |
double xprime0=jes*(jer*(x0[0]-SIGMASS)+SIGMASS);
|
| 117 |
int bin1=HISTCDF->GetXaxis()->FindBin(xprimef);
|
| 118 |
int bin2=HISTCDF->GetXaxis()->FindBin(xprime0);
|
| 119 |
if(bin1<1) bin1=1;
|
| 120 |
if(bin1>HISTCDF->GetNbinsX()) bin1=HISTCDF->GetNbinsX();
|
| 121 |
if(bin2<1) bin1=1;
|
| 122 |
if(bin2>HISTCDF->GetNbinsX()) bin2=HISTCDF->GetNbinsX();
|
| 123 |
double sig=xs*lumi*(HISTCDF->GetBinContent(bin1)-HISTCDF->GetBinContent(bin2));
|
| 124 |
|
| 125 |
return bkg+sig;
|
| 126 |
}
|
| 127 |
|
| 128 |
////////////////////////////////////////////////////////////////////////////////
|
| 129 |
// main function
|
| 130 |
////////////////////////////////////////////////////////////////////////////////
|
| 131 |
|
| 132 |
int main(int argc, char* argv[])
|
| 133 |
{
|
| 134 |
if(argc<=1) {
|
| 135 |
std::cout << "Usage: stats signalmass" << std::endl;
|
| 136 |
return 0;
|
| 137 |
}
|
| 138 |
|
| 139 |
// setup the signal histogram
|
| 140 |
TFile* histfile=new TFile("Test_Resonance_Shapes.root");
|
| 141 |
histfile->cd();
|
| 142 |
SIGMASS = std::atof(argv[1]);
|
| 143 |
int masspoint = static_cast<int>(SIGMASS);
|
| 144 |
std::ostringstream histname, cdfname;
|
| 145 |
histname << "h_qstar_" << masspoint;
|
| 146 |
cdfname << "h_qstar_" << masspoint << "_cdf";
|
| 147 |
HIST=dynamic_cast<TH1D*>(gROOT->FindObject(histname.str().c_str()));
|
| 148 |
HISTCDF=dynamic_cast<TH1D*>(gROOT->FindObject(cdfname.str().c_str()));
|
| 149 |
assert(HIST && HISTCDF && SIGMASS>0);
|
| 150 |
HIST->Scale(5.); // proper normalization
|
| 151 |
|
| 152 |
// create the output file
|
| 153 |
TFile* rootfile=new TFile(OUTPUTFILE, "RECREATE"); rootfile->cd();
|
| 154 |
|
| 155 |
// get the data
|
| 156 |
TH1D* data=getData(INPUTFILE, "data_0", NBINS-1, BOUNDARIES);
|
| 157 |
|
| 158 |
// setup an initial fitter to perform a background-only fit
|
| 159 |
Fitter initfit(data, INTEGRAL);
|
| 160 |
for(int i=0; i<NPARS; i++) initfit.defineParameter(i, PAR_NAMES[i], PAR_GUESSES[i], PAR_ERR[i], PAR_MIN[i], PAR_MAX[i], PAR_NUIS[i]);
|
| 161 |
|
| 162 |
// do an initial background-only fit, first
|
| 163 |
for(int i=0; i<NPARS; i++) if(PAR_TYPE[i]==1) initfit.fixParameter(i);
|
| 164 |
initfit.setParameter(POIINDEX, 0.0); // set the POI value to 0
|
| 165 |
initfit.doFit();
|
| 166 |
initfit.calcPull("pull_bkg_init")->Write();
|
| 167 |
initfit.calcDiff("diff_bkg_init")->Write();
|
| 168 |
initfit.write("fit_bkg_init");
|
| 169 |
|
| 170 |
// setup the limit values
|
| 171 |
double observedLowerBound, observedUpperBound;
|
| 172 |
std::vector<double> expectedLowerBounds;
|
| 173 |
std::vector<double> expectedUpperBounds;
|
| 174 |
|
| 175 |
// perform the PEs (0 = data)
|
| 176 |
for(int pe=0; pe<=NPES; ++pe) {
|
| 177 |
|
| 178 |
std::cout << "*********** pe=" << pe << " ***********" << std::endl;
|
| 179 |
std::ostringstream pestr;
|
| 180 |
pestr << "_" << pe;
|
| 181 |
|
| 182 |
// setup the fitter with the input from the background-only fit
|
| 183 |
TH1D* hist = (pe==0) ? data : initfit.makePseudoData((std::string("data")+pestr.str()).c_str());
|
| 184 |
Fitter fit(hist, INTEGRAL);
|
| 185 |
fit.setPOIIndex(POIINDEX);
|
| 186 |
fit.setPrintLevel(0);
|
| 187 |
for(int i=0; i<NPARS; i++) fit.defineParameter(i, PAR_NAMES[i], initfit.getParameter(i), PAR_ERR[i], PAR_MIN[i], PAR_MAX[i], PAR_NUIS[i]);
|
| 188 |
|
| 189 |
// perform a background-only fit
|
| 190 |
for(int i=0; i<NPARS; i++) if(PAR_TYPE[i]==1) fit.fixParameter(i);
|
| 191 |
fit.setParameter(POIINDEX, 0.0); // set the POI value to 0
|
| 192 |
fit.doFit();
|
| 193 |
fit.calcPull((std::string("pull_bkg")+pestr.str()).c_str())->Write();
|
| 194 |
fit.calcDiff((std::string("diff_bkg")+pestr.str()).c_str())->Write();
|
| 195 |
fit.write((std::string("fit_bkg")+pestr.str()).c_str());
|
| 196 |
|
| 197 |
// fix the ranges for the background parameters before calculating the posterior
|
| 198 |
for(int i=0; i<NPARS; i++) {
|
| 199 |
if(PAR_TYPE[i]==0 && PAR_NUIS[i]==1) {
|
| 200 |
double val, err;
|
| 201 |
fit.getParameter(i, val, err);
|
| 202 |
fit.setParLimits(i, val-err, val+err);
|
| 203 |
}
|
| 204 |
}
|
| 205 |
|
| 206 |
observedLowerBound=0.;
|
| 207 |
observedUpperBound=fit.calculateUpperBoundWithCLs(NSAMPLES, ALPHA);
|
| 208 |
|
| 209 |
|
| 210 |
// put the ranges back in place
|
| 211 |
for(int i=0; i<NPARS; i++) {
|
| 212 |
if(PAR_TYPE[i]==0 && PAR_NUIS[i]==1) {
|
| 213 |
fit.setParLimits(i, PAR_MIN[i], PAR_MAX[i]);
|
| 214 |
}
|
| 215 |
}
|
| 216 |
|
| 217 |
// evaluate the limit
|
| 218 |
/* std::pair<double, double> bounds=evaluateInterval(post, ALPHA, LEFTSIDETAIL);
|
| 219 |
if(pe==0) {
|
| 220 |
observedLowerBound=bounds.first;
|
| 221 |
observedUpperBound=bounds.second;
|
| 222 |
} else {
|
| 223 |
expectedLowerBounds.push_back(bounds.first);
|
| 224 |
expectedUpperBounds.push_back(bounds.second);
|
| 225 |
}*/
|
| 226 |
}
|
| 227 |
|
| 228 |
////////////////////////////////////////////////////////////////////////////////
|
| 229 |
// print the results
|
| 230 |
////////////////////////////////////////////////////////////////////////////////
|
| 231 |
|
| 232 |
std::cout << "**********************************************************************" << std::endl;
|
| 233 |
for(unsigned int i=0; i<expectedLowerBounds.size(); i++)
|
| 234 |
std::cout << "expected bound(" << (i+1) << ") = [ " << expectedLowerBounds[i] << " , " << expectedUpperBounds[i] << " ]" << std::endl;
|
| 235 |
|
| 236 |
std::cout << "\nobserved bound = [ " << observedLowerBound << " , " << observedUpperBound << " ]" << std::endl;
|
| 237 |
|
| 238 |
if(LEFTSIDETAIL>0.0 && NPES>0) {
|
| 239 |
std::cout << "\n***** expected lower bounds *****" << std::endl;
|
| 240 |
double median;
|
| 241 |
std::pair<double, double> onesigma;
|
| 242 |
std::pair<double, double> twosigma;
|
| 243 |
getQuantiles(expectedLowerBounds, median, onesigma, twosigma);
|
| 244 |
std::cout << "median: " << median << std::endl;
|
| 245 |
std::cout << "+/-1 sigma band: [ " << onesigma.first << " , " << onesigma.second << " ] " << std::endl;
|
| 246 |
std::cout << "+/-2 sigma band: [ " << twosigma.first << " , " << twosigma.second << " ] " << std::endl;
|
| 247 |
}
|
| 248 |
|
| 249 |
if(LEFTSIDETAIL<1.0 && NPES>0) {
|
| 250 |
std::cout << "\n***** expected upper bounds *****" << std::endl;
|
| 251 |
double median;
|
| 252 |
std::pair<double, double> onesigma;
|
| 253 |
std::pair<double, double> twosigma;
|
| 254 |
getQuantiles(expectedUpperBounds, median, onesigma, twosigma);
|
| 255 |
std::cout << "median: " << median << std::endl;
|
| 256 |
std::cout << "+/-1 sigma band: [ " << onesigma.first << " , " << onesigma.second << " ] " << std::endl;
|
| 257 |
std::cout << "+/-2 sigma band: [ " << twosigma.first << " , " << twosigma.second << " ] " << std::endl;
|
| 258 |
}
|
| 259 |
|
| 260 |
// close the output file
|
| 261 |
rootfile->Close();
|
| 262 |
|
| 263 |
return 0;
|
| 264 |
}
|
