#include <Asymptotic.h>

Inheritance diagram for Asymptotic:

Public Member Functions
virtual void	applyDefaultOptions ()

virtual void	applyOptions (const boost::program_options::variables_map &vm)

	Asymptotic ()

float	findExpectedLimitFromCrossing (RooAbsReal &nll, RooRealVar *r, double rMin, double rMax, double nll0, double quantile)

virtual const std::string &	name () const

virtual bool	run (RooWorkspace w, RooStats::ModelConfig mc_s, RooStats::ModelConfig mc_b, RooAbsData &data, double &limit, double &limitErr, const double hint)

virtual bool	runLimit (RooWorkspace w, RooStats::ModelConfig mc_s, RooStats::ModelConfig mc_b, RooAbsData &data, double &limit, double &limitErr, const double hint)

std::vector< std::pair< float, float > >	runLimitExpected (RooWorkspace w, RooStats::ModelConfig mc_s, RooStats::ModelConfig mc_b, RooAbsData &data, double &limit, double &limitErr, const double hint)

Public Member Functions inherited from LimitAlgo
	LimitAlgo ()

	LimitAlgo (const char *desc)

const boost::program_options::options_description &	options () const

virtual void	setNToys (const int)

virtual void	setToyNumber (const int)

Private Member Functions
RooAbsData *	asimovDataset (RooWorkspace w, RooStats::ModelConfig mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data)

double	getCLs (RooRealVar &r, double rVal, bool getAlsoExpected=false, double limit=0, double limitErr=0)

Private Attributes
utils::CheapValueSnapshot	fitFixA_

utils::CheapValueSnapshot	fitFixD_

utils::CheapValueSnapshot	fitFreeA_

utils::CheapValueSnapshot	fitFreeD_

bool	hasFloatParams_

double	minNllA_

double	minNllD_

std::auto_ptr< RooAbsReal >	nllA_

std::auto_ptr< RooAbsReal >	nllD_

std::auto_ptr< RooArgSet >	params_

RooArgSet	snapGlobalObsAsimov

RooArgSet	snapGlobalObsData

Static Private Attributes
static std::string	minimizerAlgo_ = "Minuit2"

static int	minimizerStrategy_ = 0

static float	minimizerTolerance_ = 0.01

static std::string	minosAlgo_ = "stepping"

static bool	newExpected_ = false

static bool	noFitAsimov_ = false

static bool	picky_ = false

static bool	qtilde_ = true

static double	rAbsAccuracy_ = 0.0005

static double	rRelAccuracy_ = 0.005

static double	rValue_ = 1.0

static std::string	what_ = "both"

Additional Inherited Members
Protected Attributes inherited from LimitAlgo
boost::program_options::options_description	options_

Detailed Description

CLs asymptotic limits

Author: Giovanni Petrucciani (UCSD)

Definition at line 19 of file Asymptotic.h.

Constructor & Destructor Documentation

Asymptotic::Asymptotic ( )

Definition at line 37 of file Asymptotic.cc.

References minimizerAlgo_, minimizerStrategy_, minimizerTolerance_, minosAlgo_, LimitAlgo::options_, qtilde_, rAbsAccuracy_, rRelAccuracy_, rValue_, and what_.

                        : 
 LimitAlgo("Asymptotic specific options") {
     options_.add_options()
         ("rAbsAcc", boost::program_options::value<double>(&rAbsAccuracy_)->default_value(rAbsAccuracy_), "Absolute accuracy on r to reach to terminate the scan")
         ("rRelAcc", boost::program_options::value<double>(&rRelAccuracy_)->default_value(rRelAccuracy_), "Relative accuracy on r to reach to terminate the scan")
         ("run", boost::program_options::value<std::string>(&what_)->default_value(what_), "What to run: both (default), observed, expected, blind.")
         ("singlePoint",  boost::program_options::value<double>(&rValue_),  "Just compute CLs for the given value of r")
         ("minimizerAlgo",      boost::program_options::value<std::string>(&minimizerAlgo_)->default_value(minimizerAlgo_), "Choice of minimizer used for profiling (Minuit vs Minuit2)")
         ("minimizerTolerance", boost::program_options::value<float>(&minimizerTolerance_)->default_value(minimizerTolerance_),  "Tolerance for minimizer used for profiling")
         ("minimizerStrategy",  boost::program_options::value<int>(&minimizerStrategy_)->default_value(minimizerStrategy_),      "Stragegy for minimizer")
         ("qtilde", boost::program_options::value<bool>(&qtilde_)->default_value(qtilde_),  "Allow only non-negative signal strengths (default is true).")
         ("picky", "Abort on fit failures")
         ("noFitAsimov", "Use the pre-fit asimov dataset")
         ("newExpected", "Use the new formula for expected limits")
         ("minosAlgo", boost::program_options::value<std::string>(&minosAlgo_)->default_value(minosAlgo_), "Algorithm to use to get the median expected limit: 'minos' (fastest), 'bisection', 'stepping' (default, most robust)")
     ;
 }

