| 1 |
//system
|
| 2 |
#include <iostream>
|
| 3 |
#include <cmath>
|
| 4 |
//root
|
| 5 |
#include "TCanvas.h"
|
| 6 |
#include "TH1F.h"
|
| 7 |
#include "TF1.h"
|
| 8 |
#include "TLegend.h"
|
| 9 |
#include "TStyle.h"
|
| 10 |
#include "TRandom3.h"
|
| 11 |
#include "TGraph.h"
|
| 12 |
#include "TLatex.h"
|
| 13 |
#include "TLine.h"
|
| 14 |
#include "TLimit.h"
|
| 15 |
#include "TLimitDataSource.h"
|
| 16 |
#include "TConfidenceLevel.h"
|
| 17 |
|
| 18 |
int _plotindex=0;
|
| 19 |
TH1F * MakeBackground(int evts=1000, int bins=40, double min=100, double max=500)
|
| 20 |
{
|
| 21 |
char * name = new char[100];
|
| 22 |
sprintf(name,"backgd%d",_plotindex++);
|
| 23 |
TH1F * result = new TH1F(name,";(pseudo) MET [GeV];(pseudo) events",bins, min,max);
|
| 24 |
result->Sumw2();
|
| 25 |
TF1 *fa = new TF1("fa","1.0*exp(-0.01*x)",100,500);
|
| 26 |
for (int i=0; i<evts; ++i)
|
| 27 |
result->Fill( fa->GetRandom(), 0.1 );
|
| 28 |
return result;
|
| 29 |
}
|
| 30 |
|
| 31 |
TH1F * MakeSignal(double mass, double width, int evts=100, int bins=40, double min=100, double max=500)
|
| 32 |
{
|
| 33 |
char * name = new char[100];
|
| 34 |
sprintf(name,"signal_m%f_%d",mass,_plotindex++);
|
| 35 |
TH1F * result = new TH1F(name,";(pseudo) MET [GeV];(pseudo) events",bins, min,max);
|
| 36 |
result->Sumw2();
|
| 37 |
sprintf(name,"fb_m%f",mass);
|
| 38 |
TF1 *fb = new TF1(name,"gaus(0)",100,500);
|
| 39 |
fb->SetParameter(0,1);
|
| 40 |
fb->SetParameter(1, mass );
|
| 41 |
fb->SetParameter(2, width );
|
| 42 |
for (int i=0; i<evts; ++i)
|
| 43 |
result->Fill( fb->GetRandom(), 0.1 );
|
| 44 |
return result;
|
| 45 |
|
| 46 |
}
|
| 47 |
|
| 48 |
TH1F * MakeData(TH1F*b, TH1F*s=0)
|
| 49 |
{
|
| 50 |
char * name = new char[100];
|
| 51 |
sprintf(name,"data%d",_plotindex++);
|
| 52 |
TH1F * result = new TH1F(name,";(pseudo) MET [GeV];(pseudo) events",
|
| 53 |
b->GetNbinsX(), b->GetXaxis()->GetXmin(),b->GetXaxis()->GetXmax());
|
| 54 |
TRandom3 * rand = new TRandom3(0);
|
| 55 |
for (int i=0; i<=b->GetNbinsX(); ++i){
|
| 56 |
double entry;
|
| 57 |
if (s==0)
|
| 58 |
entry = rand->Poisson( b->GetBinContent(i) );
|
| 59 |
else
|
| 60 |
entry = rand->Poisson( b->GetBinContent(i)+s->GetBinContent(i) );
|
| 61 |
result->SetBinContent(i,entry);
|
| 62 |
result->SetBinError(i,sqrt(entry));
|
| 63 |
}
|
| 64 |
return result;
|
| 65 |
}
|
| 66 |
|
| 67 |
int main()
|
| 68 |
{
|
| 69 |
gStyle->SetHistFillColor(0);
|
| 70 |
gStyle->SetPalette(1);
|
| 71 |
//gStyle->SetFillColor(0);
|
| 72 |
//gStyle->SetOptStat(kFALSE);
|
| 73 |
gStyle->SetCanvasColor(0);
|
| 74 |
gStyle->SetCanvasBorderMode(0);
|
| 75 |
gStyle->SetPadColor(0);
|
| 76 |
gStyle->SetPadBorderMode(0);
|
| 77 |
gStyle->SetFrameBorderMode(0);
|
| 78 |
|
| 79 |
gStyle->SetTitleFillColor(0);
|
| 80 |
gStyle->SetTitleBorderSize(0);
|
| 81 |
gStyle->SetTitleX(0.10);
|
| 82 |
gStyle->SetTitleY(0.98);
|
| 83 |
gStyle->SetTitleW(0.8);
|
| 84 |
gStyle->SetTitleH(0.06);
|
| 85 |
|
| 86 |
gStyle->SetErrorX(0);
|
| 87 |
gStyle->SetStatColor(0);
|
| 88 |
gStyle->SetStatBorderSize(0);
|
| 89 |
gStyle->SetStatX(0);
|
| 90 |
gStyle->SetStatY(0);
|
| 91 |
gStyle->SetStatW(0);
|
| 92 |
gStyle->SetStatH(0);
|
| 93 |
|
| 94 |
gStyle->SetTitleFont(22);
|
| 95 |
gStyle->SetLabelFont(22,"X");
|
| 96 |
gStyle->SetLabelFont(22,"Y");
|
| 97 |
gStyle->SetLabelFont(22,"Z");
|
| 98 |
gStyle->SetLabelSize(0.03,"X");
|
| 99 |
gStyle->SetLabelSize(0.03,"Y");
|
| 100 |
gStyle->SetLabelSize(0.03,"Z");
|
| 101 |
|
| 102 |
TH1F * backgd = MakeBackground(5000);
|
| 103 |
TH1F * signal = MakeSignal(300,sqrt(300), 200);
|
| 104 |
TH1F * data = MakeData( backgd );//, signal );
|
| 105 |
|
| 106 |
//Draw MET for background, signal & data:
|
| 107 |
TCanvas * c1 = new TCanvas("","",600,600);
|
| 108 |
signal->SetLineColor(2);
|
| 109 |
data->SetMarkerStyle(8);
|
| 110 |
backgd->Draw("h");
|
| 111 |
signal->Draw("h,same");
|
| 112 |
data->Draw("pe,same");
|
| 113 |
c1->SetLogy(1);
|
| 114 |
c1->SaveAs("backgd.eps");
|
| 115 |
|
| 116 |
//cout the CL's:
|
| 117 |
int fNMC = 50000;
|
| 118 |
TLimitDataSource* mydatasource = new TLimitDataSource(signal,backgd,data);
|
| 119 |
TConfidenceLevel *myconfidence = TLimit::ComputeLimit(mydatasource,fNMC);
|
| 120 |
std::cout << " CLs : " << myconfidence->CLs() << "\n"
|
| 121 |
<< " CLsb : " << myconfidence->CLsb() << "\n"
|
| 122 |
<< " CLb : " << myconfidence->CLb() << "\n"
|
| 123 |
<< "< CLs > : " << myconfidence->GetExpectedCLs_b() << "\n"
|
| 124 |
<< "< CLsb > : " << myconfidence->GetExpectedCLsb_b() << "\n"
|
| 125 |
<< "< CLb > : " << myconfidence->GetExpectedCLb_b() << "\n"
|
| 126 |
<< " -2 ln Q : " << myconfidence->GetStatistic() << "\n"
|
| 127 |
<< std::endl;
|
| 128 |
|
| 129 |
//Draw the expected and observed test statistics:
|
| 130 |
TH1F h("TConfidenceLevel_Draw","",50,0,0);
|
| 131 |
Int_t i;
|
| 132 |
double * fTSB = myconfidence->GetTestStatistic_B();
|
| 133 |
double * fTSS = myconfidence->GetTestStatistic_SB();
|
| 134 |
for (i=0; i<fNMC; i++) {
|
| 135 |
h.Fill(-2*(fTSB[i]-myconfidence->GetStot() ));
|
| 136 |
h.Fill(-2*(fTSS[i]-myconfidence->GetStot() ));
|
| 137 |
}
|
| 138 |
TH1F* b_hist = new TH1F("b_hist", ";-2 ln Q;normalized p.d.f.",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
|
| 139 |
TH1F* sb_hist = new TH1F("sb_hist",";-2 ln Q;normalized p.d.f.",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
|
| 140 |
for (i=0; i<fNMC; i++) {
|
| 141 |
b_hist->Fill(-2*(fTSB[i]-myconfidence->GetStot() ));
|
| 142 |
sb_hist->Fill(-2*(fTSS[i]-myconfidence->GetStot() ));
|
| 143 |
}
|
| 144 |
b_hist->Scale(1.0/b_hist->Integral());
|
| 145 |
sb_hist->Scale(1.0/sb_hist->Integral());
|
| 146 |
b_hist->GetYaxis()->SetTitleOffset(1.3);
|
| 147 |
b_hist->SetMinimum(0.001);
|
| 148 |
b_hist->SetMaximum(3*b_hist->GetMaximum());
|
| 149 |
b_hist->SetLineWidth(3);
|
| 150 |
sb_hist->SetLineWidth(3);
|
| 151 |
b_hist->SetLineColor(kBlue);
|
| 152 |
sb_hist->SetLineColor( 28 );
|
| 153 |
TLine * line = new TLine(myconfidence->GetStatistic(),b_hist->GetMinimum(),
|
| 154 |
myconfidence->GetStatistic(),0.6*b_hist->GetMaximum() );
|
| 155 |
line->SetLineWidth(3);
|
| 156 |
TLegend * leg = new TLegend(0.11,0.78,0.48,0.89);
|
| 157 |
leg->SetBorderSize(0);
|
| 158 |
leg->SetFillColor(0);
|
| 159 |
leg->AddEntry(line,"Observed ","l");
|
| 160 |
leg->AddEntry(b_hist,"Expected background-only ","l");
|
| 161 |
leg->AddEntry(sb_hist,"Expected signal+background ","l");
|
| 162 |
b_hist->Draw("h");
|
| 163 |
line->Draw();
|
| 164 |
leg->Draw("same");
|
| 165 |
sb_hist->Draw("h,same");
|
| 166 |
c1->SaveAs("conf.eps");
|
| 167 |
|
| 168 |
//CLs vs mass:
|
| 169 |
int nmass = 40;
|
| 170 |
TH1F * cls_obs = new TH1F("cls_obs", ";(pseudo) mass [GeV];CLs",nmass,100,500);
|
| 171 |
TH1F * cls_exp = new TH1F("cls_exp", ";(pseudo) mass [GeV];CLs",nmass,100,500);
|
| 172 |
TH1F * lnQ_obs = new TH1F("lnQ_obs", ";(pseudo) mass [GeV];-2 ln Q",nmass,100,500);
|
| 173 |
TH1F * lnQ_exp_b = new TH1F("lnQ_exp_b", ";(pseudo) mass [GeV];-2 ln Q",nmass,100,500);
|
| 174 |
TH1F * lnQ_exp_sb = new TH1F("lnQ_exp_sb", ";(pseudo) mass [GeV];-2 ln Q",nmass,100,500);
|
| 175 |
double y_cls_exp1[nmass*2];
|
| 176 |
double x_cls_exp[nmass*2];
|
| 177 |
double y_cls_exp2[nmass*2];
|
| 178 |
double y_lnQ_exp1[nmass*2];
|
| 179 |
double y_lnQ_exp2[nmass*2];
|
| 180 |
fNMC = 2000;
|
| 181 |
|
| 182 |
for (int i=0; i<=nmass; ++i){
|
| 183 |
|
| 184 |
TH1F sig( *signal );
|
| 185 |
sig.Scale( 4*exp(-i/8.) );
|
| 186 |
std::cout << i*10 << " GeV: n_signal = "<<sig.Integral()<<std::endl;
|
| 187 |
|
| 188 |
TLimitDataSource * source = new TLimitDataSource(&sig,backgd,data);
|
| 189 |
TConfidenceLevel * conf = TLimit::ComputeLimit(source,fNMC);
|
| 190 |
|
| 191 |
cls_obs->SetBinContent(i, conf->CLs() );
|
| 192 |
cls_exp->SetBinContent(i, conf->GetExpectedCLs_b() );
|
| 193 |
y_cls_exp1[i] = conf->GetExpectedCLs_b( 1);
|
| 194 |
y_cls_exp1[nmass*2-i] = conf->GetExpectedCLs_b(-1);
|
| 195 |
y_cls_exp2[i] = conf->GetExpectedCLs_b( 2);
|
| 196 |
y_cls_exp2[nmass*2-i] = conf->GetExpectedCLs_b(-2);
|
| 197 |
x_cls_exp[i] = 100+10*i;
|
| 198 |
x_cls_exp[nmass*2-i] = 100+10*i;
|
| 199 |
|
| 200 |
lnQ_obs->SetBinContent( i, conf->GetStatistic() );
|
| 201 |
lnQ_exp_b->SetBinContent( i, conf->GetExpectedStatistic_b() );
|
| 202 |
lnQ_exp_sb->SetBinContent( i, conf->GetExpectedStatistic_sb() );
|
| 203 |
y_lnQ_exp1[i] = conf->GetExpectedStatistic_b(1);
|
| 204 |
y_lnQ_exp1[nmass*2-i] = conf->GetExpectedStatistic_b(-1);
|
| 205 |
y_lnQ_exp2[i] = conf->GetExpectedStatistic_b(2);
|
| 206 |
y_lnQ_exp2[nmass*2-i] = conf->GetExpectedStatistic_b(-2);
|
| 207 |
|
| 208 |
delete conf;
|
| 209 |
delete source;
|
| 210 |
}
|
| 211 |
TGraph * cls_exp1 = new TGraph(2*nmass,x_cls_exp,y_cls_exp1);
|
| 212 |
TGraph * cls_exp2 = new TGraph(2*nmass,x_cls_exp,y_cls_exp2);
|
| 213 |
cls_exp2->SetLineColor(5);
|
| 214 |
cls_exp2->SetFillColor(5);
|
| 215 |
cls_exp1->SetLineColor(3);
|
| 216 |
cls_exp1->SetFillColor(3);
|
| 217 |
cls_obs->SetLineColor(kRed);
|
| 218 |
cls_exp->SetLineColor(kBlue);
|
| 219 |
cls_exp->SetLineStyle(3);
|
| 220 |
cls_exp->SetLineWidth(3);
|
| 221 |
cls_obs->SetLineWidth(3);
|
| 222 |
cls_obs->SetMinimum(0.00001);
|
| 223 |
TLegend * lg = new TLegend(0.5,0.28,0.89,0.4);
|
| 224 |
lg->SetBorderSize(0);
|
| 225 |
lg->SetFillColor(0);
|
| 226 |
lg->AddEntry(cls_obs,"Observed ","l");
|
| 227 |
lg->AddEntry(cls_exp,"Expected background-only ","l");
|
| 228 |
cls_obs->Draw("c");
|
| 229 |
lg->Draw("same");
|
| 230 |
cls_exp2->Draw("lf,same");
|
| 231 |
cls_exp1->Draw("lf,same");
|
| 232 |
cls_exp->Draw("c,same");
|
| 233 |
cls_obs->Draw("c,same");
|
| 234 |
line->SetLineWidth(2);
|
| 235 |
line->DrawLine(100,0.05,500,0.05);
|
| 236 |
c1->SaveAs("cls.eps");
|
| 237 |
|
| 238 |
|
| 239 |
//-2 ln Q vs mass:
|
| 240 |
TGraph * lnQ_exp1 = new TGraph(2*nmass,x_cls_exp,y_lnQ_exp1);
|
| 241 |
TGraph * lnQ_exp2 = new TGraph(2*nmass,x_cls_exp,y_lnQ_exp2);
|
| 242 |
lnQ_exp2->SetLineColor(5);
|
| 243 |
lnQ_exp2->SetFillColor(5);
|
| 244 |
lnQ_exp1->SetLineColor(3);
|
| 245 |
lnQ_exp1->SetFillColor(3);
|
| 246 |
c1->SetLogy(0);
|
| 247 |
lnQ_obs->SetLineColor(kRed);
|
| 248 |
lnQ_obs->SetLineWidth(3);
|
| 249 |
lnQ_exp_b->SetLineColor(kBlue);
|
| 250 |
lnQ_exp_b->SetLineStyle(3);
|
| 251 |
lnQ_exp_b->SetLineWidth(3);
|
| 252 |
lnQ_exp_sb->SetLineColor( 28 );
|
| 253 |
lnQ_exp_sb->SetLineStyle(3);
|
| 254 |
lnQ_exp_sb->SetLineWidth(3);
|
| 255 |
lnQ_obs->SetMinimum(-20);
|
| 256 |
line->SetLineWidth(1);
|
| 257 |
lnQ_obs->Draw("c");
|
| 258 |
lnQ_exp2->Draw("lf,same");
|
| 259 |
lnQ_exp1->Draw("lf,same");
|
| 260 |
lnQ_obs->Draw("c,same");
|
| 261 |
line->DrawLine(100,0.0,500,0.0);
|
| 262 |
lnQ_exp_b->Draw("c,same");
|
| 263 |
lnQ_exp_sb->Draw("c,same");
|
| 264 |
c1->SaveAs("lnQ.eps");
|
| 265 |
|
| 266 |
return 0;
|
| 267 |
}
|
