| 1 |
algomez |
1.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 |
|
|
#include <TVectorD.h>
|
| 14 |
|
|
#include <TMatrixD.h>
|
| 15 |
|
|
#include <TMatrixDSym.h>
|
| 16 |
|
|
#include <TMatrixDSymEigen.h>
|
| 17 |
|
|
|
| 18 |
|
|
#include "binneddata.hh"
|
| 19 |
|
|
#include "fit.hh"
|
| 20 |
|
|
#include "statistics.hh"
|
| 21 |
|
|
|
| 22 |
|
|
using namespace std;
|
| 23 |
|
|
|
| 24 |
|
|
//##################################################################################################################################
|
| 25 |
|
|
// User Section 1
|
| 26 |
|
|
//
|
| 27 |
|
|
// (Change the code outside the User Sections only if you know what you are doing)
|
| 28 |
|
|
//
|
| 29 |
|
|
|
| 30 |
|
|
////////////////////////////////////////////////////////////////////////////////
|
| 31 |
|
|
// magic numbers
|
| 32 |
|
|
////////////////////////////////////////////////////////////////////////////////
|
| 33 |
|
|
|
| 34 |
|
|
// number of pseudoexperiments (when greater than 0, expected limit with +/- 1 and 2 sigma bands is calculated)
|
| 35 |
|
|
const int NPES=200; // 200 (the more pseudo-experiments, the better. However, 200 is a reasonable choice)
|
| 36 |
|
|
|
| 37 |
|
|
// number of samples of nuisance parameters for Bayesian MC integration (when greater than 0, systematic uncertanties are included in the limit calculation)
|
| 38 |
|
|
const int NSAMPLES=0; // 5000 (larger value is better but it also slows down the code. 5000 is a reasonable compromise between the speed and precision)
|
| 39 |
|
|
|
| 40 |
|
|
// alpha (1-alpha=confidence interval)
|
| 41 |
|
|
const double ALPHA=0.05;
|
| 42 |
|
|
|
| 43 |
|
|
// left side tail
|
| 44 |
|
|
const double LEFTSIDETAIL=0.0;
|
| 45 |
|
|
|
| 46 |
|
|
// center-of-mass energy
|
| 47 |
|
|
const double sqrtS = 8000.;
|
| 48 |
|
|
|
| 49 |
|
|
// histogram binning
|
| 50 |
|
|
//const int NBINS=38;
|
| 51 |
|
|
/*double BOUNDARIES[NBINS+1] = { 890, 944, 1000, 1058, 1118, 1181, 1246, 1313, 1383, 1455, 1530,
|
| 52 |
|
|
1607, 1687, 1770, 1856, 1945, 2037, 2132, 2231, 2332, 2438, 2546,
|
| 53 |
|
|
2659, 2775, 2895, 3019, 3147, 3279, 3416, 3558, 3704, 3854, 4010,
|
| 54 |
|
|
4171, 4337, 4509, 4686, 4869, 5058 };*/
|
| 55 |
|
|
/* STOP 1 */
|
| 56 |
|
|
/*const int NBINS=93;
|
| 57 |
|
|
double BOUNDARIES[NBINS+1] = { 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200,
|
| 58 |
|
|
210, 220, 230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, 400,
|
| 59 |
|
|
410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600,
|
| 60 |
|
|
610, 620, 630, 640, 650, 660, 670, 680, 690, 700, 710, 720, 730, 740, 750, 760, 770, 780, 790, 800,
|
| 61 |
|
|
810, 820, 830, 840, 850, 860, 870, 880, 890, 900, 910, 920, 930, 940, 950, 960, 970, 980, 990, 1000 }; */
|
| 62 |
|
|
|
| 63 |
|
|
|
| 64 |
|
|
const int NBINS=60;
|
| 65 |
|
|
double BOUNDARIES[NBINS+1] = { 400, 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600,
|
| 66 |
|
|
610, 620, 630, 640, 650, 660, 670, 680, 690, 700, 710, 720, 730, 740, 750, 760, 770, 780, 790, 800,
|
| 67 |
|
|
810, 820, 830, 840, 850, 860, 870, 880, 890, 900, 910, 920, 930, 940, 950, 960, 970, 980, 990, 1000 };
|
| 68 |
|
|
/* 1010, 1020, 1030, 1040, 1050, 1060, 1070, 1080, 1090, 1100, 1110, 1120, 1130, 1140, 1150, 1160, 1170, 1180, 1190, 1200,
|
| 69 |
|
|
1210, 1220, 1230, 1240, 1250, 1260, 1270, 1280, 1290, 1300, 1310, 1320, 1330, 1340, 1350, 1360, 1370, 1380, 1390, 1400,
|
| 70 |
|
|
1410, 1420, 1430, 1440, 1450, 1460, 1470, 1480, 1490, 1500, 1510, 1520, 1530, 1540, 1550, 1560, 1570, 1580, 1590, 1600,
|
| 71 |
|
|
1610, 1620, 1630, 1640, 1650, 1660, 1670, 1680, 1690, 1700, 1710, 1720, 1730, 1740, 1750, 1760, 1770, 1780, 1790, 1800,
|
| 72 |
|
|
1810, 1820, 1830, 1840, 1850, 1860, 1870, 1880, 1890, 1900, 1910, 1920, 1930, 1940, 1950, 1960, 1970, 1980, 1990, 2000 }; */
|
| 73 |
|
|
|
| 74 |
|
|
// parameters
|
| 75 |
|
|
double SIGMASS=0;
|
| 76 |
|
|
const int NPARS=12;
|
| 77 |
|
|
const int NBKGPARS=4;
|
| 78 |
|
|
const int POIINDEX=0; // which parameter is "of interest" "cross section"
|
| 79 |
|
|
const char* PAR_NAMES[NPARS] = { "xs", "lumi", "jes", "jer", "p0", "p1", "p2", "p3", "n0", "n1", "n2", "n3" };
|
| 80 |
|
|
//double PAR_GUESSES[NPARS] = { 1E-6, 4976., 1.0, 1.0, 3.27713e-01, 8.33753e+00, 5.37123e+00, 4.05975e-02, 0, 0, 0, 0 };
|
| 81 |
|
|
//double PAR_GUESSES[NPARS] = { 1E-6, 20000., 1.0, 1.0, 2.18518e-04, -1.44935e-01, 6.91955e+00, 7.29084e-01, 0, 0, 0, 0 }; /// STOP1
|
| 82 |
|
|
//double PAR_GUESSES[NPARS] = { 1E-6, 20000., 1.0, 1.0, 2.77210e-06, 5.59013e+00, 8.86246e+00, 1.06444e+00, 0, 0, 0, 0 }; /// STOP1
|
| 83 |
|
|
double PAR_GUESSES[NPARS] = { 1E-6, 19500., 1.0, 1.0, 2.10567e+06, 8.16039e+01, -1.43295e+00, 2.43521e-01, 0, 0, 0, 0 };
|
| 84 |
|
|
double PAR_MIN[NPARS] = { 0, 0.0, 0.0, 0.0, -9999, -9999, -9999, -9999, -100, -100, -100, -100 };
|
| 85 |
|
|
const double PAR_MAX[NPARS] = { 1E2, 6000., 2.0, 2.0, 9999, 9999, 9999, 9999, 100, 100, 100, 100 };
|
| 86 |
|
|
double PAR_ERR[NPARS] = { 1E-6, 110., 0.0125, 0.10, 1e-02, 1e-01, 1e-01, 1e-02, 1, 1, 1, 1 };
|
| 87 |
|
|
const int PAR_TYPE[NPARS] = { 1, 2, 2, 2, 0, 0, 0, 0, 3, 3, 3, 3 }; // // 1,2 = signal (2 not used in the fit); 0,3 = background (3 not used in the fit)
|
| 88 |
|
|
//const int PAR_NUIS[NPARS] = { 0, 1, 1, 1, 0, 0, 0, 0, 4, 4, 4, 4 }; // 0 = not varied, >=1 = nuisance parameters with different priors (1 = Lognormal, 2 = Gaussian, 3 = Gamma, >=4 = Uniform)
|
| 89 |
|
|
|
| 90 |
|
|
//const int PAR_NUIS[NPARS] = { 0, 1, 1, 1, 0, 0, 0, 0, 4, 4, 4, 4 }; // all (same as above)
|
| 91 |
|
|
//const int PAR_NUIS[NPARS] = { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; // lumi only
|
| 92 |
|
|
//const int PAR_NUIS[NPARS] = { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; // jes only
|
| 93 |
|
|
//const int PAR_NUIS[NPARS] = { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }; // jer only
|
| 94 |
|
|
const int PAR_NUIS[NPARS] = { 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 4, 4 }; // background only
|
| 95 |
|
|
|
| 96 |
|
|
//
|
| 97 |
|
|
// End of User Section 1
|
| 98 |
|
|
//##################################################################################################################################
|
| 99 |
|
|
|
| 100 |
|
|
// input files vector
|
| 101 |
|
|
vector<string> INPUTFILES;
|
| 102 |
|
|
|
| 103 |
|
|
// covariance matrix
|
| 104 |
|
|
double COV_MATRIX[NPARS][NPARS];
|
| 105 |
|
|
TMatrixDSym covMatrix = TMatrixDSym(NBKGPARS);
|
| 106 |
|
|
TVectorD eigenValues = TVectorD(NBKGPARS);
|
| 107 |
|
|
TMatrixD eigenVectors = TMatrixD(NBKGPARS,NBKGPARS);
|
| 108 |
|
|
|
| 109 |
|
|
// constrain S to be positive in the S+B fit
|
| 110 |
|
|
const bool posS = 0;
|
| 111 |
|
|
|
| 112 |
|
|
// use B-only fit with fixed but non-zero signal when calculating the covariance matrix used for background systematics
|
| 113 |
|
|
const bool BonlyFitForSyst = 1;
|
| 114 |
|
|
|
| 115 |
|
|
// shift in the counter used to extract the covariance matrix
|
| 116 |
|
|
int shift = 1;
|
| 117 |
|
|
|
| 118 |
|
|
// branching fraction
|
| 119 |
|
|
double BR = 1.;
|
| 120 |
|
|
|
| 121 |
|
|
// resonance shape type
|
| 122 |
|
|
//string ResShapeType = "gg";
|
| 123 |
|
|
|
| 124 |
|
|
TH1D* HISTCDF=0; // signal CDF
|
| 125 |
|
|
|
| 126 |
|
|
|
| 127 |
|
|
////////////////////////////////////////////////////////////////////////////////
|
| 128 |
|
|
// function integral
|
| 129 |
|
|
////////////////////////////////////////////////////////////////////////////////
|
| 130 |
|
|
double INTEGRAL(double *x0, double *xf, double *par){
|
| 131 |
|
|
double xs=par[0];
|
| 132 |
|
|
double lumi=par[1];
|
| 133 |
|
|
double jes=par[2];
|
| 134 |
|
|
double jer=par[3];
|
| 135 |
|
|
double p0=par[4];
|
| 136 |
|
|
double p1=par[5];
|
| 137 |
|
|
double p2=par[6];
|
| 138 |
|
|
double p3=par[7];
|
| 139 |
|
|
double n[NBKGPARS] = {0.};
|
| 140 |
|
|
n[0]=par[8];
|
| 141 |
|
|
n[1]=par[9];
|
| 142 |
|
|
n[2]=par[10];
|
| 143 |
|
|
n[3]=par[11];
|
| 144 |
|
|
|
| 145 |
|
|
if( COV_MATRIX[0+shift][0+shift]>0. && (n[0]!=0. || n[1]!=0. || n[2]!=0. || n[3]!=0.) ) {
|
| 146 |
|
|
double g[NBKGPARS] = {0.};
|
| 147 |
|
|
for(int v=0; v<NBKGPARS; ++v) {
|
| 148 |
|
|
for(int k=0; k<NBKGPARS; ++k) g[k]=n[v]*eigenValues(v)*eigenVectors[k][v];
|
| 149 |
|
|
p0 += g[0];
|
| 150 |
|
|
p1 += g[1];
|
| 151 |
|
|
p2 += g[2];
|
| 152 |
|
|
p3 += g[3];
|
| 153 |
|
|
}
|
| 154 |
|
|
}
|
| 155 |
|
|
|
| 156 |
|
|
// uses Simpson's 3/8th rule to compute the background integral over a short interval
|
| 157 |
|
|
// also use a power series expansion to determine the intermediate intervals since the pow() call is expensive
|
| 158 |
|
|
|
| 159 |
|
|
double dx=(xf[0]-x0[0])/3./sqrtS;
|
| 160 |
|
|
double x=x0[0]/sqrtS;
|
| 161 |
|
|
double logx=log(x);
|
| 162 |
|
|
|
| 163 |
|
|
double a=pow(1-x,p1)/pow(x,p2+p3*logx);
|
| 164 |
|
|
double b=dx*a/x/(x-1)*(p2+p1*x-p2*x-2*p3*(x-1)*logx);
|
| 165 |
|
|
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 );
|
| 166 |
|
|
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 );
|
| 167 |
|
|
|
| 168 |
|
|
double bkg=(xf[0]-x0[0])*p0*(a+0.375*(b+c+d)+0.375*(2*b+4*c+8*d)+0.125*(3*b+9*c+27*d));
|
| 169 |
|
|
if(bkg<0.) bkg=1e-7;
|
| 170 |
|
|
|
| 171 |
|
|
if(xs==0.0) return bkg;
|
| 172 |
|
|
|
| 173 |
|
|
double xprimef=jes*(jer*(xf[0]-SIGMASS)+SIGMASS);
|
| 174 |
|
|
double xprime0=jes*(jer*(x0[0]-SIGMASS)+SIGMASS);
|
| 175 |
|
|
int bin1=HISTCDF->GetXaxis()->FindBin(xprimef);
|
| 176 |
|
|
int bin2=HISTCDF->GetXaxis()->FindBin(xprime0);
|
| 177 |
|
|
if(bin1<1) bin1=1;
|
| 178 |
|
|
if(bin1>HISTCDF->GetNbinsX()) bin1=HISTCDF->GetNbinsX();
|
| 179 |
|
|
if(bin2<1) bin2=1;
|
| 180 |
|
|
if(bin2>HISTCDF->GetNbinsX()) bin2=HISTCDF->GetNbinsX();
|
| 181 |
|
|
double sig=xs*lumi*(HISTCDF->GetBinContent(bin1)-HISTCDF->GetBinContent(bin2));
|
| 182 |
|
|
|
| 183 |
|
|
return bkg+sig;
|
| 184 |
|
|
}
|
| 185 |
|
|
|
| 186 |
|
|
////////////////////////////////////////////////////////////////////////////////
|
| 187 |
|
|
// main function
|
| 188 |
|
|
////////////////////////////////////////////////////////////////////////////////
|
| 189 |
|
|
|
| 190 |
|
|
int main(int argc, char* argv[]) {
|
| 191 |
|
|
|
| 192 |
|
|
if(argc<=1) {
|
| 193 |
|
|
//cout << "Usage: stats signalmass [BR ResShapeType]" << endl;
|
| 194 |
|
|
cout << "Usage: stats signalmass " << endl;
|
| 195 |
|
|
return 0;
|
| 196 |
|
|
}
|
| 197 |
|
|
|
| 198 |
|
|
string st2mass = argv[1];
|
| 199 |
|
|
SIGMASS = atof(argv[1]);
|
| 200 |
|
|
//int masspoint = int(SIGMASS);
|
| 201 |
|
|
/*if(argc>2) BR = atof(argv[2]);
|
| 202 |
|
|
if(argc>3) ResShapeType = argv[3];*/
|
| 203 |
|
|
BR = 1;
|
| 204 |
|
|
|
| 205 |
|
|
//##################################################################################################################################
|
| 206 |
|
|
// User Section 2
|
| 207 |
|
|
//
|
| 208 |
|
|
// (Change the code outside the User Sections only if you know what you are doing)
|
| 209 |
|
|
//
|
| 210 |
|
|
|
| 211 |
|
|
// input data file
|
| 212 |
|
|
//INPUTFILES.push_back("Data_and_ResonanceShapes/Final__histograms_CSVL_0Tag_WideJets.root");
|
| 213 |
|
|
//INPUTFILES.push_back("/uscms_data/d3/algomez/files/Stops/Results/data_4jet80_plots.root");
|
| 214 |
|
|
INPUTFILES.push_back("/uscms_data/d3/algomez/files/Stops/Results/data_4jet80_6jet60_plots.root");
|
| 215 |
|
|
|
| 216 |
|
|
// data histogram name
|
| 217 |
|
|
//string datahistname = "DATA__cutHisto_allPreviousCuts________DijetMass_pretag";
|
| 218 |
|
|
//string datahistname = "step2plots1D/massdijetWORecoBjetsCSVM"; /// STOP1
|
| 219 |
|
|
string datahistname = "step3plots1D/massRecoDiBjetDiJet_cutDiagStop2bbjj0";
|
| 220 |
|
|
|
| 221 |
|
|
// input signal files with resonance shapes
|
| 222 |
|
|
//string filename1 = "Data_and_ResonanceShapes/Resonance_Shapes_WideJets_" + ResShapeType + ".root";
|
| 223 |
|
|
//string filename2 = "Data_and_ResonanceShapes/Resonance_Shapes_WideJets_" + ResShapeType + ".root";
|
| 224 |
|
|
//string filename1 = "/uscms_data/d3/algomez/files/Stops/Results/st2_h_bb_st1_jj_450_200_4jet80_plots.root";
|
| 225 |
|
|
//string filename2 = "/uscms_data/d3/algomez/files/Stops/Results/st2_h_bb_st1_jj_450_300_4jet80_plots.root";
|
| 226 |
|
|
string filename1 = "/uscms_data/d3/algomez/files/Stops/Results/st2_h_bb_st1_jj_" + st2mass + "_200_4jet80_6jet60_plots.root";
|
| 227 |
|
|
//string filename2 = "/uscms_data/d3/algomez/files/Stops/Results/st2_h_bb_st1_jj_450_300_4jet80_6jet60_plots.root";
|
| 228 |
|
|
|
| 229 |
|
|
// signal histogram names
|
| 230 |
|
|
ostringstream histname1, histname2;
|
| 231 |
|
|
//histname1 << "h_" << ResShapeType << "_" << masspoint;
|
| 232 |
|
|
//histname2 << "h_" << ResShapeType << "_" << masspoint;
|
| 233 |
|
|
//histname1 << "step2plots1D/massdijetWORecoBjetsCSVM"; //// STOP1
|
| 234 |
|
|
//histname2 << "step2plots1D/massdijetWORecoBjetsCSVM"; //// STOP1
|
| 235 |
|
|
histname1 << "step3plots1D/massRecoDiBjetDiJet_cutDiagStop2bbjj0";
|
| 236 |
|
|
//histname2 << "step3plots1D/massRecoDiBjetDiJet_cutDiagStop2bbjj0";
|
| 237 |
|
|
|
| 238 |
|
|
// output file name
|
| 239 |
|
|
const string OUTPUTFILE="/uscms_data/d3/algomez/files/Stops/Limits/resultsStats_jj_" + st2mass + "_200.root";
|
| 240 |
|
|
|
| 241 |
|
|
//
|
| 242 |
|
|
// End of User Section 2
|
| 243 |
|
|
//##################################################################################################################################
|
| 244 |
|
|
|
| 245 |
|
|
if(BonlyFitForSyst) shift = 0;
|
| 246 |
|
|
|
| 247 |
|
|
if(!posS) PAR_MIN[0] = -PAR_MAX[0];
|
| 248 |
|
|
|
| 249 |
|
|
// initialize the covariance matrix
|
| 250 |
|
|
for(int i = 0; i<NPARS; ++i) { for(int j = 0; j<NPARS; ++j) COV_MATRIX[i][j]=0.; }
|
| 251 |
|
|
|
| 252 |
|
|
//HISTCDF=getSignalCDF(filename1.c_str(), histname1.str().c_str(), filename2.c_str(), histname2.str().c_str(), BR, 1., 1.);
|
| 253 |
|
|
HISTCDF=getSignalCDF(filename1.c_str(), histname1.str().c_str(), filename1.c_str(), histname1.str().c_str(), BR, 1., 1.);
|
| 254 |
|
|
|
| 255 |
|
|
//assert(HISTCDF && SIGMASS>0);
|
| 256 |
|
|
assert(HISTCDF);
|
| 257 |
|
|
|
| 258 |
|
|
// get the data
|
| 259 |
|
|
TH1D* data=getData(INPUTFILES, datahistname.c_str(), NBINS, BOUNDARIES);
|
| 260 |
|
|
|
| 261 |
|
|
// create the output file
|
| 262 |
|
|
ostringstream outputfile;
|
| 263 |
|
|
//outputfile << OUTPUTFILE.substr(0,OUTPUTFILE.find(".root")) << "_" << masspoint << "_" << BR << "_" << ResShapeType << ".root";
|
| 264 |
|
|
outputfile << OUTPUTFILE;
|
| 265 |
|
|
TFile* rootfile=new TFile(outputfile.str().c_str(), "RECREATE"); rootfile->cd();
|
| 266 |
|
|
|
| 267 |
|
|
// xs value
|
| 268 |
|
|
double XSval;
|
| 269 |
|
|
|
| 270 |
|
|
// setup an initial fitter to perform a signal+background fit
|
| 271 |
|
|
Fitter initfit(data, INTEGRAL);
|
| 272 |
|
|
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]);
|
| 273 |
|
|
|
| 274 |
|
|
// do an initial signal+background fit first
|
| 275 |
|
|
for(int i=0; i<NPARS; i++) if(PAR_TYPE[i]>=2 || PAR_MIN[i]==PAR_MAX[i]) initfit.fixParameter(i);
|
| 276 |
|
|
initfit.doFit();
|
| 277 |
|
|
XSval = initfit.getParameter(0); // get the xs value for later use
|
| 278 |
|
|
initfit.fixParameter(0); // a parameter needs to be fixed before its value can be changed
|
| 279 |
|
|
initfit.setParameter(0, 0.0); // set the xs value to 0 to get the B component of the S+B fit (for calculating pulls and generating pseudo-data)
|
| 280 |
|
|
initfit.setPrintLevel(0); // printlevel = -1 quiet (also suppresse all warnings), = 0 normal, = 1 verbose
|
| 281 |
|
|
initfit.calcPull("pull_bkg_init")->Write();
|
| 282 |
|
|
initfit.calcDiff("diff_bkg_init")->Write();
|
| 283 |
|
|
initfit.write("fit_bkg_init");
|
| 284 |
|
|
initfit.setParameter(0, XSval);
|
| 285 |
|
|
|
| 286 |
|
|
// setup the limit values
|
| 287 |
|
|
double observedLowerBound, observedUpperBound;
|
| 288 |
|
|
vector<double> expectedLowerBounds;
|
| 289 |
|
|
vector<double> expectedUpperBounds;
|
| 290 |
|
|
|
| 291 |
|
|
cout << "********************** pe=0 (data) **********************" << endl;
|
| 292 |
|
|
|
| 293 |
|
|
// setup the fitter with the input from the signal+background fit
|
| 294 |
|
|
Fitter fit_data(data, INTEGRAL);
|
| 295 |
|
|
fit_data.setPOIIndex(POIINDEX);
|
| 296 |
|
|
//fit_data.setPrintLevel(0);
|
| 297 |
|
|
for(int i=0; i<NPARS; i++) fit_data.defineParameter(i, PAR_NAMES[i], initfit.getParameter(i), PAR_ERR[i], PAR_MIN[i], PAR_MAX[i], PAR_NUIS[i]);
|
| 298 |
|
|
|
| 299 |
|
|
// perform a signal+background fit possibly followed by a background-only fit with a fixed but non-zero signal
|
| 300 |
|
|
for(int i=0; i<NPARS; i++) if(PAR_TYPE[i]>=2 || PAR_MIN[i]==PAR_MAX[i]) fit_data.fixParameter(i);
|
| 301 |
|
|
if(BonlyFitForSyst) { fit_data.doFit(); if(fit_data.getFitStatus().find("CONVERGED")==string::npos) { fit_data.fixParameter(0); fit_data.setParameter(0, 0.0); } else fit_data.fixParameter(0); }
|
| 302 |
|
|
fit_data.doFit(&COV_MATRIX[0][0], NPARS);
|
| 303 |
|
|
cout << "Data fit status: " << fit_data.getFitStatus() << endl;
|
| 304 |
|
|
fit_data.fixParameter(0); // a parameter needs to be fixed before its value can be changed
|
| 305 |
|
|
fit_data.setParameter(0, 0.0); // set the xs value to 0 to get the B component of the S+B fit (for calculating pulls and generating pseudo-data)
|
| 306 |
|
|
fit_data.setPrintLevel(0);
|
| 307 |
|
|
fit_data.calcPull("pull_bkg_0")->Write();
|
| 308 |
|
|
fit_data.calcDiff("diff_bkg_0")->Write();
|
| 309 |
|
|
fit_data.write("fit_bkg_0");
|
| 310 |
|
|
|
| 311 |
|
|
// calculate eigenvalues and eigenvectors
|
| 312 |
|
|
for(int i = 0; i<NBKGPARS; ++i) { for(int j = 0; j<NBKGPARS; ++j) { covMatrix(i,j)=COV_MATRIX[i+shift][j+shift]; } }
|
| 313 |
|
|
//covMatrix.Print();
|
| 314 |
|
|
const TMatrixDSymEigen eigen_data(covMatrix);
|
| 315 |
|
|
eigenValues = eigen_data.GetEigenValues();
|
| 316 |
|
|
eigenValues.Sqrt();
|
| 317 |
|
|
//eigenValues.Print();
|
| 318 |
|
|
eigenVectors = eigen_data.GetEigenVectors();
|
| 319 |
|
|
//eigenVectors.Print();
|
| 320 |
|
|
|
| 321 |
|
|
fit_data.setParLimits(0, 0.0, PAR_MAX[0]); // for the posterior calculation, the signal xs has to be positive
|
| 322 |
|
|
TGraph* post_data=fit_data.calculatePosterior(NSAMPLES);
|
| 323 |
|
|
post_data->Write("post_0");
|
| 324 |
|
|
cout << "Call limit reached: " << (fit_data.callLimitReached() ? "True" : "False") << endl;
|
| 325 |
|
|
|
| 326 |
|
|
// evaluate the limit
|
| 327 |
|
|
pair<double, double> bounds_data=evaluateInterval(post_data, ALPHA, LEFTSIDETAIL);
|
| 328 |
|
|
observedLowerBound=bounds_data.first;
|
| 329 |
|
|
observedUpperBound=bounds_data.second;
|
| 330 |
|
|
|
| 331 |
|
|
// reset the covariance matrix
|
| 332 |
|
|
for(int i = 0; i<NPARS; ++i) { for(int j = 0; j<NPARS; ++j) COV_MATRIX[i][j]=0.; }
|
| 333 |
|
|
|
| 334 |
|
|
// perform the PEs
|
| 335 |
|
|
for(int pe=1; pe<=NPES; ++pe) {
|
| 336 |
|
|
|
| 337 |
|
|
cout << "********************** pe=" << pe << " **********************" << endl;
|
| 338 |
|
|
ostringstream pestr;
|
| 339 |
|
|
pestr << "_" << pe;
|
| 340 |
|
|
|
| 341 |
|
|
// setup the fitter with the input from the signal+background fit
|
| 342 |
|
|
fit_data.setParameter(0, 0.0); // set the xs value to 0 to get the B component of the S+B fit (for calculating pulls and generating pseudo-data)
|
| 343 |
|
|
TH1D* hist = fit_data.makePseudoData((string("data")+pestr.str()).c_str());
|
| 344 |
|
|
fit_data.setParameter(0, PAR_GUESSES[0]);
|
| 345 |
|
|
|
| 346 |
|
|
Fitter fit(hist, INTEGRAL);
|
| 347 |
|
|
fit.setPOIIndex(POIINDEX);
|
| 348 |
|
|
fit.setPrintLevel(0);
|
| 349 |
|
|
for(int i=0; i<NPARS; i++) fit.defineParameter(i, PAR_NAMES[i], fit_data.getParameter(i), PAR_ERR[i], PAR_MIN[i], PAR_MAX[i], PAR_NUIS[i]);
|
| 350 |
|
|
|
| 351 |
|
|
// perform a signal+background fit possibly followed by a background-only fit with a fixed but non-zero signal
|
| 352 |
|
|
for(int i=0; i<NPARS; i++) if(PAR_TYPE[i]>=2 || PAR_MIN[i]==PAR_MAX[i]) fit.fixParameter(i);
|
| 353 |
|
|
if(BonlyFitForSyst) { fit.doFit(); if(fit.getFitStatus().find("CONVERGED")==string::npos) { fit.fixParameter(0); fit.setParameter(0, 0.0); } else fit.fixParameter(0); }
|
| 354 |
|
|
fit.doFit(&COV_MATRIX[0][0], NPARS);
|
| 355 |
|
|
if(fit.getFitStatus().find("CONVERGED")==string::npos) continue; // skip the PE if the fit did not converge
|
| 356 |
|
|
fit.fixParameter(0); // a parameter needs to be fixed before its value can be changed
|
| 357 |
|
|
fit.setParameter(0, 0.0); // set the xs value to 0 to get the B component of the S+B fit (for calculating pulls and generating pseudo-data)
|
| 358 |
|
|
fit.calcPull((string("pull_bkg")+pestr.str()).c_str())->Write();
|
| 359 |
|
|
fit.calcDiff((string("diff_bkg")+pestr.str()).c_str())->Write();
|
| 360 |
|
|
fit.write((string("fit_bkg")+pestr.str()).c_str());
|
| 361 |
|
|
|
| 362 |
|
|
// calculate eigenvalues and eigenvectors
|
| 363 |
|
|
for(int i = 0; i<NBKGPARS; ++i) { for(int j = 0; j<NBKGPARS; ++j) { covMatrix(i,j)=COV_MATRIX[i+shift][j+shift]; } }
|
| 364 |
|
|
const TMatrixDSymEigen eigen(covMatrix);
|
| 365 |
|
|
eigenValues = eigen.GetEigenValues();
|
| 366 |
|
|
bool hasNegativeElement = false;
|
| 367 |
|
|
for(int i = 0; i<NBKGPARS; ++i) { if(eigenValues(i)<0.) hasNegativeElement = true; }
|
| 368 |
|
|
if(hasNegativeElement) continue; // this is principle should never happen. However, if it does, skip the PE
|
| 369 |
|
|
eigenValues.Sqrt();
|
| 370 |
|
|
eigenVectors = eigen.GetEigenVectors();
|
| 371 |
|
|
|
| 372 |
|
|
fit.setParLimits(0, 0.0, PAR_MAX[0]); // for the posterior calculation, the signal xs has to be positive
|
| 373 |
|
|
TGraph* post=fit.calculatePosterior(NSAMPLES);
|
| 374 |
|
|
post->Write((string("post")+pestr.str()).c_str());
|
| 375 |
|
|
cout << "Call limit reached in pe=" << pe << ": " << (fit.callLimitReached() ? "True" : "False") << endl;
|
| 376 |
|
|
|
| 377 |
|
|
// evaluate the limit
|
| 378 |
|
|
pair<double, double> bounds=evaluateInterval(post, ALPHA, LEFTSIDETAIL);
|
| 379 |
|
|
if(bounds.first==0. && bounds.second>0.) {
|
| 380 |
|
|
expectedLowerBounds.push_back(bounds.first);
|
| 381 |
|
|
expectedUpperBounds.push_back(bounds.second);
|
| 382 |
|
|
}
|
| 383 |
|
|
|
| 384 |
|
|
// reset the covariance matrix
|
| 385 |
|
|
for(int i = 0; i<NPARS; ++i) { for(int j = 0; j<NPARS; ++j) COV_MATRIX[i][j]=0.; }
|
| 386 |
|
|
}
|
| 387 |
|
|
|
| 388 |
|
|
////////////////////////////////////////////////////////////////////////////////
|
| 389 |
|
|
// print the results
|
| 390 |
|
|
////////////////////////////////////////////////////////////////////////////////
|
| 391 |
|
|
|
| 392 |
|
|
cout << "**********************************************************************" << endl;
|
| 393 |
|
|
for(unsigned int i=0; i<expectedLowerBounds.size(); i++)
|
| 394 |
|
|
cout << "expected bound(" << (i+1) << ") = [ " << expectedLowerBounds[i] << " , " << expectedUpperBounds[i] << " ]" << endl;
|
| 395 |
|
|
|
| 396 |
|
|
cout << "\nobserved bound = [ " << observedLowerBound << " , " << observedUpperBound << " ]" << endl;
|
| 397 |
|
|
|
| 398 |
|
|
if(LEFTSIDETAIL>0.0 && NPES>0) {
|
| 399 |
|
|
cout << "\n***** expected lower bounds *****" << endl;
|
| 400 |
|
|
double median;
|
| 401 |
|
|
pair<double, double> onesigma;
|
| 402 |
|
|
pair<double, double> twosigma;
|
| 403 |
|
|
getQuantiles(expectedLowerBounds, median, onesigma, twosigma);
|
| 404 |
|
|
cout << "median: " << median << endl;
|
| 405 |
|
|
cout << "+/-1 sigma band: [ " << onesigma.first << " , " << onesigma.second << " ] " << endl;
|
| 406 |
|
|
cout << "+/-2 sigma band: [ " << twosigma.first << " , " << twosigma.second << " ] " << endl;
|
| 407 |
|
|
}
|
| 408 |
|
|
|
| 409 |
|
|
if(LEFTSIDETAIL<1.0 && NPES>0) {
|
| 410 |
|
|
cout << "\n***** expected upper bounds *****" << endl;
|
| 411 |
|
|
double median;
|
| 412 |
|
|
pair<double, double> onesigma;
|
| 413 |
|
|
pair<double, double> twosigma;
|
| 414 |
|
|
getQuantiles(expectedUpperBounds, median, onesigma, twosigma);
|
| 415 |
|
|
cout << "median: " << median << endl;
|
| 416 |
|
|
cout << "+/-1 sigma band: [ " << onesigma.first << " , " << onesigma.second << " ] " << endl;
|
| 417 |
|
|
cout << "+/-2 sigma band: [ " << twosigma.first << " , " << twosigma.second << " ] " << endl;
|
| 418 |
|
|
}
|
| 419 |
|
|
|
| 420 |
|
|
// close the output file
|
| 421 |
|
|
rootfile->Close();
|
| 422 |
|
|
|
| 423 |
|
|
return 0;
|
| 424 |
|
|
}
|