Member Function Documentation

void Asymptotic::applyDefaultOptions ( )

virtual

Reimplemented from LimitAlgo.

Definition at line 70 of file Asymptotic.cc.

References noFitAsimov_, and what_.

                                      { 
     what_ = "observed"; noFitAsimov_ = true; // faster
 }

void Asymptotic::applyOptions ( const boost::program_options::variables_map & vm )

virtual

Reimplemented from LimitAlgo.

Definition at line 55 of file Asymptotic.cc.

References gather_cfg::cout, newExpected_, noFitAsimov_, picky_, and what_.

                                                                          {
     if (vm.count("singlePoint") && !vm["singlePoint"].defaulted()) {
         if (!vm["run"].defaulted()) throw std::invalid_argument("Asymptotic: when using --singlePoint you can't use --run (at least for now)");
         what_ = "singlePoint";
     } else {
         if (what_ != "observed" && what_ != "expected" && what_ != "both" && what_ != "blind") 
             throw std::invalid_argument("Asymptotic: option 'run' can only be 'observed', 'expected' or 'both' (the default) or 'blind' (a-priori expected)");
     }
     picky_ = vm.count("picky");
     noFitAsimov_ = vm.count("noFitAsimov");
     newExpected_ = vm.count("newExpected");
     if (what_ == "blind") { what_ = "expected"; noFitAsimov_ = true; } 
     if (noFitAsimov_) std::cout << "Will use a-priori expected background instead of a-posteriori one." << std::endl; 
 }

RooAbsData * Asymptotic::asimovDataset	(	RooWorkspace *	w,
		RooStats::ModelConfig *	mc_s,
		RooStats::ModelConfig *	mc_b,
		RooAbsData &	data
	)

private

Definition at line 542 of file Asymptotic.cc.

References asimovutils::asimovDatasetNominal(), asimovutils::asimovDatasetWithFit(), noFitAsimov_, utils::setAllConstant(), snapGlobalObsAsimov, snapGlobalObsData, and withSystematics.

Referenced by runLimit(), and runLimitExpected().

                                                                                                                               {
     // Do this only once
     if (w->data("_Asymptotic_asimovDataset_") != 0) {
         return w->data("_Asymptotic_asimovDataset_");
     }
     // snapshot data global observables
     RooArgSet gobs;
     if (withSystematics && mc_s->GetGlobalObservables()) {
         gobs.add(*mc_s->GetGlobalObservables());
         snapGlobalObsData.removeAll();
         utils::setAllConstant(gobs, true);
         gobs.snapshot(snapGlobalObsData);
     }
     // get asimov dataset and global observables
     RooAbsData *asimovData = (noFitAsimov_  ? asimovutils::asimovDatasetNominal(mc_s, 0.0, verbose) :
                                               asimovutils::asimovDatasetWithFit(mc_s, data, snapGlobalObsAsimov, 0.0, verbose));
     asimovData->SetName("_Asymptotic_asimovDataset_");
     w->import(*asimovData); // I'm assuming the Workspace takes ownership. Might be false.
     delete asimovData;      //  ^^^^^^^^----- now assuming that the workspace clones.
     return w->data("_Asymptotic_asimovDataset_");
 }

float Asymptotic::findExpectedLimitFromCrossing	(	RooAbsReal &	nll,
		RooRealVar *	r,
		double	rMin,
		double	rMax,
		double	nll0,
		double	quantile
	)

Definition at line 349 of file Asymptotic.cc.

References cl, CloseCoutSentry::clear(), CascadeMinimizer::Constrained, CascadeMinimizer::improve(), max(), min, CascadeMinimizer::minimize(), CascadeMinimizer::minimizer(), minimizerAlgo_, minimizerStrategy_, minimizerTolerance_, minosAlgo_, N, convertSQLiteXML::ok, picky_, funct::pow(), rAbsAccuracy_, rRelAccuracy_, CascadeMinimizer::setErrorLevel(), CascadeMinimizer::setStrategy(), dtDQMClient_cfg::threshold, CascadeMinimizer::Unconstrained, and detailsBasic3DVector::y.

