Logo ROOT  
Reference Guide
limit.C File Reference

Detailed Description

View in nbviewer Open in SWAN This program demonstrates the computation of 95 % C.L.

limits. It uses a set of randomly created histograms.

Computing limits...
CLs : 0.0179535
CLsb : 0.00930313
CLb : 0.51818
< CLs > : 0.0165344
< CLsb > : 0.00826754
< CLb > : 0.50002
Computing limits with stat systematics...
CLs : 0.0189723
CLsb : 0.00989482
CLb : 0.52154
< CLs > : 0.0172591
< CLsb > : 0.00862989
< CLb > : 0.50002
Computing limits with systematics...
CLs : 0.0352037
CLsb : 0.0184327
CLb : 0.5236
< CLs > : 0.0328702
< CLsb > : 0.0164358
< CLb > : 0.50002
#include <iostream>
#include "TH1.h"
#include "THStack.h"
#include "TCanvas.h"
#include "TFrame.h"
#include "TRandom2.h"
#include "TSystem.h"
#include "TVector.h"
#include "TObjArray.h"
#include "TLimit.h"
using std::cout;
using std::endl;
void limit() {
// Create a new canvas.
TCanvas *c1 = new TCanvas("c1","Dynamic Filling Example",200,10,700,500);
c1->SetFillColor(42);
// Create some histograms
TH1D* backgroundHist = new TH1D("background","The expected background",30,-4,4);
TH1D* signalHist = new TH1D("signal","the expected signal",30,-4,4);
TH1D* dataHist = new TH1D("data","some fake data points",30,-4,4);
backgroundHist->SetFillColor(48);
signalHist->SetFillColor(41);
dataHist->SetMarkerStyle(21);
dataHist->SetMarkerColor(kBlue);
backgroundHist->Sumw2(); // needed for stat uncertainty
signalHist->Sumw2(); // needed for stat uncertainty
// Fill histograms randomly
Float_t bg,sig,dt;
for (Int_t i = 0; i < 25000; i++) {
bg = r.Gaus(0,1);
sig = r.Gaus(1,.2);
backgroundHist->Fill(bg,0.02);
signalHist->Fill(sig,0.001);
}
for (Int_t i = 0; i < 500; i++) {
dt = r.Gaus(0,1);
dataHist->Fill(dt);
}
THStack *hs = new THStack("hs","Signal and background compared to data...");
hs->Add(backgroundHist);
hs->Add(signalHist);
hs->Draw("hist");
dataHist->Draw("PE1,Same");
c1->Modified();
c1->Update();
c1->GetFrame()->SetFillColor(21);
c1->GetFrame()->SetBorderSize(6);
c1->GetFrame()->SetBorderMode(-1);
c1->Modified();
c1->Update();
// Compute the limits
cout << "Computing limits... " << endl;
TLimitDataSource* mydatasource = new TLimitDataSource(signalHist,backgroundHist,dataHist);
TConfidenceLevel *myconfidence = TLimit::ComputeLimit(mydatasource,50000);
cout << "CLs : " << myconfidence->CLs() << endl;
cout << "CLsb : " << myconfidence->CLsb() << endl;
cout << "CLb : " << myconfidence->CLb() << endl;
cout << "< CLs > : " << myconfidence->GetExpectedCLs_b() << endl;
cout << "< CLsb > : " << myconfidence->GetExpectedCLsb_b() << endl;
cout << "< CLb > : " << myconfidence->GetExpectedCLb_b() << endl;
// Add stat uncertainty
cout << endl << "Computing limits with stat systematics... " << endl;
TConfidenceLevel *mystatconfidence = TLimit::ComputeLimit(mydatasource,50000,true);
cout << "CLs : " << mystatconfidence->CLs() << endl;
cout << "CLsb : " << mystatconfidence->CLsb() << endl;
cout << "CLb : " << mystatconfidence->CLb() << endl;
cout << "< CLs > : " << mystatconfidence->GetExpectedCLs_b() << endl;
cout << "< CLsb > : " << mystatconfidence->GetExpectedCLsb_b() << endl;
cout << "< CLb > : " << mystatconfidence->GetExpectedCLb_b() << endl;
// Add some systematics
cout << endl << "Computing limits with systematics... " << endl;
TVectorD errorb(2);
TVectorD errors(2);
TObjArray* names = new TObjArray();
TObjString name1("bg uncertainty");
TObjString name2("sig uncertainty");
names->AddLast(&name1);
names->AddLast(&name2);
errorb[0]=0.05; // error source 1: 5%
errorb[1]=0; // error source 2: 0%
errors[0]=0; // error source 1: 0%
errors[1]=0.01; // error source 2: 1%
TLimitDataSource* mynewdatasource = new TLimitDataSource();
mynewdatasource->AddChannel(signalHist,backgroundHist,dataHist,&errors,&errorb,names);
TConfidenceLevel *mynewconfidence = TLimit::ComputeLimit(mynewdatasource,50000,true);
cout << "CLs : " << mynewconfidence->CLs() << endl;
cout << "CLsb : " << mynewconfidence->CLsb() << endl;
cout << "CLb : " << mynewconfidence->CLb() << endl;
cout << "< CLs > : " << mynewconfidence->GetExpectedCLs_b() << endl;
cout << "< CLsb > : " << mynewconfidence->GetExpectedCLsb_b() << endl;
cout << "< CLb > : " << mynewconfidence->GetExpectedCLb_b() << endl;
// show canonical -2lnQ plots in a new canvas
// - The histogram of -2lnQ for background hypothesis (full)
// - The histogram of -2lnQ for signal and background hypothesis (dashed)
TCanvas *c2 = new TCanvas("c2");
myconfidence->Draw();
// clean up (except histograms and canvas)
delete myconfidence;
delete mydatasource;
delete mystatconfidence;
delete mynewconfidence;
delete mynewdatasource;
}
ROOT::R::TRInterface & r
Definition: Object.C:4
int Int_t
Definition: RtypesCore.h:43
float Float_t
Definition: RtypesCore.h:55
@ kBlue
Definition: Rtypes.h:64
R__EXTERN TSystem * gSystem
Definition: TSystem.h:556
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
The Canvas class.
Definition: TCanvas.h:27
Class to compute 95% CL limits.
Double_t GetExpectedCLs_b(Int_t sigma=0) const
Double_t GetExpectedCLb_b(Int_t sigma=0) const
Get the expected Confidence Level for the background only if there is only background.
Double_t CLsb(bool use_sMC=kFALSE) const
Get the Confidence Level for the signal plus background hypothesis.
Double_t CLb(bool use_sMC=kFALSE) const
Get the Confidence Level for the background only.
Double_t GetExpectedCLsb_b(Int_t sigma=0) const
Get the expected Confidence Level for the signal plus background hypothesis if there is only backgrou...
Double_t CLs(bool use_sMC=kFALSE) const
Get the Confidence Level defined by CLs = CLsb/CLb.
void Draw(const Option_t *option="")
Display sort of a "canonical" -2lnQ plot.
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8476
The Histogram stack class.
Definition: THStack.h:31
virtual void Draw(Option_t *chopt="")
Draw this multihist with its current attributes.
Definition: THStack.cxx:445
virtual void Add(TH1 *h, Option_t *option="")
add a new histogram to the list Only 1-d and 2-d histograms currently supported.
Definition: THStack.cxx:359
This class serves as interface to feed data into the TLimit routines.
virtual void AddChannel(TH1 *, TH1 *, TH1 *)
static TConfidenceLevel * ComputeLimit(TLimitDataSource *data, Int_t nmc=50000, bool stat=false, TRandom *generator=0)
Definition: TLimit.cxx:103
An array of TObjects.
Definition: TObjArray.h:37
virtual void AddLast(TObject *obj)
Add object in the next empty slot in the array.
Definition: TObjArray.cxx:178
Collectable string class.
Definition: TObjString.h:28
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:27
virtual Bool_t ProcessEvents()
Process pending events (GUI, timers, sockets).
Definition: TSystem.cxx:414
return c1
Definition: legend1.C:41
return c2
Definition: legend2.C:14
Author
Christophe Delaere

Definition in file limit.C.