Logo ROOT  
Reference Guide
limit.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_math
3/// \notebooks
4/// This program demonstrates the computation of 95 % C.L. limits.
5/// It uses a set of randomly created histograms.
6///
7/// \macro_image
8/// \macro_output
9/// \macro_code
10///
11/// \author Christophe Delaere
12
13#include <iostream>
14#include "TH1.h"
15#include "THStack.h"
16#include "TCanvas.h"
17#include "TFrame.h"
18#include "TRandom2.h"
19#include "TSystem.h"
20#include "TVector.h"
21#include "TObjArray.h"
22#include "TLimit.h"
23#include "TLimitDataSource.h"
24#include "TConfidenceLevel.h"
25
26using std::cout;
27using std::endl;
28
29void limit() {
30 // Create a new canvas.
31 TCanvas *c1 = new TCanvas("c1","Dynamic Filling Example",200,10,700,500);
32 c1->SetFillColor(42);
33
34 // Create some histograms
35 TH1D* backgroundHist = new TH1D("background","The expected background",30,-4,4);
36 TH1D* signalHist = new TH1D("signal","the expected signal",30,-4,4);
37 TH1D* dataHist = new TH1D("data","some fake data points",30,-4,4);
38 backgroundHist->SetFillColor(48);
39 signalHist->SetFillColor(41);
40 dataHist->SetMarkerStyle(21);
41 dataHist->SetMarkerColor(kBlue);
42 backgroundHist->Sumw2(); // needed for stat uncertainty
43 signalHist->Sumw2(); // needed for stat uncertainty
44
45 // Fill histograms randomly
46 TRandom2 r;
47 Float_t bg,sig,dt;
48 for (Int_t i = 0; i < 25000; i++) {
49 bg = r.Gaus(0,1);
50 sig = r.Gaus(1,.2);
51 backgroundHist->Fill(bg,0.02);
52 signalHist->Fill(sig,0.001);
53 }
54 for (Int_t i = 0; i < 500; i++) {
55 dt = r.Gaus(0,1);
56 dataHist->Fill(dt);
57 }
58 THStack *hs = new THStack("hs","Signal and background compared to data...");
59 hs->Add(backgroundHist);
60 hs->Add(signalHist);
61 hs->Draw("hist");
62 dataHist->Draw("PE1,Same");
63 c1->Modified();
64 c1->Update();
65 c1->GetFrame()->SetFillColor(21);
66 c1->GetFrame()->SetBorderSize(6);
67 c1->GetFrame()->SetBorderMode(-1);
68 c1->Modified();
69 c1->Update();
71
72 // Compute the limits
73 cout << "Computing limits... " << endl;
74 TLimitDataSource* mydatasource = new TLimitDataSource(signalHist,backgroundHist,dataHist);
75 TConfidenceLevel *myconfidence = TLimit::ComputeLimit(mydatasource,50000);
76 cout << "CLs : " << myconfidence->CLs() << endl;
77 cout << "CLsb : " << myconfidence->CLsb() << endl;
78 cout << "CLb : " << myconfidence->CLb() << endl;
79 cout << "< CLs > : " << myconfidence->GetExpectedCLs_b() << endl;
80 cout << "< CLsb > : " << myconfidence->GetExpectedCLsb_b() << endl;
81 cout << "< CLb > : " << myconfidence->GetExpectedCLb_b() << endl;
82
83 // Add stat uncertainty
84 cout << endl << "Computing limits with stat systematics... " << endl;
85 TConfidenceLevel *mystatconfidence = TLimit::ComputeLimit(mydatasource,50000,true);
86 cout << "CLs : " << mystatconfidence->CLs() << endl;
87 cout << "CLsb : " << mystatconfidence->CLsb() << endl;
88 cout << "CLb : " << mystatconfidence->CLb() << endl;
89 cout << "< CLs > : " << mystatconfidence->GetExpectedCLs_b() << endl;
90 cout << "< CLsb > : " << mystatconfidence->GetExpectedCLsb_b() << endl;
91 cout << "< CLb > : " << mystatconfidence->GetExpectedCLb_b() << endl;
92
93 // Add some systematics
94 cout << endl << "Computing limits with systematics... " << endl;
95 TVectorD errorb(2);
96 TVectorD errors(2);
97 TObjArray* names = new TObjArray();
98 TObjString name1("bg uncertainty");
99 TObjString name2("sig uncertainty");
100 names->AddLast(&name1);
101 names->AddLast(&name2);
102 errorb[0]=0.05; // error source 1: 5%
103 errorb[1]=0; // error source 2: 0%
104 errors[0]=0; // error source 1: 0%
105 errors[1]=0.01; // error source 2: 1%
106 TLimitDataSource* mynewdatasource = new TLimitDataSource();
107 mynewdatasource->AddChannel(signalHist,backgroundHist,dataHist,&errors,&errorb,names);
108 TConfidenceLevel *mynewconfidence = TLimit::ComputeLimit(mynewdatasource,50000,true);
109 cout << "CLs : " << mynewconfidence->CLs() << endl;
110 cout << "CLsb : " << mynewconfidence->CLsb() << endl;
111 cout << "CLb : " << mynewconfidence->CLb() << endl;
112 cout << "< CLs > : " << mynewconfidence->GetExpectedCLs_b() << endl;
113 cout << "< CLsb > : " << mynewconfidence->GetExpectedCLsb_b() << endl;
114 cout << "< CLb > : " << mynewconfidence->GetExpectedCLb_b() << endl;
115
116 // show canonical -2lnQ plots in a new canvas
117 // - The histogram of -2lnQ for background hypothesis (full)
118 // - The histogram of -2lnQ for signal and background hypothesis (dashed)
119 TCanvas *c2 = new TCanvas("c2");
120 myconfidence->Draw();
121
122 // clean up (except histograms and canvas)
123 delete myconfidence;
124 delete mydatasource;
125 delete mystatconfidence;
126 delete mynewconfidence;
127 delete mynewdatasource;
128}
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