Referenced by runLimitExpected().

                                                                                                                                  {
     // EQ 37 of CMS NOTE 2011-005:
     //   mu_N = sigma * ( normal_quantile_c( (1-cl) * normal_cdf(N) ) + N )
     // --> (mu_N/sigma) = N + normal_quantile_c( (1-cl) * clb )
     // but qmu = (mu_N/sigma)^2
     // --> qmu = [ N + normal_quantile_c( (1-cl)*CLb ) ]^2
     // remember that qmu = 2*nll
     double N = ROOT::Math::normal_quantile(clb, 1.0);
     double errorlevel = 0.5 * pow(N+ROOT::Math::normal_quantile_c(clb*(1-cl),1.0), 2);
     int minosStat = -1;
     if (minosAlgo_ == "minos") {
         double rMax0 = r->getMax();
         // Have to repeat the fit, but I'm already at the minimum
         CascadeMinimizer minim(nll, CascadeMinimizer::Unconstrained, r);
         minim.setStrategy(minimizerStrategy_);
         minim.setErrorLevel(errorlevel); 
         CloseCoutSentry sentry(verbose < 3);
         minim.minimize(verbose-2);
         sentry.clear();
         for (int tries = 0; tries < 3; ++tries) {
             minosStat = minim.minimizer().minos(RooArgSet(*r));
             if (minosStat != -1) {
                 while ((minosStat != -1) && (r->getVal()+r->getAsymErrorHi())/r->getMax() > 0.9) {
                     if (r->getMax() >= 100*rMax0) { minosStat = -1; break; }
                     r->setMax(2*r->getMax());
                     CascadeMinimizer minim2(nll, CascadeMinimizer::Unconstrained, r);
                     minim2.setStrategy(minimizerStrategy_);
                     minim2.setErrorLevel(errorlevel); 
                     minim2.minimize(verbose-2);
                     minosStat = minim2.minimizer().minos(RooArgSet(*r));
                 }
                 break;
             }
             minim.setStrategy(2);
             if (tries == 1) { 
                 if (minimizerAlgo_.find("Minuit2") != std::string::npos) {
                     minim.minimizer().minimize("Minuit","minimize");
                 } else {
                     minim.minimizer().minimize("Minuit2","minmize");
                 }
             }
         }
         if (minosStat != -1) return r->getVal()+r->getAsymErrorHi();
     } else {
         double threshold = nll0 + errorlevel;
         double rCross = 0.5*(rMin+rMax), rErr = 0.5*(rMax-rMin);
         r->setVal(rCross); r->setConstant(true);
         CascadeMinimizer minim2(nll, CascadeMinimizer::Constrained);
         minim2.setStrategy(minimizerStrategy_);
         if (minosAlgo_ == "bisection") {
             if (verbose > 1) printf("Will search for NLL crossing by bisection\n");
             while (rErr > std::max(rRelAccuracy_*rCross, rAbsAccuracy_)) {
                 if (rCross >= r->getMax()) r->setMax(rCross*1.1);
                 r->setVal(rCross);
                 bool ok = true;
                 { 
                     CloseCoutSentry sentry2(verbose < 3);
                     ok = minim2.improve(verbose-2);
                 }
                 if (!ok && picky_) break; else minosStat = 0;
                 double here = nll.getVal();
                 if (verbose > 1) printf("At %s = %f:\tdelta(nll) = %.5f\n", r->GetName(), rCross, here-nll0);
                 if (fabs(here - threshold) < 0.05*minimizerTolerance_) break;
                 if (here < threshold) rMin = rCross; else rMax = rCross;
                 rCross = 0.5*(rMin+rMax); rErr = 0.5*(rMax-rMin);
             } 
         } else if (minosAlgo_ == "stepping") {
             if (verbose > 1) printf("Will search for NLL crossing by stepping.\n");
             rCross = 0.05 * rMax; rErr = rMax; 
             double stride = rCross; bool overstepped = false;
             while (rErr > std::max(rRelAccuracy_*rCross, rAbsAccuracy_)) {
                 if (rCross >= r->getMax()) r->setMax(rCross*1.1);
                 double there = nll.getVal();
                 r->setVal(rCross);
                 bool ok = true;
                 { 
                     CloseCoutSentry sentry2(verbose < 3);
                     ok = minim2.improve(verbose-2);
                 }
                 if (!ok && picky_) break; else minosStat = 0;
                 double here = nll.getVal();
                 if (verbose > 1) printf("At %s = %f:\tdelta(nll) = %.5f\n", r->GetName(), rCross, here-nll0);
                 if (fabs(here - threshold) < 0.05*minimizerTolerance_) break;
                 if (here < threshold) { 
                     if ((threshold-here) < 0.5*fabs(threshold-there)) stride *= 0.5;
                     rCross += stride; 
                 } else { 
                     stride *= 0.5; overstepped = true;
                     rCross -= stride;
                 }
                 if (overstepped) rErr = stride;
             }
         } else if (minosAlgo_ == "new") {
             if (verbose > 1) printf("Will search for NLL crossing with new algorithm.\n");
             // 
             // Let X(x,y) = (x-a*y)^2 / s^2 + y^2    be the chi-square in case of correlations
             // then yhat(x) = a*x / (a^2 + s^2)
             // and X(x, yhat(x)) = x^2 / (a^2 + s^2)
             // For an unprofiled step
             //    X(x+dx, yhat(x)) - X(x,yhat(x))     = dx^2 / s^2         + 2 * x * dx / (a^2 + s^2)
             // For a profiled step  
             //    X(x+dx, yhat(x+dx)) - X(x,yhat(x))  = dx^2 / (a^2 + s^2) + 2 * x * dx / (a^2 + s^2)
             // So,
             //     dX_prof - dX_unprof = dx^2 * a^2 / (s^2 * (a^2 + s^2) )
             // The idea is then to take this approximation
             //     X_approx(x)  =  X(x, y(x1)) - k * (x-x1)^2 
             //           with k = [ X(x1, y1(x0)) - X(y1, y1(x1) ] / (x1-x0)^2
             double r_0   = rMin; 
             r->setVal(rMin);
             double nll_0 = nll.getVal(); 
             double rMax0 = rMax*100;
             double kappa = 0;
             // part 1: try to bracket the crossing between two points that have profiled nll above & below threshold
             double rStep = 0.05 * (rMax - r_0);
             double r_1 = r_0, nll_1 = nll_0;
             do {
                 r_1 += rStep; 
                 if (r_1 >= r->getMax()) r->setMax(r_1*1.1);
                 r->setVal(r_1);
                 nll_1 = nll.getVal();
                 // we profile if the NLL changed by more than 0.5, or if we got above threshold
                 bool binNLLchange = (nll_1 < threshold && nll_1 - nll_0 > 0.5);
                 bool aboveThresh  = (nll_1 > threshold + kappa*std::pow(r_1-r_0,2));
                 if (binNLLchange || aboveThresh) { 
                     if (verbose > 1) printf("At %s = %f:\tdelta(nll unprof) = %.5f\t                         \tkappa=%.5f\n", r->GetName(), r_1, nll_1-nll0, kappa);
                     { 
                         CloseCoutSentry sentry2(verbose < 3);
                         bool ok = minim2.improve(verbose-2);
                         if (!ok && picky_) return std::numeric_limits<float>::quiet_NaN();
                     }
                     double nll_1_prof = nll.getVal();
                     kappa = (nll_1 - nll_1_prof) / std::pow(r_1 - r_0,2);
                     if (verbose > 1) printf("At %s = %f:\tdelta(nll unprof) = %.5f\tdelta(nll prof) = %.5f\tkappa=%.5f\n", r->GetName(), r_1, nll_1-nll0, nll.getVal()-nll0, kappa);
                     if (nll_1_prof > threshold) { 
                         nll_1 = nll_1_prof; 
                         break; 
                     } else {
                         r_0 = r_1; 
                         nll_0 = nll_1_prof;
                         if (aboveThresh) rStep *= 2;
                     }
                 } else {
                     if (verbose > 1) printf("At %s = %f:\tdelta(nll unprof) = %.5f\t                         \tkappa=%.5f\n", r->GetName(), r_1, nll_1-nll0, kappa);
                 }
                 if (r_1 > rMax0) return std::numeric_limits<float>::quiet_NaN();
             } while (true);
             // now crossing is bracketed, do bisection
             if (verbose > 1) printf("At %s = %f:\t                         \tdelta(nll prof) = %.5f\tkappa=%.5f\n", r->GetName(), r_0, nll_0-nll0, kappa);
             if (verbose > 1) printf("At %s = %f:\t                         \tdelta(nll prof) = %.5f\tkappa=%.5f\n", r->GetName(), r_1, nll_1-nll0, kappa);
             minosStat = 0;
             do {
                // LOOP PRECONDITIONS:
                //   - r_0 and r_1 have profiled nll values on the two sides of the threshold
                //   - nuisance parameters have been profiled at r_1 
                double rEps = 0.2*std::max(rRelAccuracy_*r_1, rAbsAccuracy_);
                // bisection loop to find point with the right nll_approx
                double r_lo = std::min(r_0,r_1), r_hi = std::max(r_1,r_0);
                while (r_hi - r_lo > rEps) {
                    double r_2 = 0.5*(r_hi+r_lo); 
                    r->setVal(r_2);
                    double y0 = nll.getVal(), y = y0 - kappa*std::pow(r_2-r_1,2);
                    if (verbose > 1) printf("At %s = %f:\tdelta(nll unprof) = %.5f\tdelta(nll appr) = %.5f\tkappa=%.5f\n", r->GetName(), r_2, y0-nll0, y-nll0, kappa);
                    if (y < threshold) { r_lo = r_2; } else { r_hi = r_2; }
                } 
                // profile at that point
                rCross = r->getVal(); 
                double nll_unprof = nll.getVal();
                bool ok = true;
                { 
                    CloseCoutSentry sentry2(verbose < 3);
                    ok = minim2.improve(verbose-2);
                }
                if (!ok && picky_) return std::numeric_limits<float>::quiet_NaN();
                double nll_prof = nll.getVal();
                if (verbose > 1) printf("At %s = %f:\tdelta(nll unprof) = %.5f\tdelta(nll prof) = %.5f\tdelta(nll appr) = %.5f\n", r->GetName(), rCross, nll_unprof-nll0, nll_prof-nll0, nll_unprof-nll0 - kappa*std::pow(rCross-r_1,2));
                if (fabs(nll_prof - threshold) < 0.1*minimizerTolerance_) { break; }
                // not yet bang on, so update r_0, kappa
                kappa = (nll_unprof - nll_prof)/std::pow(rCross-r_1,2);
                // (r0 or r1) --> r0, and rCross --> r1;  
                if ((nll_prof < threshold) == (nll_0 < threshold)) { // if rCross is on the same side of r_0
                    r_0 = r_1;   
                    nll_0 = nll_1; 
                } else {
                    // stay with r_0 as is
                }
                r_1 = rCross; nll_1 = nll_prof;
             } while (fabs(r_1-r_0) > std::max(rRelAccuracy_*rCross, rAbsAccuracy_));
         }
         if (minosStat != -1) return rCross;
     }
     return std::numeric_limits<float>::quiet_NaN();
 }

