CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_1/src/PhysicsTools/RooStatsCms/src/BinomialInterval.cc

Go to the documentation of this file.
00001 #include <string>
00002 #include <cstdio>
00003 #include "Math/PdfFuncMathCore.h"
00004 #include "Math/QuantFuncMathCore.h"
00005 
00006 #if (defined (STANDALONE) or defined (__CINT__) )
00007 #include "BinomialInterval.h"
00008 
00009 ClassImp(BinomialInterval)
00010 #else
00011 #include "PhysicsTools/RooStatsCms/interface/BinomialInterval.h"
00012 #endif
00013 
00014 void BinomialInterval::init(const double alpha, const tail_type type) {
00015   alpha_     = alpha;
00016   type_      = type;
00017   alpha_min_ = type_ == equal_tailed ? alpha_/2 : alpha_;
00018   kappa_     = ROOT::Math::normal_quantile(1 - alpha/2, 1);
00019   kappa2_    = kappa_*kappa_;
00020 }
00021 
00022 bool BinomialInterval::contains(double p) {
00023   if (type_ == upper_tailed)
00024     return p <= upper_;
00025   else if (type_ == lower_tailed)
00026     return p >= lower_;
00027   else //if (type_ == equal_tailed)
00028     return p >= lower_ && p <= upper_;
00029 }
00030 
00031 double BinomialInterval::coverage_prob(const double p, const int trials) {
00032   double prob = 0;
00033 
00034   for (int X = 0; X <= trials; ++X) {
00035     calculate(X, trials);
00036 
00037     if (contains(p))
00038       prob += ROOT::Math::binomial_pdf(X, p, trials);
00039   }
00040   
00041   return prob;
00042 }
00043 
00044 void BinomialInterval::scan_rho(const int ntot, const int nrho, double* rho, double* prob) {
00045   for (int i = 0; i < nrho; ++i) {
00046     rho[i]  = double(i)/nrho;
00047     prob[i] = coverage_prob(rho[i], ntot);
00048   }
00049 }
00050 
00051 void BinomialInterval::scan_ntot(const double rho, const int ntot_min, const int ntot_max,
00052                                   double* ntot, double* prob) {
00053   for (int i = 0; i < ntot_max - ntot_min + 1; ++i) {
00054     int nt = i + ntot_min;
00055     ntot[i] = nt;
00056     prob[i] = coverage_prob(rho, nt);
00057   }
00058 }
00059 
00060 void BinomialInterval::dump(const int trials_min, const int trials_max) {
00061   const std::string fn = std::string("table.") + name() + std::string(".txt");
00062   FILE* fdump = fopen(fn.c_str(), "wt");
00063 
00064   for (int n = trials_min; n <= trials_max; ++n) {
00065     for (int X = 0; X <= n; X++) {
00066       calculate(X, n);
00067       fprintf(fdump, "%i %i %f %f\n", X, n, lower_, upper_);
00068     }
00069     fprintf(fdump, "\n");
00070   }
00071 
00072   fclose(fdump);
00073 }