#include "TConfidenceLevel.h"
#include "TH1F.h"
#include "TMath.h"
#include "Riostream.h"
ClassImp(TConfidenceLevel)
Double_t const TConfidenceLevel::fgMCLM2S = 0.025;
Double_t const TConfidenceLevel::fgMCLM1S = 0.16;
Double_t const TConfidenceLevel::fgMCLMED = 0.5;
Double_t const TConfidenceLevel::fgMCLP1S = 0.84;
Double_t const TConfidenceLevel::fgMCLP2S = 0.975;
Double_t const TConfidenceLevel::fgMCL3S1S = 2.6998E-3;
Double_t const TConfidenceLevel::fgMCL5S1S = 5.7330E-7;
Double_t const TConfidenceLevel::fgMCL3S2S = 1.349898E-3;
Double_t const TConfidenceLevel::fgMCL5S2S = 2.866516E-7;
TConfidenceLevel::TConfidenceLevel()
{
fStot = 0;
fBtot = 0;
fDtot = 0;
fTSD = 0;
fTSB = 0;
fTSS = 0;
fLRS = 0;
fLRB = 0;
fNMC = 0;
fNNMC = 0;
fISS = 0;
fISB = 0;
fMCL3S = fgMCL3S1S;
fMCL5S = fgMCL5S1S;
}
TConfidenceLevel::TConfidenceLevel(Int_t mc, bool onesided)
{
fStot = 0;
fBtot = 0;
fDtot = 0;
fTSD = 0;
fTSB = 0;
fTSS = 0;
fLRS = 0;
fLRB = 0;
fNMC = mc;
fNNMC = mc;
fISS = new Int_t[mc];
fISB = new Int_t[mc];
fMCL3S = onesided ? fgMCL3S1S : fgMCL3S2S;
fMCL5S = onesided ? fgMCL5S1S : fgMCL5S2S;
}
TConfidenceLevel::~TConfidenceLevel()
{
if (fISS)
delete[]fISS;
if (fISB)
delete[]fISB;
if (fTSB)
delete[]fTSB;
if (fTSS)
delete[]fTSS;
if (fLRS)
delete[]fLRS;
if (fLRB)
delete[]fLRB;
}
Double_t TConfidenceLevel::GetExpectedStatistic_b(Int_t sigma) const
{
switch (sigma) {
case -2:
return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
case -1:
return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
case 0:
return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
case 1:
return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
case 2:
return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
default:
return 0;
}
}
Double_t TConfidenceLevel::GetExpectedStatistic_sb(Int_t sigma) const
{
switch (sigma) {
case -2:
return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
case -1:
return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
case 0:
return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
case 1:
return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
case 2:
return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
default:
return 0;
}
}
Double_t TConfidenceLevel::CLb(bool use_sMC) const
{
Double_t result = 0;
if (use_sMC) {
for (Int_t i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] < fTSD)
result += (1 / (fLRS[fISS[i]] * fNMC));
} else {
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] < fTSD)
result = (Double_t(i + 1)) / fNMC;
}
return result;
}
Double_t TConfidenceLevel::CLsb(bool use_sMC) const
{
Double_t result = 0;
if (use_sMC) {
for (Int_t i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSD)
result = i / fNMC;
} else {
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSD)
result += (fLRB[fISB[i]]) / fNMC;
}
return result;
}
Double_t TConfidenceLevel::CLs(bool use_sMC) const
{
Double_t clb = CLb(kFALSE);
Double_t clsb = CLsb(use_sMC);
if(clb==0) { cout << "Warning: clb = 0 !" << endl; return 0;}
else return clsb/clb;
}
Double_t TConfidenceLevel::GetExpectedCLsb_b(Int_t sigma) const
{
Double_t result = 0;
switch (sigma) {
case -2:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
case -1:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
case 0:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
case 1:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
case 2:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
default:
return 0;
}
}
Double_t TConfidenceLevel::GetExpectedCLb_sb(Int_t sigma) const
{
Double_t result = 0;
switch (sigma) {
case 2:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
case 1:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
case 0:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
case -1:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
case -2:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
default:
return 0;
}
}
Double_t TConfidenceLevel::GetExpectedCLb_b(Int_t sigma) const
{
Double_t result = 0;
switch (sigma) {
case 2:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
result = (i + 1) / double (fNMC);
return result;
}
case 1:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
result = (i + 1) / double (fNMC);
return result;
}
case 0:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
result = (i + 1) / double (fNMC);
return result;
}
case -1:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
result = (i + 1) / double (fNMC);
return result;
}
case -2:
{
for (Int_t i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
result = (i + 1) / double (fNMC);
return result;
}
}
return result;
}
Double_t TConfidenceLevel::GetAverageCLsb() const
{
Double_t result = 0;
Double_t psumsb = 0;
for (Int_t i = 0; i < fNMC; i++) {
psumsb += fLRB[fISB[i]] / fNMC;
result += psumsb / fNMC;
}
return result;
}
Double_t TConfidenceLevel::GetAverageCLs() const
{
Double_t result = 0;
Double_t psumsb = 0;
for (Int_t i = 0; i < fNMC; i++) {
psumsb += fLRB[fISB[i]] / fNMC;
result += ((psumsb / fNMC) / ((i + 1) / fNMC));
}
return result;
}
Double_t TConfidenceLevel::Get3sProbability() const
{
Double_t result = 0;
Double_t psumbs = 0;
for (Int_t i = 0; i < fNMC; i++) {
psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
if (psumbs <= fMCL3S)
result = i / fNMC;
}
return result;
}
Double_t TConfidenceLevel::Get5sProbability() const
{
Double_t result = 0;
Double_t psumbs = 0;
for (Int_t i = 0; i < fNMC; i++) {
psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
if (psumbs <= fMCL5S)
result = i / fNMC;
}
return result;
}
void TConfidenceLevel::Draw(const Option_t*)
{
TH1F h("TConfidenceLevel_Draw","",50,0,0);
Int_t i;
for (i=0; i<fNMC; i++) {
h.Fill(-2*(fTSB[i]-fStot));
h.Fill(-2*(fTSS[i]-fStot));
}
TH1F* b_hist = new TH1F("b_hist", "-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
TH1F* sb_hist = new TH1F("sb_hist","-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
for (i=0; i<fNMC; i++) {
b_hist->Fill(-2*(fTSB[i]-fStot));
sb_hist->Fill(-2*(fTSS[i]-fStot));
}
b_hist->Draw();
sb_hist->Draw("Same");
sb_hist->SetLineStyle(3);
}
void TConfidenceLevel::SetTSB(Double_t * in)
{
fTSB = in;
TMath::Sort(fNNMC, fTSB, fISB, 0);
}
void TConfidenceLevel::SetTSS(Double_t * in)
{
fTSS = in;
TMath::Sort(fNNMC, fTSS, fISS, 0);
}