double Asymptotic::getCLs	(	RooRealVar &	r,
		double	rVal,
		bool	getAlsoExpected = `false`,
		double *	limit = `0`,
		double *	limitErr = `0`
	)

private

Definition at line 203 of file Asymptotic.cc.

References CloseCoutSentry::clear(), Combine::commitPoint(), CascadeMinimizer::Constrained, utils::CheapValueSnapshot::empty(), fitFixA_, fitFixD_, fitFreeA_, fitFreeD_, hasFloatParams_, CascadeMinimizer::improve(), minimizerStrategy_, minNllA_, minNllD_, N, nllA_, nllD_, params_, picky_, utils::CheapValueSnapshot::Print(), qtilde_, utils::CheapValueSnapshot::readFrom(), CascadeMinimizer::setStrategy(), snapGlobalObsAsimov, snapGlobalObsData, and mathSSE::sqrt().

Referenced by runLimit().

                                                                                                            {
   r.setMax(1.1 * rVal);
   r.setConstant(true);
 
   CloseCoutSentry sentry(verbose < 3);
 
   CascadeMinimizer minimD(*nllD_, CascadeMinimizer::Constrained, &r);
   minimD.setStrategy(minimizerStrategy_);  
 
   (!fitFixD_.empty() ? fitFixD_ : fitFreeD_).writeTo(*params_);
   *params_ = snapGlobalObsData;
   r.setVal(rVal);
   r.setConstant(true);
   if (hasFloatParams_) {
       if (!minimD.improve(verbose-2) && picky_) return -999;
       fitFixD_.readFrom(*params_);
       if (verbose >= 2) fitFixD_.Print("V");
   }
   double qmu = 2*(nllD_->getVal() - minNllD_); if (qmu < 0) qmu = 0;
 
   CascadeMinimizer minimA(*nllA_, CascadeMinimizer::Constrained, &r);
   minimA.setStrategy(minimizerStrategy_); 
 
   (!fitFixA_.empty() ? fitFixA_ : fitFreeA_).writeTo(*params_);
   *params_ = snapGlobalObsAsimov;
   r.setVal(rVal);
   r.setConstant(true);
   if (hasFloatParams_) {
       if (!minimA.improve(verbose-2) && picky_) return -999;
       fitFixA_.readFrom(*params_);
       if (verbose >= 2) fitFixA_.Print("V");
   }
   double qA  = 2*(nllA_->getVal() - minNllA_); if (qA < 0) qA = 0;
 
   double CLsb = ROOT::Math::normal_cdf_c(sqrt(qmu));
   double CLb  = ROOT::Math::normal_cdf(sqrt(qA)-sqrt(qmu));
   if (qtilde_ && qmu > qA) {
     // In this region, things are tricky
     double mos = sqrt(qA); // mu/sigma
     CLsb = ROOT::Math::normal_cdf_c( (qmu + qA)/(2*mos) );
     CLb  = ROOT::Math::normal_cdf_c( (qmu - qA)/(2*mos) );
   }
   double CLs  = (CLb == 0 ? 0 : CLsb/CLb);
   sentry.clear();
   if (verbose > 0) printf("At %s = %f:\tq_mu = %.5f\tq_A  = %.5f\tCLsb = %7.5f\tCLb  = %7.5f\tCLs  = %7.5f\n", r.GetName(), rVal, qmu, qA, CLsb, CLb, CLs);
 
   if (getAlsoExpected) {
     const double quantiles[5] = { 0.025, 0.16, 0.50, 0.84, 0.975 };
     for (int iq = 0; iq < 5; ++iq) {
         double N = ROOT::Math::normal_quantile(quantiles[iq], 1.0);
         double clb = quantiles[iq];
         double clsplusb = ROOT::Math::normal_cdf_c( sqrt(qA) - N, 1.);
         *limit = (clb != 0 ? clsplusb/clb : 0); *limitErr = 0;
         Combine::commitPoint(true, quantiles[iq]);
         if (verbose > 0) printf("Expected %4.1f%%: CLsb = %.5f  CLb = %.5f   CLs = %.5f\n", quantiles[iq]*100, clsplusb, clb, clsplusb/clb);
     }
   }
   return CLs; 
 }   

