#include "TLimit.h"
#include "TArrayD.h"
#include "TVectorD.h"
#include "TOrdCollection.h"
#include "TConfidenceLevel.h"
#include "TLimitDataSource.h"
#include "TRandom3.h"
#include "TH1.h"
#include "TObjArray.h"
#include "TMath.h"
#include "TIterator.h"
#include "TObjString.h"
#include "TClassTable.h"
#include "Riostream.h"
ClassImp(TLimit)
TArrayD *TLimit::fgTable = new TArrayD(0);
TOrdCollection *TLimit::fgSystNames = new TOrdCollection();
TConfidenceLevel *TLimit::ComputeLimit(TLimitDataSource * data,
Int_t nmc, bool stat,
TRandom * generator)
{
/*
<FONT SIZE=+0>
<p>Supposing that there is a plotfile.root file containing 3 histograms
(signal, background and data), you can imagine doing things like:</p>
<p>
<BLOCKQUOTE><PRE>
TFile* infile=new TFile("plotfile.root","READ");
infile->cd();
TH1* sh=(TH1*)infile->Get("signal");
TH1* bh=(TH1*)infile->Get("background");
TH1* dh=(TH1*)infile->Get("data");
TLimitDataSource* mydatasource = new TLimitDataSource(sh,bh,dh);
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;
delete myconfidence;
delete mydatasource;
infile->Close();
</PRE></BLOCKQUOTE></p>
<p></p>
<p>More information can still be found on
<a HREF="http://cern.ch/aleph-proj-alphapp/doc/tlimit.html">this</a> page.</p>
</FONT>
*/
//End_Html
TConfidenceLevel *result = new TConfidenceLevel(nmc);
TRandom *myrandom = generator ? generator : new TRandom3;
Int_t maxbins = 0;
Double_t nsig = 0;
Double_t nbg = 0;
Int_t ncand = 0;
Int_t i;
for (i = 0; i <= data->GetSignal()->GetLast(); i++) {
maxbins = (((TH1 *) (data->GetSignal()->At(i)))->GetNbinsX() + 2) > maxbins ?
(((TH1 *) (data->GetSignal()->At(i)))->GetNbinsX() + 2) : maxbins;
nsig += ((TH1 *) (data->GetSignal()->At(i)))->Integral();
nbg += ((TH1 *) (data->GetBackground()->At(i)))->Integral();
ncand += (Int_t) ((TH1 *) (data->GetCandidates()->At(i)))->Integral();
}
result->SetBtot(nbg);
result->SetStot(nsig);
result->SetDtot(ncand);
Double_t buffer = 0;
fgTable->Set(maxbins * (data->GetSignal()->GetLast() + 1));
for (Int_t channel = 0; channel <= data->GetSignal()->GetLast(); channel++)
for (Int_t bin = 0;
bin <= ((TH1 *) (data->GetSignal()->At(channel)))->GetNbinsX()+1;
bin++) {
Double_t s = (Double_t) ((TH1 *) (data->GetSignal()->At(channel)))->GetBinContent(bin);
Double_t b = (Double_t) ((TH1 *) (data->GetBackground()->At(channel)))->GetBinContent(bin);
Double_t d = (Double_t) ((TH1 *) (data->GetCandidates()->At(channel)))->GetBinContent(bin);
if ((b == 0) && (s > 0)) {
cout << "WARNING: Ignoring bin " << bin << " of channel "
<< channel << " which has s=" << s << " but b=" << b << endl;
cout << " Maybe the MC statistic has to be improved..." << endl;
}
if ((s > 0) && (b > 0))
buffer += LogLikelihood(s, b, b, d);
if ((s > 0) && (b > 0))
fgTable->AddAt(LogLikelihood(s, b, b, 1), (channel * maxbins) + bin);
else if ((s > 0) && (b == 0))
fgTable->AddAt(20, (channel * maxbins) + bin);
}
result->SetTSD(buffer);
Double_t *tss = new Double_t[nmc];
Double_t *tsb = new Double_t[nmc];
Double_t *lrs = new Double_t[nmc];
Double_t *lrb = new Double_t[nmc];
TLimitDataSource* tmp_src = (TLimitDataSource*)(data->Clone());
TLimitDataSource* tmp_src2 = (TLimitDataSource*)(data->Clone());
for (i = 0; i < nmc; i++) {
tss[i] = 0;
tsb[i] = 0;
lrs[i] = 0;
lrb[i] = 0;
TLimitDataSource* fluctuated = Fluctuate(data, tmp_src, !i, myrandom, stat) ? tmp_src : data;
TLimitDataSource* fluctuated2 = Fluctuate(data, tmp_src2, false, myrandom, stat) ? tmp_src2 : data;
for (Int_t channel = 0;
channel <= fluctuated->GetSignal()->GetLast(); channel++) {
for (Int_t bin = 0;
bin <=((TH1 *) (fluctuated->GetSignal()->At(channel)))->GetNbinsX()+1;
bin++) {
if ((Double_t) ((TH1 *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin) != 0) {
Double_t rate = (Double_t) ((TH1 *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin) +
(Double_t) ((TH1 *) (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) ((TH1 *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin);
Double_t s2= (Double_t) ((TH1 *) (fluctuated2->GetSignal()->At(channel)))->GetBinContent(bin);
Double_t b = (Double_t) ((TH1 *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
Double_t b2= (Double_t) ((TH1 *) (fluctuated2->GetBackground()->At(channel)))->GetBinContent(bin);
if ((s > 0) && (b2 > 0))
lrs[i] += LogLikelihood(s, b, b2, rand) - s - b + b2;
else if ((s > 0) && (b2 == 0))
lrs[i] += 20 * rand - s;
rate = (Double_t) ((TH1 *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
rand = myrandom->Poisson(rate);
tsb[i] += rand * fgTable->At((channel * maxbins) + bin);
if ((s2 > 0) && (b > 0))
lrb[i] += LogLikelihood(s2, b2, b, rand) - s2 - b2 + b;
else if ((s > 0) && (b == 0))
lrb[i] += 20 * rand - s;
}
}
}
lrs[i] = TMath::Exp(lrs[i] < 710 ? lrs[i] : 710);
lrb[i] = TMath::Exp(lrb[i] < 710 ? lrb[i] : 710);
}
delete tmp_src;
delete tmp_src2;
result->SetTSS(tss);
result->SetTSB(tsb);
result->SetLRS(lrs);
result->SetLRB(lrb);
if (!generator)
delete myrandom;
return result;
}
bool TLimit::Fluctuate(TLimitDataSource * input, TLimitDataSource * output,
bool init, TRandom * generator, bool stat)
{
if (init) {
TIterator *errornames = input->GetErrorNames()->MakeIterator();
TObjArray *listofnames = 0;
delete fgSystNames;
fgSystNames = new TOrdCollection();
while ((listofnames = ((TObjArray *) errornames->Next()))) {
TObjString *name = 0;
TIterator *loniter = listofnames->MakeIterator();
while ((name = (TObjString *) (loniter->Next())))
if ((fgSystNames->IndexOf(name)) < 0)
fgSystNames->AddLast(name);
}
fgSystNames->Sort();
}
if (!output)
output = (TLimitDataSource*)(input->Clone());
if ((fgSystNames->GetSize() <= 0)&&(!stat)) {
return 0;
}
if (fgSystNames->GetSize() <= 0) {
output->SetOwner();
for (Int_t channel = 0; channel <= input->GetSignal()->GetLast(); channel++) {
TH1 *newsignal = (TH1*)(output->GetSignal()->At(channel));
TH1 *oldsignal = (TH1*)(input->GetSignal()->At(channel));
if(stat)
for(int i=1; i<=newsignal->GetNbinsX(); i++) {
newsignal->SetBinContent(i,oldsignal->GetBinContent(i)+generator->Gaus(0,oldsignal->GetBinError(i)));
}
newsignal->SetDirectory(0);
TH1 *newbackground = (TH1*)(output->GetBackground()->At(channel));
TH1 *oldbackground = (TH1*)(input->GetBackground()->At(channel));
if(stat)
for(int i=1; i<=newbackground->GetNbinsX(); i++)
newbackground->SetBinContent(i,oldbackground->GetBinContent(i)+generator->Gaus(0,oldbackground->GetBinError(i)));
newbackground->SetDirectory(0);
}
return 1;
}
Bool_t retoss = kTRUE;
Double_t *serrf = new Double_t[(input->GetSignal()->GetLast()) + 1];
Double_t *berrf = new Double_t[(input->GetSignal()->GetLast()) + 1];
do {
Double_t *toss = new Double_t[fgSystNames->GetSize()];
for (Int_t i = 0; i < fgSystNames->GetSize(); i++)
toss[i] = generator->Gaus(0, 1);
retoss = kFALSE;
for (Int_t channel = 0;
channel <= input->GetSignal()->GetLast();
channel++) {
serrf[channel] = 0;
berrf[channel] = 0;
for (Int_t bin = 0;
bin <((TVectorD *) (input->GetErrorOnSignal()->At(channel)))->GetNrows();
bin++) {
serrf[channel] += ((TVectorD *) (input->GetErrorOnSignal()->At(channel)))->operator[](bin) *
toss[fgSystNames->BinarySearch((TObjString*) (((TObjArray *) (input->GetErrorNames()->At(channel)))->At(bin)))];
berrf[channel] += ((TVectorD *) (input->GetErrorOnBackground()->At(channel)))->operator[](bin) *
toss[fgSystNames->BinarySearch((TObjString*) (((TObjArray *) (input->GetErrorNames()->At(channel)))->At(bin)))];
}
if ((serrf[channel] < -1.0) || (berrf[channel] < -0.9)) {
retoss = kTRUE;
continue;
}
}
delete[]toss;
} while (retoss);
output->SetOwner();
for (Int_t channel = 0; channel <= input->GetSignal()->GetLast();
channel++) {
TH1 *newsignal = (TH1*)(output->GetSignal()->At(channel));
TH1 *oldsignal = (TH1*)(input->GetSignal()->At(channel));
if(stat)
for(int i=1; i<=newsignal->GetNbinsX(); i++)
newsignal->SetBinContent(i,oldsignal->GetBinContent(i)+generator->Gaus(0,oldsignal->GetBinError(i)));
else
for(int i=1; i<=newsignal->GetNbinsX(); i++)
newsignal->SetBinContent(i,oldsignal->GetBinContent(i));
newsignal->Scale(1 + serrf[channel]);
newsignal->SetDirectory(0);
TH1 *newbackground = (TH1*)(output->GetBackground()->At(channel));
TH1 *oldbackground = (TH1*)(input->GetBackground()->At(channel));
if(stat)
for(int i=1; i<=newbackground->GetNbinsX(); i++)
newbackground->SetBinContent(i,oldbackground->GetBinContent(i)+generator->Gaus(0,oldbackground->GetBinError(i)));
else
for(int i=1; i<=newbackground->GetNbinsX(); i++)
newbackground->SetBinContent(i,oldbackground->GetBinContent(i));
newbackground->Scale(1 + berrf[channel]);
newbackground->SetDirectory(0);
}
delete[] serrf;
delete[] berrf;
return 1;
}
TConfidenceLevel *TLimit::ComputeLimit(TH1* s, TH1* b, TH1* d,
Int_t nmc, bool stat,
TRandom * generator)
{
TLimitDataSource* lds = new TLimitDataSource(s,b,d);
TConfidenceLevel* out = ComputeLimit(lds,nmc,stat,generator);
delete lds;
return out;
}
TConfidenceLevel *TLimit::ComputeLimit(TH1* s, TH1* b, TH1* d,
TVectorD* se, TVectorD* be, TObjArray* l,
Int_t nmc, bool stat,
TRandom * generator)
{
TLimitDataSource* lds = new TLimitDataSource(s,b,d,se,be,l);
TConfidenceLevel* out = ComputeLimit(lds,nmc,stat,generator);
delete lds;
return out;
}
TConfidenceLevel *TLimit::ComputeLimit(Double_t s, Double_t b, Int_t d,
Int_t nmc,
bool stat,
TRandom * generator)
{
TH1D* sh = new TH1D("__sh","__sh",1,0,2);
sh->Fill(1,s);
TH1D* bh = new TH1D("__bh","__bh",1,0,2);
bh->Fill(1,b);
TH1D* dh = new TH1D("__dh","__dh",1,0,2);
dh->Fill(1,d);
TLimitDataSource* lds = new TLimitDataSource(sh,bh,dh);
TConfidenceLevel* out = ComputeLimit(lds,nmc,stat,generator);
delete lds;
delete sh;
delete bh;
delete dh;
return out;
}
TConfidenceLevel *TLimit::ComputeLimit(Double_t s, Double_t b, Int_t d,
TVectorD* se, TVectorD* be, TObjArray* l,
Int_t nmc,
bool stat,
TRandom * generator)
{
TH1D* sh = new TH1D("__sh","__sh",1,0,2);
sh->Fill(1,s);
TH1D* bh = new TH1D("__bh","__bh",1,0,2);
bh->Fill(1,b);
TH1D* dh = new TH1D("__dh","__dh",1,0,2);
dh->Fill(1,d);
TLimitDataSource* lds = new TLimitDataSource(sh,bh,dh,se,be,l);
TConfidenceLevel* out = ComputeLimit(lds,nmc,stat,generator);
delete lds;
delete sh;
delete bh;
delete dh;
return out;
}
Double_t TLimit::LogLikelihood(Double_t s, Double_t b, Double_t b2, Double_t d)
{
return d*TMath::Log((s+b)/b2);
}