From e4db8ed8aa2440f81d20b051d1c1713f798e9be4 Mon Sep 17 00:00:00 2001
From: Rene Brun <Rene.Brun@cern.ch>
Date: Fri, 21 Mar 2003 14:53:50 +0000
Subject: [PATCH] Fix from Christophe Delaere. All TH1F and TArrayF are
replaced by TH1D and TArrayD respectively. This solves some numerical
precision problems observed with the float implementation
git-svn-id: http://root.cern.ch/svn/root/trunk@6337 27541ba8-7e3a-0410-8455-c3a389f83636
---
hist/inc/TLimit.h | 8 ++---
hist/inc/TLimitDataSource.h | 14 ++++-----
hist/src/TLimit.cxx | 58 +++++++++++++++++------------------
hist/src/TLimitDataSource.cxx | 28 ++++++++---------
tutorials/limit.C | 10 +++---
5 files changed, 59 insertions(+), 59 deletions(-)
diff --git a/hist/inc/TLimit.h b/hist/inc/TLimit.h
index 2925b548f5d..f7d7ac4c50e 100644
--- a/hist/inc/TLimit.h
+++ b/hist/inc/TLimit.h
@@ -1,4 +1,4 @@
-// @(#)root/hist:$Name: $:$Id: TLimit.h,v 1.34 2002/08/16 21:16:00 brun Exp $
+// @(#)root/hist:$Name: $:$Id: TLimit.h,v 1.1 2002/09/06 19:57:59 brun Exp $
// Author: Christophe.Delaere@cern.ch 21/08/2002
#ifndef ROOT_TLimit
@@ -14,7 +14,7 @@
class TConfidenceLevel;
class TRandom;
class TLimitDataSource;
-class TArrayF;
+class TArrayD;
class TOrdCollection;
//____________________________________________________________________
@@ -41,9 +41,9 @@ class TLimit {
static TLimitDataSource *Fluctuate(TLimitDataSource * input, bool init,TRandom *);
inline static Double_t LogLikelihood(Double_t s, Double_t b, Double_t d) { return d * TMath::Log(1 + (s / b)); }
private:
- static TArrayF *fgTable; // a log table... just to speed up calculation
+ static TArrayD *fgTable; // a log table... just to speed up calculation
static TOrdCollection *fgSystNames; // Collection of systematics names
- ClassDef(TLimit, 1) // Class to compute 95% CL limits
+ ClassDef(TLimit, 2) // Class to compute 95% CL limits
};
#endif
diff --git a/hist/inc/TLimitDataSource.h b/hist/inc/TLimitDataSource.h
index dab5f06a1f5..5bb9cc4586b 100644
--- a/hist/inc/TLimitDataSource.h
+++ b/hist/inc/TLimitDataSource.h
@@ -1,4 +1,4 @@
-// @(#)root/hist:$Name: $:$Id: TLimitDataSource.h,v 1.34 2002/08/16 21:16:00 brun Exp $
+// @(#)root/hist:$Name: $:$Id: TLimitDataSource.h,v 1.1 2002/09/06 19:57:59 brun Exp $
// Author: Christophe.Delaere@cern.ch 21/08/2002
#ifndef ROOT_TLimitDataSource
@@ -8,7 +8,7 @@
#include "TObjArray.h"
#endif
-class TH1F;
+class TH1D;
//_______________________________________________________________________
//
@@ -25,9 +25,9 @@ class TLimitDataSource {
public:
TLimitDataSource();
virtual ~TLimitDataSource() {}
- TLimitDataSource(TH1F* s,TH1F* b,TH1F* d);
- virtual void AddChannel(TH1F*,TH1F*,TH1F*);
- virtual void AddChannel(TH1F*,TH1F*,TH1F*,TH1F*, TH1F*, TObjArray*);
+ TLimitDataSource(TH1D* s,TH1D* b,TH1D* d);
+ virtual void AddChannel(TH1D*,TH1D*,TH1D*);
+ virtual void AddChannel(TH1D*,TH1D*,TH1D*,TH1D*, TH1D*, TObjArray*);
inline virtual TObjArray* GetSignal() { return &fSignal;}
inline virtual TObjArray* GetBackground() { return &fBackground;}
inline virtual TObjArray* GetCandidates() { return &fCandidates;}
@@ -44,10 +44,10 @@ private:
TObjArray fErrorOnBackground; //packed error sources for background
TObjArray fIds; //packed IDs for the different error sources
// some dummy objects that the class will use and delete
- TObjArray fDummyTH1F; //array of dummy object (used for bookeeping)
+ TObjArray fDummyTH1D; //array of dummy object (used for bookeeping)
TObjArray fDummyIds; //array of dummy object (used for bookeeping)
- ClassDef(TLimitDataSource, 1 ) // input for TLimit routines
+ ClassDef(TLimitDataSource, 2 ) // input for TLimit routines
};
#endif
diff --git a/hist/src/TLimit.cxx b/hist/src/TLimit.cxx
index 1f0289d19b2..c89ede1c35d 100644
--- a/hist/src/TLimit.cxx
+++ b/hist/src/TLimit.cxx
@@ -1,4 +1,4 @@
-// @(#)root/hist:$Name: $:$Id: TLimit.cxx,v 1.3 2002/09/06 21:41:23 brun Exp $
+// @(#)root/hist:$Name: $:$Id: TLimit.cxx,v 1.4 2002/09/10 14:59:15 rdm Exp $
// Author: Christophe.Delaere@cern.ch 21/08/2002
///////////////////////////////////////////////////////////////////////////
@@ -16,7 +16,7 @@
*************************************************************************/
#include "TLimit.h"
-#include "TArrayF.h"
+#include "TArrayD.h"
#include "TOrdCollection.h"
#include "TConfidenceLevel.h"
#include "TLimitDataSource.h"
@@ -31,7 +31,7 @@
ClassImp(TLimit)
-TArrayF *TLimit::fgTable = new TArrayF(0);
+TArrayD *TLimit::fgTable = new TArrayD(0);
TOrdCollection *TLimit::fgSystNames = new TOrdCollection();
TConfidenceLevel *TLimit::ComputeLimit(TLimitDataSource * data,
@@ -78,9 +78,9 @@ TConfidenceLevel *TLimit::ComputeLimit(TLimitDataSource * data,
<BLOCKQUOTE><PRE>
TFile* infile=new TFile("plotfile.root","READ");
infile->cd();
- TH1F* sh=(TH1F*)infile->Get("signal");
- TH1F* bh=(TH1F*)infile->Get("background");
- TH1F* dh=(TH1F*)infile->Get("data");
+ TH1D* sh=(TH1D*)infile->Get("signal");
+ TH1D* bh=(TH1D*)infile->Get("background");
+ TH1D* dh=(TH1D*)infile->Get("data");
TLimitDataSource* mydatasource = new TLimitDataSource(sh,bh,dh);
TConfidenceLevel *myconfidence = TLimit::ComputeLimit(mydatasource,50000);
cout << " CLs : " << myconfidence->CLs() << endl;
@@ -112,12 +112,12 @@ TConfidenceLevel *TLimit::ComputeLimit(TLimitDataSource * data,
Int_t ncand = 0;
Int_t i;
for (i = 0; i <= data->GetSignal()->GetLast(); i++) {
- nbins += ((TH1F *) (data->GetSignal()->At(i)))->GetNbinsX();
- maxbins = ((TH1F *) (data->GetSignal()->At(i)))->GetNbinsX() > maxbins ?
- ((TH1F *) (data->GetSignal()->At(i)))->GetNbinsX() + 1 : maxbins;
- nsig += ((TH1F *) (data->GetSignal()->At(i)))->Integral();
- nbg += ((TH1F *) (data->GetBackground()->At(i)))->Integral();
- ncand += (Int_t) ((TH1F *) (data->GetCandidates()->At(i)))->Integral();
+ nbins += ((TH1D *) (data->GetSignal()->At(i)))->GetNbinsX();
+ maxbins = ((TH1D *) (data->GetSignal()->At(i)))->GetNbinsX() > maxbins ?
+ ((TH1D *) (data->GetSignal()->At(i)))->GetNbinsX() + 1 : maxbins;
+ nsig += ((TH1D *) (data->GetSignal()->At(i)))->Integral();
+ nbg += ((TH1D *) (data->GetBackground()->At(i)))->Integral();
+ ncand += (Int_t) ((TH1D *) (data->GetCandidates()->At(i)))->Integral();
}
result->SetBtot(nbg);
result->SetStot(nsig);
@@ -126,11 +126,11 @@ TConfidenceLevel *TLimit::ComputeLimit(TLimitDataSource * data,
fgTable->Set(maxbins * (data->GetSignal()->GetLast() + 1));
for (Int_t channel = 0; channel <= data->GetSignal()->GetLast(); channel++)
for (Int_t bin = 0;
- bin <= ((TH1F *) (data->GetSignal()->At(channel)))->GetNbinsX();
+ bin <= ((TH1D *) (data->GetSignal()->At(channel)))->GetNbinsX();
bin++) {
- Double_t s = (Double_t) ((TH1F *) (data->GetSignal()->At(channel)))->GetBinContent(bin);
- Double_t b = (Double_t) ((TH1F *) (data->GetBackground()->At(channel)))->GetBinContent(bin);
- Double_t d = (Double_t) ((TH1F *) (data->GetCandidates()->At(channel)))->GetBinContent(bin);
+ Double_t s = (Double_t) ((TH1D *) (data->GetSignal()->At(channel)))->GetBinContent(bin);
+ Double_t b = (Double_t) ((TH1D *) (data->GetBackground()->At(channel)))->GetBinContent(bin);
+ Double_t d = (Double_t) ((TH1D *) (data->GetCandidates()->At(channel)))->GetBinContent(bin);
// Compute the value of the "-2lnQ" for the actual data
if ((b == 0) && (s > 0)) {
cout << "WARNING: Ignoring bin " << bin << " of channel "
@@ -168,22 +168,22 @@ TConfidenceLevel *TLimit::ComputeLimit(TLimitDataSource * data,
for (Int_t channel = 0;
channel <= fluctuated->GetSignal()->GetLast(); channel++) {
for (Int_t bin = 0;
- bin <=((TH1F *) (fluctuated->GetSignal()->At(channel)))->GetNbinsX();
+ bin <=((TH1D *) (fluctuated->GetSignal()->At(channel)))->GetNbinsX();
bin++) {
- if ((Double_t) ((TH1F *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin) != 0) {
+ if ((Double_t) ((TH1D *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin) != 0) {
// s+b hypothesis
- Double_t rate = (Double_t) ((TH1F *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin) +
- (Double_t) ((TH1F *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
+ Double_t rate = (Double_t) ((TH1D *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin) +
+ (Double_t) ((TH1D *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
Double_t rand = myrandom->Poisson(rate);
tss[i] += rand * fgTable->At((channel * maxbins) + bin);
- Double_t s = (Double_t) ((TH1F *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin);
- Double_t b = (Double_t) ((TH1F *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
+ Double_t s = (Double_t) ((TH1D *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin);
+ Double_t b = (Double_t) ((TH1D *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
if ((s > 0) && (b > 0))
lrs[i] += statistic(s, b, rand) - s;
else if ((s > 0) && (b == 0))
lrs[i] += 20 * rand - s;
// b hypothesis
- rate = (Double_t) ((TH1F *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
+ rate = (Double_t) ((TH1D *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
rand = myrandom->Poisson(rate);
tsb[i] += rand * fgTable->At((channel * maxbins) + bin);
if ((s > 0) && (b > 0))
@@ -254,11 +254,11 @@ TLimitDataSource *TLimit::Fluctuate(TLimitDataSource * input, bool init,
serrf[channel] = 0;
berrf[channel] = 0;
for (Int_t bin = 0;
- bin <=((TH1F *) (input->GetErrorOnSignal()->At(channel)))->GetNbinsX();
+ bin <=((TH1D *) (input->GetErrorOnSignal()->At(channel)))->GetNbinsX();
bin++) {
- serrf[channel] += ((TH1F *) (input->GetErrorOnSignal()->At(channel)))->GetBinContent(bin) *
+ serrf[channel] += ((TH1D *) (input->GetErrorOnSignal()->At(channel)))->GetBinContent(bin) *
toss[fgSystNames->BinarySearch((TObjString*) (((TObjArray *) (input->GetErrorNames()->At(channel)))->At(bin)))];
- berrf[channel] += ((TH1F *) (input->GetErrorOnBackground()->At(channel)))->GetBinContent(bin) *
+ berrf[channel] += ((TH1D *) (input->GetErrorOnBackground()->At(channel)))->GetBinContent(bin) *
toss[fgSystNames->BinarySearch((TObjString*) (((TObjArray *) (input->GetErrorNames()->At(channel)))->At(bin)))];
}
if ((serrf[channel] < -1.0) || (berrf[channel] < -0.9)) {
@@ -274,13 +274,13 @@ TLimitDataSource *TLimit::Fluctuate(TLimitDataSource * input, bool init,
result->SetOwner();
for (Int_t channel = 0; channel <= input->GetSignal()->GetLast();
channel++) {
- TH1F *newsignal = new TH1F(*(TH1F *) (input->GetSignal()->At(channel)));
+ TH1D *newsignal = new TH1D(*(TH1D *) (input->GetSignal()->At(channel)));
newsignal->Scale(1 + serrf[channel]);
newsignal->SetDirectory(0);
- TH1F *newbackground = new TH1F(*(TH1F *) (input->GetBackground()->At(channel)));
+ TH1D *newbackground = new TH1D(*(TH1D *) (input->GetBackground()->At(channel)));
newbackground->Scale(1 + berrf[channel]);
newbackground->SetDirectory(0);
- TH1F *newcandidates = new TH1F(*(TH1F *) (input->GetCandidates()));
+ TH1D *newcandidates = new TH1D(*(TH1D *) (input->GetCandidates()));
newcandidates->SetDirectory(0);
result->AddChannel(newsignal, newbackground, newcandidates);
}
diff --git a/hist/src/TLimitDataSource.cxx b/hist/src/TLimitDataSource.cxx
index 0455cb2e60f..fdf6126a442 100644
--- a/hist/src/TLimitDataSource.cxx
+++ b/hist/src/TLimitDataSource.cxx
@@ -1,4 +1,4 @@
-// @(#)root/hist:$Name: $:$Id: TLimitDataSource.cxx,v 1.1 2002/09/06 19:58:00 brun Exp $
+// @(#)root/hist:$Name: $:$Id: TLimitDataSource.cxx,v 1.2 2002/09/13 15:23:56 brun Exp $
// Author: Christophe.Delaere@cern.ch 21/08/2002
///////////////////////////////////////////////////////////////////////////
@@ -19,50 +19,50 @@ ClassImp(TLimitDataSource)
TLimitDataSource::TLimitDataSource()
{
// Default constructor
- fDummyTH1F.SetOwner();
+ fDummyTH1D.SetOwner();
fDummyIds.SetOwner();
}
-TLimitDataSource::TLimitDataSource(TH1F * s, TH1F * b, TH1F * d)
+TLimitDataSource::TLimitDataSource(TH1D * s, TH1D * b, TH1D * d)
{
// Another constructor, directly adds one channel
// with signal, background and data given as input.
- fDummyTH1F.SetOwner();
+ fDummyTH1D.SetOwner();
fDummyIds.SetOwner();
AddChannel(s, b, d);
}
-void TLimitDataSource::AddChannel(TH1F * s, TH1F * b, TH1F * d)
+void TLimitDataSource::AddChannel(TH1D * s, TH1D * b, TH1D * d)
{
// Adds a channel with signal, background and data given as input.
- TH1F *empty;
+ TH1D *empty;
TRandom3 generator;
fSignal.AddLast(s);
fBackground.AddLast(b);
fCandidates.AddLast(d);
char rndname[20];
sprintf(rndname, "rndname%f", generator.Rndm());
- empty = new TH1F(rndname, "", s->GetSize(), 0, 1);
+ empty = new TH1D(rndname, "", s->GetSize(), 0, 1);
empty->SetDirectory(0);
fErrorOnSignal.AddLast(empty);
- fDummyTH1F.AddLast(empty);
+ fDummyTH1D.AddLast(empty);
sprintf(rndname, "rndname%f", generator.Rndm());
- empty = new TH1F(rndname, "", s->GetSize(), 0, 1);
+ empty = new TH1D(rndname, "", s->GetSize(), 0, 1);
empty->SetDirectory(0);
fErrorOnBackground.AddLast(empty);
- fDummyTH1F.AddLast(empty);
+ fDummyTH1D.AddLast(empty);
TObjArray *dummy = new TObjArray(0);
fIds.AddLast(dummy);
fDummyIds.AddLast(dummy);
}
-void TLimitDataSource::AddChannel(TH1F * s, TH1F * b, TH1F * d, TH1F * es,
- TH1F * eb, TObjArray * names)
+void TLimitDataSource::AddChannel(TH1D * s, TH1D * b, TH1D * d, TH1D * es,
+ TH1D * eb, TObjArray * names)
{
// Adds a channel with signal, background and data given as input.
// In addition, error sources are defined.
- // TH1F are here used for convenience: each bin has to be seen as
+ // TH1D are here used for convenience: each bin has to be seen as
// an error source (relative).
// names is an array of strings containing the names of the sources.
// Sources with the same name are correlated.
@@ -87,7 +87,7 @@ void TLimitDataSource::SetOwner(bool swtch)
fErrorOnSignal.SetOwner(swtch);
fErrorOnBackground.SetOwner(swtch);
fIds.SetOwner(swtch);
- fDummyTH1F.SetOwner(!swtch);
+ fDummyTH1D.SetOwner(!swtch);
fDummyIds.SetOwner(!swtch);
}
diff --git a/tutorials/limit.C b/tutorials/limit.C
index 50cd833a987..ae9e40c462d 100644
--- a/tutorials/limit.C
+++ b/tutorials/limit.C
@@ -13,9 +13,9 @@ void limit() {
c1->SetFillColor(42);
// Create some histograms
- TH1F* background = new TH1F("background","The expected background",30,-4,4);
- TH1F* signal = new TH1F("signal","the expected signal",30,-4,4);
- TH1F* data = new TH1F("data","some fake data points",30,-4,4);
+ TH1D* background = new TH1D("background","The expected background",30,-4,4);
+ TH1D* signal = new TH1D("signal","the expected signal",30,-4,4);
+ TH1D* data = new TH1D("data","some fake data points",30,-4,4);
background->SetFillColor(48);
signal->SetFillColor(41);
data->SetMarkerStyle(21);
@@ -61,8 +61,8 @@ void limit() {
// Add some systematics
cout << endl << "Computing limits with systematics... " << endl;
- TH1F* errorb = new TH1F("errorb","errors on background",1,0,1);
- TH1F* errors = new TH1F("errors","errors on signal",1,0,1);
+ TH1D* errorb = new TH1D("errorb","errors on background",1,0,1);
+ TH1D* errors = new TH1D("errors","errors on signal",1,0,1);
TObjArray* names = new TObjArray();
TObjString name1("bg uncertainty");
TObjString name2("sig uncertainty");
--
GitLab