virtual const std::string& Asymptotic::name ( void ) const

inlinevirtual

Implements LimitAlgo.

Definition at line 31 of file Asymptotic.h.

Referenced by BeautifulSoup.Tag::_invert().

31 { static std::string name_ = "Asymptotic"; return name_; }

bool Asymptotic::run	(	RooWorkspace *	w,
		RooStats::ModelConfig *	mc_s,
		RooStats::ModelConfig *	mc_b,
		RooAbsData &	data,
		double &	limit,
		double &	limitErr,
		const double *	hint
	)

virtual

Implements LimitAlgo.

Definition at line 74 of file Asymptotic.cc.

References gather_cfg::cout, DEBUG, minimizerAlgo_, minimizerStrategy_, minimizerTolerance_, run_regression::ret, runLimit(), runLimitExpected(), rValue_, dqm::qstatus::WARNING, and what_.

                                                                                                                                                                  {
     RooFitGlobalKillSentry silence(verbose <= 1 ? RooFit::WARNING : RooFit::DEBUG);
     ProfileLikelihood::MinimizerSentry minimizerConfig(minimizerAlgo_, minimizerTolerance_);
     if (verbose > 0) std::cout << "Will compute " << what_ << " limit(s) using minimizer " << minimizerAlgo_ 
                         << " with strategy " << minimizerStrategy_ << " and tolerance " << minimizerTolerance_ << std::endl;
     
     bool ret = false; 
     std::vector<std::pair<float,float> > expected;
     if (what_ == "both" || what_ == "expected") expected = runLimitExpected(w, mc_s, mc_b, data, limit, limitErr, hint);
     if (what_ != "expected") ret = runLimit(w, mc_s, mc_b, data, limit, limitErr, hint);
 
     if (verbose >= 0) {
         const char *rname = mc_s->GetParametersOfInterest()->first()->GetName();
         std::cout << "\n -- Asymptotic -- " << "\n";
         if (what_ == "singlePoint") {
             printf("Observed CLs for %s = %.1f: %6.4f \n", rname, rValue_, limit);
         } else if (ret && what_ != "expected") {
             printf("Observed Limit: %s < %6.4f\n", rname, limit);
         }
         for (std::vector<std::pair<float,float> >::const_iterator it = expected.begin(), ed = expected.end(); it != ed; ++it) {
             printf("Expected %4.1f%%: %s < %6.4f\n", it->first*100, rname, it->second);
         }
         std::cout << std::endl;
     }
 
     // note that for expected we have to return FALSE even if we succeed because otherwise it goes into the observed limit as well
     return ret;
 }

