| 1 |
#include "../cl95cms/cl95cms.C"
|
| 2 |
#include "cuts.C"
|
| 3 |
#include "TFile.h"
|
| 4 |
#include "TH1.h"
|
| 5 |
#include "TCanvas.h"
|
| 6 |
#include "TGraphErrors.h"
|
| 7 |
#include <map>
|
| 8 |
#include <iomanip>
|
| 9 |
#include <utility>
|
| 10 |
|
| 11 |
#define PRODUCT(a,b,ea,eb) a*b,sqrt(ea*b*ea*b+a*eb*a*eb)
|
| 12 |
//const double lumi = 34.18, lumiErr = 0.052;
|
| 13 |
const double lumi = 34.74, lumiErr = 0.11;
|
| 14 |
|
| 15 |
const double bgExpect[7] = {20., 7.3, 3.0, 1.3, 0.57, 0.26, 0.11};
|
| 16 |
const double bgUncert[7] = {10.1, 5.2, 2.6, 1.25, 0.58, 0.26, 0.12};
|
| 17 |
const int bgReal [7] = {16, 9, 5, 1, 1, 1, 1};
|
| 18 |
|
| 19 |
// Old backgrounds:
|
| 20 |
//const double bgExpect[7] = {29., 12.0, 5.3, 2.4, 1.2, 0.54, 0.24};
|
| 21 |
//const double bgUncert[7] = {9., 4.7, 2.7, 1.4, 0.7, 0.40, 0.18};
|
| 22 |
|
| 23 |
using namespace std;
|
| 24 |
|
| 25 |
map<string, map<int, pair<double,double> > > gpCuts;
|
| 26 |
|
| 27 |
double lim(double eff, double effErr, double _bgExpect, double _bgExpectErr, int bg){
|
| 28 |
double actualLimitSigma = CL95(lumi, lumi*lumiErr, eff, effErr, _bgExpect, _bgExpectErr, bg);
|
| 29 |
cout<<"Actual limit on sigma="<<actualLimitSigma<<" (pb)"<<endl;
|
| 30 |
return actualLimitSigma;
|
| 31 |
}
|
| 32 |
|
| 33 |
double run(double eff, double effErr, double _bgExpect, double _bgExpectErr){
|
| 34 |
double expectedLimitSigma = CLA (lumi, lumi*lumiErr, eff, effErr, _bgExpect, _bgExpectErr);
|
| 35 |
cout<<"Expected limit on sigma="<<expectedLimitSigma<<" (pb)"<<endl;
|
| 36 |
return expectedLimitSigma;
|
| 37 |
}
|
| 38 |
|
| 39 |
void calcEff(const char *signalNtupleName=0){
|
| 40 |
if( !signalNtupleName ){
|
| 41 |
cout<<"calcEff: give me ntuple"<<endl;
|
| 42 |
return;
|
| 43 |
}
|
| 44 |
|
| 45 |
cout<<"Sample: "<<signalNtupleName<<endl;
|
| 46 |
|
| 47 |
double limits[7];
|
| 48 |
|
| 49 |
// Run the threshold
|
| 50 |
for(int threshold=200,thrNumber=0; threshold<550; threshold+=50,thrNumber++){
|
| 51 |
cout<<"Calculating limit for the threshold "<<threshold<<endl;
|
| 52 |
cout<<"Expect "<<bgExpect[thrNumber]<<" +- "<<bgUncert[thrNumber]<<" events of background"<<endl;
|
| 53 |
|
| 54 |
double eff = gpCuts[signalNtupleName][threshold].first; // qT threshold efficiency
|
| 55 |
double effErr = gpCuts[signalNtupleName][threshold].second; // and the absolute statistical uncertainty
|
| 56 |
|
| 57 |
// double eff = cuts(signalNtupleName,threshold); // qT threshold efficiency
|
| 58 |
// double effErr = cuts_error; // and the absolute statistical uncertainty
|
| 59 |
cout<<"Using eff="<<eff<<" and effErr="<<effErr<<endl;
|
| 60 |
limits[thrNumber] = run(eff,effErr,bgExpect[thrNumber],bgUncert[thrNumber]);
|
| 61 |
cout<<endl<<endl;
|
| 62 |
}
|
| 63 |
|
| 64 |
int min=0;
|
| 65 |
double minLimit=100;
|
| 66 |
for(int thrNumber=0; thrNumber<7; thrNumber++)
|
| 67 |
if( limits[thrNumber]<minLimit ){
|
| 68 |
minLimit = limits[thrNumber];
|
| 69 |
min = thrNumber;
|
| 70 |
}
|
| 71 |
|
| 72 |
cout<<"For "<<signalNtupleName<<" the best limit is "<<minLimit<<" is at pT threshold: "<<(200+min*50)<<endl;
|
| 73 |
}
|
| 74 |
|
| 75 |
void limitOptimization(void){
|
| 76 |
map<int,pair<double,double> > CI05;
|
| 77 |
/* Old eff:
|
| 78 |
CI05[200] = pair<double,double>(0.59976*0.5364,0.0076); //2501/4170
|
| 79 |
CI05[250] = pair<double,double>(0.59424*0.4351,0.0083); //2084/3507
|
| 80 |
CI05[300] = pair<double,double>(0.58871*0.3466,0.0092); //1689/2869
|
| 81 |
CI05[350] = pair<double,double>(0.56217*0.2796,0.0113); //1076/1914
|
| 82 |
CI05[400] = pair<double,double>(0.53359*0.2215,0.0127); // 826/1548
|
| 83 |
CI05[450] = pair<double,double>(0.53359*0.1724,0.0127); // 826/1548
|
| 84 |
CI05[500] = pair<double,double>(0.50552*0.1315,0.0140); // 641/1268
|
| 85 |
*/
|
| 86 |
CI05[200] = pair<double,double>(PRODUCT(0.747,0.496,0.008,0.007));
|
| 87 |
CI05[250] = pair<double,double>(PRODUCT(0.740,0.409,0.008,0.007));
|
| 88 |
CI05[300] = pair<double,double>(PRODUCT(0.730,0.331,0.009,0.007));
|
| 89 |
CI05[350] = pair<double,double>(PRODUCT(0.737,0.269,0.010,0.006));
|
| 90 |
CI05[400] = pair<double,double>(PRODUCT(0.717,0.213,0.012,0.006));
|
| 91 |
CI05[450] = pair<double,double>(PRODUCT(0.690,0.165,0.013,0.005));
|
| 92 |
CI05[500] = pair<double,double>(PRODUCT(0.649,0.125,0.015,0.005));
|
| 93 |
gpCuts["../cl95cms/tuples/old/tup_ci_M05.root"] = CI05;
|
| 94 |
|
| 95 |
// calcEff("../cl95cms/tuples/old/tup_ci_M05.root");
|
| 96 |
|
| 97 |
cout<<"Expect CI05 limit "<<run(CI05[350].first,CI05[350].second,bgExpect[3],bgUncert[3])<<endl;
|
| 98 |
cout<<"Actual CI05 limit "<<lim(CI05[350].first,CI05[350].second,bgExpect[3],bgUncert[3],bgReal[3])<<endl;
|
| 99 |
// cout<<" (old CI05 limit "<<lim(0.204147,0.007,bgExpect[2],bgUncert[2],bgReal[2])<<")"<<endl;
|
| 100 |
// cout<<" (old CI10 limit "<<lim(0.204147,0.007,3,2.6,5)<<")"<<endl;
|
| 101 |
|
| 102 |
map<int,pair<double,double> > CI10;
|
| 103 |
/* Old eff:
|
| 104 |
CI10[200] = pair<double,double>(0.567604*0.702642,0.00647); //3329/5865
|
| 105 |
CI10[250] = pair<double,double>(0.573392*0.633208,0.00670); //3129/5457
|
| 106 |
CI10[300] = pair<double,double>(0.575979*0.556075,0.00704); //2839/4929
|
| 107 |
CI10[350] = pair<double,double>(0.582302*0.478943,0.00747); //2540/4362
|
| 108 |
CI10[400] = pair<double,double>(0.574853*0.398943,0.00807); //2158/3754
|
| 109 |
CI10[450] = pair<double,double>(0.561973*0.330113,0.00877); //1800/3203
|
| 110 |
CI10[500] = pair<double,double>(0.542021*0.263245,0.00959); //1464/2701
|
| 111 |
*/
|
| 112 |
CI10[200] = pair<double,double>(PRODUCT(0.734,0.664,0.007,0.007));
|
| 113 |
CI10[250] = pair<double,double>(PRODUCT(0.737,0.621,0.007,0.007));
|
| 114 |
CI10[300] = pair<double,double>(PRODUCT(0.735,0.562,0.007,0.007));
|
| 115 |
CI10[350] = pair<double,double>(PRODUCT(0.726,0.495,0.008,0.007));
|
| 116 |
CI10[400] = pair<double,double>(PRODUCT(0.712,0.419,0.008,0.007));
|
| 117 |
CI10[450] = pair<double,double>(PRODUCT(0.700,0.350,0.009,0.007));
|
| 118 |
CI10[500] = pair<double,double>(PRODUCT(0.664,0.280,0.010,0.006));
|
| 119 |
gpCuts["../cl95cms/tuples/old/tup_ci_M10.root"] = CI10;
|
| 120 |
|
| 121 |
// calcEff("../cl95cms/tuples/old/tup_ci_M10.root");
|
| 122 |
cout<<"Expect CI10 limit "<<run(CI10[400].first,CI10[400].second,bgExpect[4],bgUncert[4])<<endl;
|
| 123 |
cout<<"Actual CI10 limit "<<lim(CI10[400].first,CI10[400].second,bgExpect[4],bgUncert[4],bgReal[4])<<endl;
|
| 124 |
|
| 125 |
map<int,pair<double,double> > CI15;
|
| 126 |
/*
|
| 127 |
CI15[200] = pair<double,double>(0.512116*0.71772,0.0064); //3149/6149
|
| 128 |
CI15[250] = pair<double,double>(0.512016*0.68105,0.0065); //3068/5992
|
| 129 |
CI15[300] = pair<double,double>(0.514795*0.64022,0.0066); //2975/5779
|
| 130 |
CI15[350] = pair<double,double>(0.513435*0.59322,0.0067); //2828/5508
|
| 131 |
CI15[400] = pair<double,double>(0.511898*0.54499,0.0069); //2689/5253
|
| 132 |
CI15[450] = pair<double,double>(0.508294*0.48413,0.0072); //2482/4883
|
| 133 |
CI15[500] = pair<double,double>(0.497284*0.424961,0.0075); //2197/4418
|
| 134 |
*/
|
| 135 |
CI15[200] = pair<double,double>(PRODUCT(0.668,0.640,0.007,0.007));
|
| 136 |
CI15[250] = pair<double,double>(PRODUCT(0.668,0.623,0.007,0.007));
|
| 137 |
CI15[300] = pair<double,double>(PRODUCT(0.668,0.604,0.007,0.007));
|
| 138 |
CI15[350] = pair<double,double>(PRODUCT(0.663,0.573,0.007,0.007));
|
| 139 |
CI15[400] = pair<double,double>(PRODUCT(0.651,0.538,0.007,0.007));
|
| 140 |
CI15[450] = pair<double,double>(PRODUCT(0.634,0.490,0.008,0.007));
|
| 141 |
CI15[500] = pair<double,double>(PRODUCT(0.627,0.439,0.008,0.007));
|
| 142 |
|
| 143 |
gpCuts["../cl95cms/tuples/old/tup_ci_M15.root"] = CI15;
|
| 144 |
|
| 145 |
// calcEff("../cl95cms/tuples/old/tup_ci_M15.root");
|
| 146 |
cout<<"Expect CI15 limit "<<run(CI15[400].first,CI15[400].second,bgExpect[4],bgUncert[4])<<endl;
|
| 147 |
cout<<"Actual CI15 limit "<<lim(CI15[400].first,CI15[400].second,bgExpect[4],bgUncert[4],bgReal[4])<<endl;
|
| 148 |
|
| 149 |
map<int,pair<double,double> > CI20;
|
| 150 |
/*
|
| 151 |
CI20[200] = pair<double,double>(0.4452*0.699049,0.0062); //2880/6469
|
| 152 |
CI20[250] = pair<double,double>(0.4467*0.670083,0.0062); //2843/6365
|
| 153 |
CI20[300] = pair<double,double>(0.4463*0.642602,0.0063); //2798/6270
|
| 154 |
CI20[350] = pair<double,double>(0.4448*0.614825,0.0063); //2730/6137
|
| 155 |
CI20[400] = pair<double,double>(0.4434*0.581699,0.0064); //2667/6015
|
| 156 |
CI20[450] = pair<double,double>(0.4381*0.548425,0.0065); //2551/5823
|
| 157 |
CI20[500] = pair<double,double>(0.4335*0.511586,0.0066); //2436/5620
|
| 158 |
*/
|
| 159 |
CI20[200] = pair<double,double>(PRODUCT(0.592,0.576,0.007,0.007));
|
| 160 |
CI20[250] = pair<double,double>(PRODUCT(0.593,0.568,0.007,0.007));
|
| 161 |
CI20[300] = pair<double,double>(PRODUCT(0.592,0.558,0.007,0.007));
|
| 162 |
CI20[350] = pair<double,double>(PRODUCT(0.588,0.544,0.007,0.007));
|
| 163 |
CI20[400] = pair<double,double>(PRODUCT(0.583,0.529,0.007,0.007));
|
| 164 |
CI20[450] = pair<double,double>(PRODUCT(0.575,0.506,0.007,0.007));
|
| 165 |
CI20[500] = pair<double,double>(PRODUCT(0.567,0.482,0.008,0.007));
|
| 166 |
|
| 167 |
gpCuts["../cl95cms/tuples/old/tup_ci_M20.root"] = CI20;
|
| 168 |
|
| 169 |
// calcEff("../cl95cms/tuples/old/tup_ci_M20.root");
|
| 170 |
cout<<"Expect CI20 limit "<<run(CI20[400].first,CI20[400].second,bgExpect[4],bgUncert[4])<<endl;
|
| 171 |
cout<<"Actual CI20 limit "<<lim(CI20[400].first,CI20[400].second,bgExpect[4],bgUncert[4],bgReal[4])<<endl;
|
| 172 |
|
| 173 |
map<int,pair<double,double> > GI05;
|
| 174 |
/*
|
| 175 |
GI05[200] = pair<double,double>(0.7042*0.434767,0.0086); //1995/2833
|
| 176 |
GI05[250] = pair<double,double>(0.6761*0.128997,0.0167); // 528/781
|
| 177 |
GI05[300] = pair<double,double>(0.6613*0.046490,0.0269); // 205/310
|
| 178 |
GI05[350] = pair<double,double>(0.6609*0.024807,0.0359); // 115/174
|
| 179 |
GI05[400] = pair<double,double>(0.7079*0.014149,0.0482); // 63/89
|
| 180 |
GI05[450] = pair<double,double>(0.775*0.0088203, 0.066); // 31/40
|
| 181 |
GI05[500] = pair<double,double>(0.7391*0.004778,0.0916); // 17/23
|
| 182 |
*/
|
| 183 |
GI05[200] = pair<double,double>(PRODUCT(0.735,0.3852,0.009,0.007));
|
| 184 |
GI05[250] = pair<double,double>(PRODUCT(0.803,0.1102,0.015,0.004));
|
| 185 |
GI05[300] = pair<double,double>(PRODUCT(0.785,0.0402,0.026,0.003));
|
| 186 |
GI05[350] = pair<double,double>(PRODUCT(0.784,0.0218,0.035,0.002));
|
| 187 |
GI05[400] = pair<double,double>(PRODUCT(0.822,0.0120,0.045,0.002));
|
| 188 |
GI05[450] = pair<double,double>(PRODUCT(0.973,0.0072,0.027,0.001));
|
| 189 |
GI05[500] = pair<double,double>(PRODUCT(0.818,0.0036,0.082,0.001));
|
| 190 |
|
| 191 |
gpCuts["../cl95cms/tuples/old/tup_qgf_M05.root"] = GI05;
|
| 192 |
|
| 193 |
// calcEff("../cl95cms/tuples/old/tup_qgf_M05.root");
|
| 194 |
cout<<"Expect GI05 limit "<<run(GI05[200].first,GI05[200].second,bgExpect[0],bgUncert[0])<<endl;
|
| 195 |
// cout<<" (old GI05 limit "<<run(.306075968,0.009,bgExpect[0],bgUncert[0])<<")"<<endl;
|
| 196 |
cout<<"Actual GI05 limit "<<lim(GI05[200].first,GI05[200].second,bgExpect[0],bgUncert[0],bgReal[0])<<endl;
|
| 197 |
// cout<<" (old GI05 limit "<<lim(.306075968,0.009,bgExpect[0],bgUncert[0],bgReal[0])<<")"<<endl;
|
| 198 |
|
| 199 |
map<int,pair<double,double> > GI10;
|
| 200 |
/*
|
| 201 |
GI10[200] = pair<double,double>(0.696193*0.792589,0.0063); //3694/5306
|
| 202 |
GI10[250] = pair<double,double>(0.698324*0.729971,0.0065); //3500/5012
|
| 203 |
GI10[300] = pair<double,double>(0.70368 *0.654486,0.0067); //3232/4593
|
| 204 |
GI10[350] = pair<double,double>(0.701044*0.552582,0.0072); //2821/4024
|
| 205 |
GI10[400] = pair<double,double>(0.703646*0.429233,0.0080); //2277/3236
|
| 206 |
GI10[450] = pair<double,double>(0.698768*0.277063,0.0098); //1531/2191
|
| 207 |
GI10[500] = pair<double,double>(0.664491*0.122319,0.0171); // 509/766
|
| 208 |
*/
|
| 209 |
GI10[200] = pair<double,double>(PRODUCT(0.803,0.735,0.006,0.006));
|
| 210 |
GI10[250] = pair<double,double>(PRODUCT(0.807,0.695,0.006,0.007));
|
| 211 |
GI10[300] = pair<double,double>(PRODUCT(0.808,0.637,0.006,0.007));
|
| 212 |
GI10[350] = pair<double,double>(PRODUCT(0.799,0.548,0.007,0.007));
|
| 213 |
GI10[400] = pair<double,double>(PRODUCT(0.781,0.433,0.008,0.007));
|
| 214 |
GI10[450] = pair<double,double>(PRODUCT(0.749,0.280,0.010,0.006));
|
| 215 |
GI10[500] = pair<double,double>(PRODUCT(0.956,0.122,0.008,0.005));
|
| 216 |
|
| 217 |
gpCuts["../cl95cms/tuples/old/tup_qgf_M10.root"] = GI10;
|
| 218 |
|
| 219 |
// calcEff("../cl95cms/tuples/old/tup_qgf_M10.root");
|
| 220 |
cout<<"Expect GI10 limit "<<run(GI10[350].first,GI10[350].second,bgExpect[3],bgUncert[3])<<endl;
|
| 221 |
cout<<"Actual GI10 limit "<<lim(GI10[350].first,GI10[350].second,bgExpect[3],bgUncert[3],bgReal[3])<<endl;
|
| 222 |
|
| 223 |
map<int,pair<double,double> > GI15;
|
| 224 |
/*
|
| 225 |
GI15[200] = pair<double,double>(0.70088*0.768832,0.0067); //3264/4657
|
| 226 |
GI15[250] = pair<double,double>(0.70383*0.729915,0.0072); //2809/3991
|
| 227 |
GI15[300] = pair<double,double>(0.71012*0.69534, 0.0083); //2141/3015
|
| 228 |
GI15[350] = pair<double,double>(0.72122*0.654585,0.0112); //1159/1607
|
| 229 |
GI15[400] = pair<double,double>(0.72072*0.607149,0.0213); // 320/444
|
| 230 |
GI15[450] = pair<double,double>(0.68692*0.556539,0.0317); // 147/214
|
| 231 |
GI15[500] = pair<double,double>(0.64394*0.490062,0.0417); // 85/132
|
| 232 |
*/
|
| 233 |
GI15[200] = pair<double,double>(PRODUCT(0.774,0.690,0.006,0.007));
|
| 234 |
GI15[250] = pair<double,double>(PRODUCT(0.770,0.669,0.007,0.007));
|
| 235 |
GI15[300] = pair<double,double>(PRODUCT(0.759,0.649,0.008,0.007));
|
| 236 |
GI15[350] = pair<double,double>(PRODUCT(0.763,0.619,0.011,0.007));
|
| 237 |
GI15[400] = pair<double,double>(PRODUCT(0.956,0.585,0.010,0.007));
|
| 238 |
GI15[450] = pair<double,double>(PRODUCT(0.918,0.547,0.020,0.007));
|
| 239 |
GI15[500] = pair<double,double>(PRODUCT(0.847,0.485,0.034,0.007));
|
| 240 |
|
| 241 |
gpCuts["../cl95cms/tuples/old/tup_qgf_M15.root"] = GI15;
|
| 242 |
|
| 243 |
// calcEff("../cl95cms/tuples/old/tup_qgf_M15.root");
|
| 244 |
cout<<"Expect GI15 limit "<<run(GI15[400].first,GI15[400].second,bgExpect[4],bgUncert[4])<<endl;
|
| 245 |
cout<<"Actual GI15 limit "<<lim(GI15[400].first,GI15[400].second,bgExpect[4],bgUncert[4],bgReal[4])<<endl;
|
| 246 |
|
| 247 |
// for(map<int,pair<double,double> >::const_iterator iter=GI15.begin();
|
| 248 |
// iter!=GI15.end(); iter++) //cout<<"pT>"<<iter->first<<" "<<iter->second.first<<" +- "<<iter->second.second<<endl;
|
| 249 |
//cout<<setw(4)<<iter->second.first*100<<"$\\pm$"<<setw(3)<<iter->second.second*100<<"\\%"<<endl;
|
| 250 |
}
|
