ROOT_528-00b_version: hist/hist/src/TConfidenceLevel.cxx Source File

00001 // @(#)root/hist:$Id: TConfidenceLevel.cxx 36150 2010-10-07 13:42:41Z rdm $
00002 // Author: Christophe.Delaere@cern.ch   21/08/2002
00003 
00004 ///////////////////////////////////////////////////////////////////////////
00005 //
00006 // TConfidenceLevel
00007 //
00008 // Class to compute 95% CL limits
00009 //
00010 ///////////////////////////////////////////////////////////////////////////
00011 
00012 /*************************************************************************
00013  * C.Delaere                                                             *
00014  * adapted from the mclimit code from Tom Junk                           *
00015  * see http://cern.ch/thomasj/searchlimits/ecl.html                      *
00016  *************************************************************************/
00017 
00018 #include "TConfidenceLevel.h"
00019 #include "TH1F.h"
00020 #include "TMath.h"
00021 #include "Riostream.h"
00022 
00023 ClassImp(TConfidenceLevel)
00024 
00025 Double_t const TConfidenceLevel::fgMCLM2S = 0.025;
00026 Double_t const TConfidenceLevel::fgMCLM1S = 0.16;
00027 Double_t const TConfidenceLevel::fgMCLMED = 0.5;
00028 Double_t const TConfidenceLevel::fgMCLP1S = 0.84;
00029 Double_t const TConfidenceLevel::fgMCLP2S = 0.975;
00030 // LHWG "one-sided" definition
00031 Double_t const TConfidenceLevel::fgMCL3S1S = 2.6998E-3;
00032 Double_t const TConfidenceLevel::fgMCL5S1S = 5.7330E-7;
00033 // the other definition (not chosen by the LHWG)
00034 Double_t const TConfidenceLevel::fgMCL3S2S = 1.349898E-3;
00035 Double_t const TConfidenceLevel::fgMCL5S2S = 2.866516E-7;
00036 
00037 
00038 //______________________________________________________________________________
00039 TConfidenceLevel::TConfidenceLevel()
00040 {
00041    // Default constructor
00042 
00043    fStot = 0;
00044    fBtot = 0;
00045    fDtot = 0;
00046    fTSD  = 0;
00047    fTSB  = 0;
00048    fTSS  = 0;
00049    fLRS  = 0;
00050    fLRB  = 0;
00051    fNMC  = 0;
00052    fNNMC = 0;
00053    fISS  = 0;
00054    fISB  = 0;
00055    fMCL3S = fgMCL3S1S;
00056    fMCL5S = fgMCL5S1S;
00057 }
00058 
00059 
00060 //______________________________________________________________________________
00061 TConfidenceLevel::TConfidenceLevel(Int_t mc, bool onesided)
00062 {
00063    // a constructor that fix some conventions:
00064    // mc is the number of Monte Carlo experiments
00065    // while onesided specifies if the intervals are one-sided or not.
00066 
00067    fStot = 0;
00068    fBtot = 0;
00069    fDtot = 0;
00070    fTSD  = 0;
00071    fTSB  = 0;
00072    fTSS  = 0;
00073    fLRS  = 0;
00074    fLRB  = 0;
00075    fNMC  = mc;
00076    fNNMC = mc;
00077    fISS  = new Int_t[mc];
00078    fISB  = new Int_t[mc];
00079    fMCL3S = onesided ? fgMCL3S1S : fgMCL3S2S;
00080    fMCL5S = onesided ? fgMCL5S1S : fgMCL5S2S;
00081 }
00082 
00083 
00084 //______________________________________________________________________________
00085 TConfidenceLevel::~TConfidenceLevel()
00086 {
00087    // The destructor
00088 
00089    if (fISS)
00090       delete[]fISS;
00091    if (fISB)
00092       delete[]fISB;
00093    if (fTSB)
00094       delete[]fTSB;
00095    if (fTSS)
00096       delete[]fTSS;
00097    if (fLRS)
00098       delete[]fLRS;
00099    if (fLRB)
00100       delete[]fLRB;
00101 }
00102 
00103 
00104 //______________________________________________________________________________
00105 Double_t TConfidenceLevel::GetExpectedStatistic_b(Int_t sigma) const
00106 {
00107    // Get the expected statistic value in the background only hypothesis
00108 
00109    switch (sigma) {
00110    case -2:
00111       return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
00112    case -1:
00113       return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
00114    case 0:
00115       return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
00116    case 1:
00117       return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
00118    case 2:
00119       return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
00120    default:
00121       return 0;
00122    }
00123 }
00124 
00125 
00126 //______________________________________________________________________________
00127 Double_t TConfidenceLevel::GetExpectedStatistic_sb(Int_t sigma) const
00128 {
00129    // Get the expected statistic value in the signal plus background hypothesis
00130 
00131    switch (sigma) {
00132    case -2:
00133       return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
00134    case -1:
00135       return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
00136    case 0:
00137       return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
00138    case 1:
00139       return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
00140    case 2:
00141       return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
00142    default:
00143       return 0;
00144    }
00145 }
00146 
00147 
00148 //______________________________________________________________________________
00149 Double_t TConfidenceLevel::CLb(bool use_sMC) const
00150 {
00151    // Get the Confidence Level for the background only
00152 
00153    Double_t result = 0;
00154    if (use_sMC) {
00155       for (Int_t i = 0; i < fNMC; i++)
00156          if (fTSS[fISS[i]] < fTSD)
00157             result += (1 / (fLRS[fISS[i]] * fNMC));
00158    } else {
00159       for (Int_t i = 0; i < fNMC; i++)
00160          if (fTSB[fISB[i]] < fTSD)
00161             result = (Double_t(i + 1)) / fNMC;
00162    }
00163    return result;
00164 }
00165 
00166 
00167 //______________________________________________________________________________
00168 Double_t TConfidenceLevel::CLsb(bool use_sMC) const
00169 {
00170    // Get the Confidence Level for the signal plus background hypothesis
00171 
00172    Double_t result = 0;
00173    if (use_sMC) {
00174       for (Int_t i = 0; i < fNMC; i++)
00175          if (fTSS[fISS[i]] <= fTSD)
00176             result = i / fNMC;
00177    } else {
00178       for (Int_t i = 0; i < fNMC; i++)
00179          if (fTSB[fISB[i]] <= fTSD)
00180             result += (fLRB[fISB[i]]) / fNMC;
00181    }
00182    return result;
00183 }
00184 
00185 
00186 //______________________________________________________________________________
00187 Double_t TConfidenceLevel::CLs(bool use_sMC) const
00188 {
00189    // Get the Confidence Level defined by CLs = CLsb/CLb.
00190    // This quantity is stable w.r.t. background fluctuations.
00191 
00192    Double_t clb = CLb(kFALSE);
00193    Double_t clsb = CLsb(use_sMC);
00194    if(clb==0) { cout << "Warning: clb = 0 !" << endl; return 0;}
00195    else return clsb/clb;
00196 }
00197 
00198 
00199 //______________________________________________________________________________
00200 Double_t TConfidenceLevel::GetExpectedCLsb_b(Int_t sigma) const
00201 {
00202    // Get the expected Confidence Level for the signal plus background hypothesis
00203    // if there is only background.
00204 
00205    Double_t result = 0;
00206    switch (sigma) {
00207    case -2:
00208       {
00209          for (Int_t i = 0; i < fNMC; i++)
00210             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
00211                result += fLRB[fISB[i]] / fNMC;
00212          return result;
00213       }
00214    case -1:
00215       {
00216          for (Int_t i = 0; i < fNMC; i++)
00217             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
00218                result += fLRB[fISB[i]] / fNMC;
00219          return result;
00220       }
00221    case 0:
00222       {
00223          for (Int_t i = 0; i < fNMC; i++)
00224             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
00225                result += fLRB[fISB[i]] / fNMC;
00226          return result;
00227       }
00228    case 1:
00229       {
00230          for (Int_t i = 0; i < fNMC; i++)
00231             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
00232                result += fLRB[fISB[i]] / fNMC;
00233          return result;
00234       }
00235    case 2:
00236       {
00237          for (Int_t i = 0; i < fNMC; i++)
00238             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
00239                result += fLRB[fISB[i]] / fNMC;
00240          return result;
00241       }
00242    default:
00243       return 0;
00244    }
00245 }
00246 
00247 
00248 //______________________________________________________________________________
00249 Double_t TConfidenceLevel::GetExpectedCLb_sb(Int_t sigma) const
00250 {
00251    // Get the expected Confidence Level for the background only
00252    // if there is signal and background.
00253 
00254    Double_t result = 0;
00255    switch (sigma) {
00256    case 2:
00257       {
00258          for (Int_t i = 0; i < fNMC; i++)
00259             if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
00260                result += fLRS[fISS[i]] / fNMC;
00261          return result;
00262       }
00263    case 1:
00264       {
00265          for (Int_t i = 0; i < fNMC; i++)
00266             if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
00267                result += fLRS[fISS[i]] / fNMC;
00268          return result;
00269       }
00270    case 0:
00271       {
00272          for (Int_t i = 0; i < fNMC; i++)
00273             if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
00274                result += fLRS[fISS[i]] / fNMC;
00275          return result;
00276       }
00277    case -1:
00278       {
00279          for (Int_t i = 0; i < fNMC; i++)
00280             if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
00281                result += fLRS[fISS[i]] / fNMC;
00282          return result;
00283       }
00284    case -2:
00285       {
00286          for (Int_t i = 0; i < fNMC; i++)
00287             if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
00288                result += fLRS[fISS[i]] / fNMC;
00289          return result;
00290       }
00291    default:
00292       return 0;
00293    }
00294 }
00295 
00296 
00297 //______________________________________________________________________________
00298 Double_t TConfidenceLevel::GetExpectedCLb_b(Int_t sigma) const
00299 {
00300    // Get the expected Confidence Level for the background only
00301    // if there is only background.
00302 
00303    Double_t result = 0;
00304    switch (sigma) {
00305    case 2:
00306       {
00307          for (Int_t i = 0; i < fNMC; i++)
00308             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
00309                result = (i + 1) / double (fNMC);
00310          return result;
00311       }
00312    case 1:
00313       {
00314          for (Int_t i = 0; i < fNMC; i++)
00315             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
00316                result = (i + 1) / double (fNMC);
00317          return result;
00318       }
00319    case 0:
00320       {
00321          for (Int_t i = 0; i < fNMC; i++)
00322             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
00323                result = (i + 1) / double (fNMC);
00324          return result;
00325       }
00326    case -1:
00327       {
00328          for (Int_t i = 0; i < fNMC; i++)
00329             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
00330                result = (i + 1) / double (fNMC);
00331          return result;
00332       }
00333    case -2:
00334       {
00335          for (Int_t i = 0; i < fNMC; i++)
00336             if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
00337                result = (i + 1) / double (fNMC);
00338          return result;
00339       }
00340    }
00341    return result;
00342 }
00343 
00344 
00345 //______________________________________________________________________________
00346 Double_t TConfidenceLevel::GetAverageCLsb() const
00347 {
00348    // Get average CLsb.
00349 
00350    Double_t result = 0;
00351    Double_t psumsb = 0;
00352    for (Int_t i = 0; i < fNMC; i++) {
00353       psumsb += fLRB[fISB[i]] / fNMC;
00354       result += psumsb / fNMC;
00355    }
00356    return result;
00357 }
00358 
00359 
00360 //______________________________________________________________________________
00361 Double_t TConfidenceLevel::GetAverageCLs() const
00362 {
00363    // Get average CLs.
00364 
00365    Double_t result = 0;
00366    Double_t psumsb = 0;
00367    for (Int_t i = 0; i < fNMC; i++) {
00368       psumsb += fLRB[fISB[i]] / fNMC;
00369       result += ((psumsb / fNMC) / ((i + 1) / fNMC));
00370    }
00371    return result;
00372 }
00373 
00374 
00375 //______________________________________________________________________________
00376 Double_t TConfidenceLevel::Get3sProbability() const
00377 {
00378    // Get 3s probability.
00379 
00380    Double_t result = 0;
00381    Double_t psumbs = 0;
00382    for (Int_t i = 0; i < fNMC; i++) {
00383       psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
00384       if (psumbs <= fMCL3S)
00385          result = i / fNMC;
00386    }
00387    return result;
00388 }
00389 
00390 
00391 //______________________________________________________________________________
00392 Double_t TConfidenceLevel::Get5sProbability() const
00393 {
00394    // Get 5s probability.
00395 
00396    Double_t result = 0;
00397    Double_t psumbs = 0;
00398    for (Int_t i = 0; i < fNMC; i++) {
00399       psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
00400       if (psumbs <= fMCL5S)
00401          result = i / fNMC;
00402    }
00403    return result;
00404 }
00405 
00406 
00407 //______________________________________________________________________________
00408 void  TConfidenceLevel::Draw(const Option_t*)
00409 {
00410    // Display sort of a "canonical" -2lnQ plot.
00411    // This results in a plot with 2 elements:
00412    // - The histogram of -2lnQ for background hypothesis (full)
00413    // - The histogram of -2lnQ for signal and background hypothesis (dashed)
00414    // The 2 histograms are respectively named b_hist and sb_hist.
00415 
00416    TH1F h("TConfidenceLevel_Draw","",50,0,0);
00417    Int_t i;
00418    for (i=0; i<fNMC; i++) {
00419       h.Fill(-2*(fTSB[i]-fStot));
00420       h.Fill(-2*(fTSS[i]-fStot));
00421    }
00422    TH1F* b_hist  = new TH1F("b_hist", "-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
00423    TH1F* sb_hist = new TH1F("sb_hist","-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
00424    for (i=0; i<fNMC; i++) {
00425       b_hist->Fill(-2*(fTSB[i]-fStot));
00426       sb_hist->Fill(-2*(fTSS[i]-fStot));
00427    }
00428    b_hist->Draw();
00429    sb_hist->Draw("Same");
00430    sb_hist->SetLineStyle(3);
00431 }
00432 
00433 
00434 //______________________________________________________________________________
00435 void  TConfidenceLevel::SetTSB(Double_t * in)
00436 {
00437    // Set the TSB.
00438    fTSB = in;
00439    TMath::Sort(fNNMC, fTSB, fISB, 0);
00440 }
00441 
00442 
00443 //______________________________________________________________________________
00444 void  TConfidenceLevel::SetTSS(Double_t * in)
00445 {
00446    // Set the TSS.
00447    fTSS = in;
00448    TMath::Sort(fNNMC, fTSS, fISS, 0);
00449 }