bool Asymptotic::runLimit	(	RooWorkspace *	w,
		RooStats::ModelConfig *	mc_s,
		RooStats::ModelConfig *	mc_b,
		RooAbsData &	data,
		double &	limit,
		double &	limitErr,
		const double *	hint
	)

virtual

Definition at line 103 of file Asymptotic.cc.

References a, asimovDataset(), dtNoiseDBValidation_cfg::cerr, cl, createBeamHaloJobs::constraints, gather_cfg::cout, funct::false, fitFreeA_, fitFreeD_, getCLs(), hasFloatParams_, MessageLogger_cff::limit, create_public_lumi_plots::log, max(), CascadeMinimizer::minimize(), minimizerStrategy_, minNllA_, minNllD_, nllA_, nllD_, params_, utils::CheapValueSnapshot::Print(), qtilde_, alignCSCRings::r, rAbsAccuracy_, utils::CheapValueSnapshot::readFrom(), rRelAccuracy_, rValue_, CascadeMinimizer::setStrategy(), snapGlobalObsAsimov, snapGlobalObsData, CascadeMinimizer::Unconstrained, what_, withSystematics, and utils::CheapValueSnapshot::writeTo().

Referenced by run().

                                                                                                                                                                       {
   RooRealVar *r = dynamic_cast<RooRealVar *>(mc_s->GetParametersOfInterest()->first()); 
   w->loadSnapshot("clean");
   RooAbsData &asimov = *asimovDataset(w, mc_s, mc_b, data);
 
   r->setConstant(false);
   r->setVal(0.1*r->getMax());
   r->setMin(qtilde_ ? 0 : -r->getMax());
  
   if (params_.get() == 0) params_.reset(mc_s->GetPdf()->getParameters(data));
 
   hasFloatParams_ = false;
   std::auto_ptr<TIterator> itparam(params_->createIterator());
   for (RooAbsArg *a = (RooAbsArg *) itparam->Next(); a != 0; a = (RooAbsArg *) itparam->Next()) {
       RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);
       if ( rrv != 0 && rrv != r && rrv->isConstant() == false ) { hasFloatParams_ = true; break; }
   }
     
 
   RooArgSet constraints; if (withSystematics) constraints.add(*mc_s->GetNuisanceParameters());
   nllD_.reset(mc_s->GetPdf()->createNLL(data,   RooFit::Constrain(constraints)));
   nllA_.reset(mc_s->GetPdf()->createNLL(asimov, RooFit::Constrain(constraints)));
 
   if (verbose > 0) std::cout << (qtilde_ ? "Restricting" : "Not restricting") << " " << r->GetName() << " to positive values." << std::endl;
   if (verbose > 1) params_->Print("V");
  
   if (verbose > 0) std::cout << "\nMake global fit of real data" << std::endl;
   {
     CloseCoutSentry sentry(verbose < 3);
     *params_ = snapGlobalObsData;
     CascadeMinimizer minim(*nllD_, CascadeMinimizer::Unconstrained, r);
     minim.setStrategy(minimizerStrategy_);
     minim.minimize(verbose-2);
     fitFreeD_.readFrom(*params_);
     minNllD_ = nllD_->getVal();
   }
   if (verbose > 0) std::cout << "NLL at global minimum of data: " << minNllD_ << " (" << r->GetName() << " = " << r->getVal() << ")" << std::endl;
   double rErr = std::max<double>(r->getError(), 0.02 * (r->getMax() - r->getMin()));
 
   r->setMin(0);
 
   if (verbose > 1) fitFreeD_.Print("V");
   if (verbose > 0) std::cout << "\nMake global fit of asimov data" << std::endl;
   {
     CloseCoutSentry sentry(verbose < 3);
     *params_ = snapGlobalObsAsimov;
     CascadeMinimizer minim(*nllA_, CascadeMinimizer::Unconstrained, r);
     minim.setStrategy(minimizerStrategy_);
     minim.minimize(verbose-2);
     fitFreeA_.readFrom(*params_);
     minNllA_ = nllA_->getVal();
   }
   if (verbose > 0) std::cout << "NLL at global minimum of asimov: " << minNllA_ << " (" << r->GetName() << " = " << r->getVal() << ")" << std::endl;
   if (verbose > 1) fitFreeA_.Print("V");
 
   fitFreeD_.writeTo(*params_);
   r->setConstant(true);
 
   if (what_ == "singlePoint") {
     limit = getCLs(*r, rValue_, true, &limit, &limitErr);
     return true;
   }
 
   double clsTarget = 1-cl;
   double rMin = std::max<double>(0, r->getVal()), rMax = rMin + 3 * rErr;
   double clsMax = 1, clsMin = 0;
   for (int tries = 0; tries < 5; ++tries) {
     double cls = getCLs(*r, rMax);
     if (cls == -999) { std::cerr << "Minimization failed in an unrecoverable way" << std::endl; break; }
     if (cls < clsTarget) { clsMin = cls; break; }
     rMax *= 2;
   }
   
   do {
     if (clsMax < 3*clsTarget && clsMin > 0.3*clsTarget) {
         double rCross = rMin + (rMax-rMin)*log(clsMax/clsTarget)/log(clsMax/clsMin);
         if ((rCross-rMin) < (rMax - rCross)) {
             limit = 0.8*rCross + 0.2*rMax;
         } else {
             limit = 0.8*rCross + 0.2*rMin;
         }
         limitErr = 0.5*(rMax - rMin);
     } else {
         limit = 0.5*(rMin + rMax); 
         limitErr = 0.5*(rMax - rMin);
     }
     double cls = getCLs(*r, limit);
     if (cls == -999) { std::cerr << "Minimization failed in an unrecoverable way" << std::endl; break; }
     if (cls > clsTarget) {
         clsMax = cls;
         rMin = limit;
     } else {
         clsMin = cls;
         rMax = limit;
     }
   } while (limitErr > std::max(rRelAccuracy_ * limit, rAbsAccuracy_));
 
   return true;
 }

std::vector< std::pair< float, float > > Asymptotic::runLimitExpected	(	RooWorkspace *	w,
		RooStats::ModelConfig *	mc_s,
		RooStats::ModelConfig *	mc_b,
		RooAbsData &	data,
		double &	limit,
		double &	limitErr,
		const double *	hint
	)

Definition at line 263 of file Asymptotic.cc.

References alpha, asimovDataset(), cl, CloseCoutSentry::clear(), Combine::commitPoint(), gather_cfg::cout, alignCSCRings::e, findExpectedLimitFromCrossing(), edm::detail::isnan(), minimizerStrategy_, minosAlgo_, N, newExpected_, params_, funct::pow(), qtilde_, alignCSCRings::r, snapGlobalObsAsimov, and CascadeMinimizer::Unconstrained.

Referenced by run().

                                                                                                                                                                                                           {
     // See equation 35-38 of AN 2011/298 and references cited therein
     //
     //   (35)  sigma^2   = mu^2 / q_mu(Asimov)
     //   (38)  mu_median = sigma * normal_quantile(1-0.5*(1-cl))
     //
     // -->  q_mu(Asimov) = pow(normal_quantile(1-0.5*(1-cl)), 2)
     //      can be solved to find mu_median
     //
     // --> then (38) gives sigma, and the quantiles are given by (37)
     //      mu_N = sigma * (normal_quantile(1 - quantile*(1-cl), 1.0) + normal_quantile(quantile));
     //
     // 1) get parameter of interest
     RooArgSet  poi(*mc_s->GetParametersOfInterest());
     RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
 
     // 2) get asimov dataset
     RooAbsData *asimov = asimovDataset(w, mc_s, mc_b, data);
 
     // 2b) load asimov global observables
     if (params_.get() == 0) params_.reset(mc_s->GetPdf()->getParameters(data));
     *params_ = snapGlobalObsAsimov;
 
     // 3) solve for q_mu
     r->setConstant(false);
     //r->setMin(0);
     r->setMin(qtilde_ ? 0 : -r->getMax()); // FIXME TEST
     r->setVal(0.01*r->getMax());
     r->setError(0.1*r->getMax());
     //r->removeMax();
     
     std::auto_ptr<RooAbsReal> nll(mc_s->GetPdf()->createNLL(*asimov, RooFit::Constrain(*mc_s->GetNuisanceParameters())));
     CascadeMinimizer minim(*nll, CascadeMinimizer::Unconstrained, r);
     minim.setStrategy(minimizerStrategy_);
     minim.setErrorLevel(0.5*pow(ROOT::Math::normal_quantile(1-0.5*(1-cl),1.0), 2)); // the 0.5 is because qmu is -2*NLL
                         // eventually if cl = 0.95 this is the usual 1.92!
     CloseCoutSentry sentry(verbose < 3);    
     minim.minimize(verbose-2);
     sentry.clear();
     if (verbose > 1) {
         std::cout << "Fit to asimov dataset:" << std::endl;
         std::auto_ptr<RooFitResult> res(minim.save());
         res->Print("V");
     }
     if (r->getVal()/r->getMax() > 1e-3) {
         if (verbose) printf("WARNING: Best fit of asimov dataset is at %s = %f (%f times %sMax), while it should be at zero\n",
                 r->GetName(), r->getVal(), r->getVal()/r->getMax(), r->GetName());
     }
 
 
     // 3) get ingredients for equation 37
     double nll0 = nll->getVal();
     double median = findExpectedLimitFromCrossing(*nll, r, r->getMin(), r->getMax(), nll0, 0.5);
     double sigma  = median / ROOT::Math::normal_quantile(1-0.5*(1-cl),1.0);
     double alpha = 1-cl;
     if (verbose > 0) { 
         std::cout << "Median for expected limits: " << median << std::endl; 
         std::cout << "Sigma  for expected limits: " << sigma  << std::endl; 
     }
 
     std::vector<std::pair<float,float> > expected;
     const double quantiles[5] = { 0.025, 0.16, 0.50, 0.84, 0.975 };
     for (int iq = 0; iq < 5; ++iq) {
         double N = ROOT::Math::normal_quantile(quantiles[iq], 1.0);
         if (newExpected_ && iq != 2) { // the median is exactly the same in the two methods
             std::string minosAlgoBackup = minosAlgo_;
             if (minosAlgo_ == "stepping") minosAlgo_ = "bisection";
             switch (iq) {
                 case 0: limit = findExpectedLimitFromCrossing(*nll, r, r->getMin(),            median,      nll0, quantiles[iq]); break;
                 case 1: limit = findExpectedLimitFromCrossing(*nll, r, expected.back().second, median,      nll0, quantiles[iq]); break;
                 case 3: limit = findExpectedLimitFromCrossing(*nll, r, expected.back().second, median+2*sigma, nll0, quantiles[iq]); break;
                 case 4: limit = findExpectedLimitFromCrossing(*nll, r, expected.back().second, median+4*sigma, nll0, quantiles[iq]); break;
             }
             minosAlgo_ = minosAlgoBackup;
             if (std::isnan(limit)) { expected.clear(); break; } 
         } else {
             limit = sigma*(ROOT::Math::normal_quantile(1 - alpha * quantiles[iq], 1.0) + N);
         }
         limitErr = 0;
         Combine::commitPoint(true, quantiles[iq]);
         expected.push_back(std::pair<float,float>(quantiles[iq], limit));
     }
     return expected;
 
 }