#include <FitterAlgoBase.h>

Inheritance diagram for FitterAlgoBase:

Public Member Functions
void	applyOptionsBase (const boost::program_options::variables_map &vm)

	FitterAlgoBase (const char *title="<FillMe> specific options")

virtual bool	run (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
virtual void	applyDefaultOptions ()

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

	LimitAlgo ()

	LimitAlgo (const char *desc)

virtual const std::string &	name () const =0

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

virtual void	setNToys (const int)

virtual void	setToyNumber (const int)

virtual	~LimitAlgo ()

Protected Member Functions
RooFitResult *	doFit (RooAbsPdf &pdf, RooAbsData &data, RooRealVar &r, const RooCmdArg &constrain, bool doHesse=true, int ndim=1, bool reuseNLL=false)

RooFitResult *	doFit (RooAbsPdf &pdf, RooAbsData &data, const RooArgList &rs, const RooCmdArg &constrain, bool doHesse=true, int ndim=1, bool reuseNLL=false)

double	findCrossing (CascadeMinimizer &minim, RooAbsReal &nll, RooRealVar &r, double level, double rStart, double rBound)

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

Protected Attributes
std::auto_ptr< RooAbsReal >	nll

Protected Attributes inherited from LimitAlgo
boost::program_options::options_description	options_

Static Protected Attributes
static bool	do95_ = false

static bool	keepFailures_ = false

static int	maxFailedSteps_ = 5

static std::string	minimizerAlgo_ = "Minuit2"

static std::string	minimizerAlgoForMinos_ = "Minuit2,simplex"

static int	minimizerStrategy_ = 1

static int	minimizerStrategyForMinos_ = 0

static float	minimizerTolerance_ = 1e-2

static float	minimizerToleranceForMinos_ = 1e-4

static float	nllValue_ = std::numeric_limits<float>::quiet_NaN()

static float	preFitValue_ = 1.0

static bool	robustFit_ = false

static bool	saveNLL_ = false

static float	stepSize_ = 0.1

Detailed Description

Do a ML fit of the data with background and signal+background hypothesis and print out diagnostics plots

Author: Giovanni Petrucciani (UCSD)

Definition at line 20 of file FitterAlgoBase.h.

Constructor & Destructor Documentation

FitterAlgoBase::FitterAlgoBase ( const char * title = "<FillMe> specific options" )

Definition at line 53 of file FitterAlgoBase.cc.

References do95_, maxFailedSteps_, minimizerAlgo_, minimizerAlgoForMinos_, minimizerStrategy_, minimizerStrategyForMinos_, minimizerTolerance_, minimizerToleranceForMinos_, LimitAlgo::options_, preFitValue_, robustFit_, and stepSize_.

                                                 :
     LimitAlgo(title)
 {
     options_.add_options()
         ("minimizerAlgo",      boost::program_options::value<std::string>(&minimizerAlgo_)->default_value(minimizerAlgo_), "Choice of minimizer (Minuit vs Minuit2)")
         ("minimizerTolerance", boost::program_options::value<float>(&minimizerTolerance_)->default_value(minimizerTolerance_),  "Tolerance for minimizer")
         ("minimizerStrategy",  boost::program_options::value<int>(&minimizerStrategy_)->default_value(minimizerStrategy_),      "Stragegy for minimizer")
         ("preFitValue",        boost::program_options::value<float>(&preFitValue_)->default_value(preFitValue_),  "Value of signal strength pre-fit")
         ("do95",       boost::program_options::value<bool>(&do95_)->default_value(do95_),  "Compute also 2-sigma interval from delta(nll) = 1.92 instead of 0.5")
         ("robustFit",  boost::program_options::value<bool>(&robustFit_)->default_value(robustFit_),  "Search manually for 1 and 2 sigma bands instead of using Minos")
         ("maxFailedSteps",  boost::program_options::value<int>(&maxFailedSteps_)->default_value(maxFailedSteps_),  "How many failed steps to retry before giving up")
         ("stepSize",        boost::program_options::value<float>(&stepSize_)->default_value(stepSize_),  "Step size for robust fits (multiplier of the range)")
         ("minimizerAlgoForMinos",      boost::program_options::value<std::string>(&minimizerAlgoForMinos_)->default_value(minimizerAlgoForMinos_), "Choice of minimizer (Minuit vs Minuit2) for profiling in robust fits")
         ("minimizerStrategyForMinos",  boost::program_options::value<int>(&minimizerStrategyForMinos_)->default_value(minimizerStrategyForMinos_),      "Stragegy for minimizer for profiling in robust fits")
         ("minimizerToleranceForMinos",  boost::program_options::value<float>(&minimizerToleranceForMinos_)->default_value(minimizerToleranceForMinos_),      "Tolerance for minimizer for profiling in robust fits")
         ("saveNLL",  "Save the negative log-likelihood at the minimum in the output tree (note: value is relative to the pre-fit state)")
         ("keepFailures",  "Save the results even if the fit is declared as failed (for NLL studies)")
     ;
 }

Member Function Documentation

void FitterAlgoBase::applyOptionsBase ( const boost::program_options::variables_map & vm )

Definition at line 73 of file FitterAlgoBase.cc.

References keepFailures_, and saveNLL_.

Referenced by ChannelCompatibilityCheck::applyOptions(), MaxLikelihoodFit::applyOptions(), and MultiDimFit::applyOptions().

 {
     saveNLL_ = vm.count("saveNLL");
     keepFailures_ = vm.count("keepFailures");
 }

RooFitResult * FitterAlgoBase::doFit	(	RooAbsPdf &	pdf,
		RooAbsData &	data,
		RooRealVar &	r,
		const RooCmdArg &	constrain,
		bool	doHesse = `true`,
		int	ndim = `1`,
		bool	reuseNLL = `false`
	)

protected

Fit data with pdf, with parameters of interest in r, and specified constraint If ndim = 1, errors on each parameter are from a 1-dim chisquare, as for a single parameter fit If ndim > 1, errors on each parameter are from a n-dim chisquare, as for a joint estimation of N parameters

Definition at line 90 of file FitterAlgoBase.cc.

Referenced by MultiDimFit::runSpecific(), MaxLikelihoodFit::runSpecific(), and ChannelCompatibilityCheck::runSpecific().

                                                                                                                                                       {
     return doFit(pdf, data, RooArgList (r), constrain, doHesse, ndim, reuseNLL);
 }

RooFitResult * FitterAlgoBase::doFit	(	RooAbsPdf &	pdf,
		RooAbsData &	data,
		const RooArgList &	rs,
		const RooCmdArg &	constrain,
		bool	doHesse = `true`,
		int	ndim = `1`,
		bool	reuseNLL = `false`
	)

protected

Definition at line 94 of file FitterAlgoBase.cc.

References CloseCoutSentry::clear(), CascadeMinimizer::Constrained, gather_cfg::cout, do95_, findCrossing(), i, CascadeMinimizer::improve(), edm::detail::isnan(), keepFailures_, max(), CascadeMinimizer::minimize(), CascadeMinimizer::minimizer(), minimizerStrategy_, minimizerStrategyForMinos_, n, nll, nllValue_, convertSQLiteXML::ok, alignCSCRings::r, run_regression::ret, robustFit_, CascadeMinimizer::save(), CascadeMinimizer::setErrorLevel(), CascadeMinimizer::setStrategy(), and CascadeMinimizer::Unconstrained.

                                                                                                                                                              {
     RooFitResult *ret = 0;
     if (reuseNLL) nll->setData(data);   // reuse nll but swap out the data
     else nll.reset(pdf.createNLL(data, constrain, RooFit::Extended(pdf.canBeExtended()))); // make a new nll
 
     double nll0 = nll->getVal();
     double delta68 = 0.5*ROOT::Math::chisquared_quantile_c(1-0.68,ndim);
     double delta95 = 0.5*ROOT::Math::chisquared_quantile_c(1-0.95,ndim);
     CascadeMinimizer minim(*nll, CascadeMinimizer::Unconstrained, rs.getSize() ? dynamic_cast<RooRealVar*>(rs.first()) : 0);
     minim.setStrategy(minimizerStrategy_);
     minim.setErrorLevel(delta68);
     CloseCoutSentry sentry(verbose < 3);    
     bool ok = minim.minimize(verbose-1);
     nllValue_ =  nll->getVal() - nll0;
     if (!ok && !keepFailures_) { std::cout << "Initial minimization failed. Aborting." << std::endl; return 0; }
     if (doHesse) minim.minimizer().hesse();
     sentry.clear();
     ret = minim.save();
     if (verbose > 1) { ret->Print("V");  }
 
     std::auto_ptr<RooArgSet> allpars(pdf.getParameters(data));
 
     for (int i = 0, n = rs.getSize(); i < n; ++i) {
         // if this is not the first fit, reset parameters  
         if (i) {
             RooArgSet oldparams(ret->floatParsFinal());
             *allpars = oldparams;
         }
     
         // get the parameter to scan, amd output variable in fit result
         RooRealVar &r = dynamic_cast<RooRealVar &>(*rs.at(i));
         RooRealVar &rf = dynamic_cast<RooRealVar &>(*ret->floatParsFinal().find(r.GetName()));
         double r0 = r.getVal(), rMin = r.getMin(), rMax = r.getMax();
 
         if (!robustFit_) {
             if (do95_) {
                 throw std::runtime_error("95% CL errors with Minos are not working at the moment.");
                 minim.setErrorLevel(delta95);
                 minim.improve(verbose-1);
                 minim.setErrorLevel(delta95);
                 if (minim.minimizer().minos(RooArgSet(r)) != -1) {
                     rf.setRange("err95", r.getVal() + r.getAsymErrorLo(), r.getVal() + r.getAsymErrorHi());
                 }
                 minim.setErrorLevel(delta68);
                 minim.improve(verbose-1);
             }
             if (minim.minimizer().minos(RooArgSet(r)) != -1) {
                 rf.setRange("err68", r.getVal() + r.getAsymErrorLo(), r.getVal() + r.getAsymErrorHi());
                 rf.setAsymError(r.getAsymErrorLo(), r.getAsymErrorHi());
             }
        } else {
             r.setVal(r0); r.setConstant(true);
             
             CascadeMinimizer minim2(*nll, CascadeMinimizer::Constrained);
             minim2.setStrategy(minimizerStrategyForMinos_);
 
             std::auto_ptr<RooArgSet> allpars(nll->getParameters((const RooArgSet *)0));
 
             double nll0 = nll->getVal();
             double threshold68 = nll0 + delta68;
             double threshold95 = nll0 + delta95;
             // search for crossings
 
             assert(!std::isnan(r0));
             // high error
             double hi68 = findCrossing(minim2, *nll, r, threshold68, r0,   rMax);
             double hi95 = do95_ ? findCrossing(minim2, *nll, r, threshold95, std::isnan(hi68) ? r0 : hi68, std::max(rMax, std::isnan(hi68*2-r0) ? r0 : hi68*2-r0)) : r0;
             // low error 
             *allpars = RooArgSet(ret->floatParsFinal()); r.setVal(r0); r.setConstant(true);
             double lo68 = findCrossing(minim2, *nll, r, threshold68, r0,   rMin); 
             double lo95 = do95_ ? findCrossing(minim2, *nll, r, threshold95, std::isnan(lo68) ? r0 : lo68, rMin) : r0;
 
             rf.setAsymError(!std::isnan(lo68) ? lo68 - r0 : 0, !std::isnan(hi68) ? hi68 - r0 : 0);
             rf.setRange("err68", !std::isnan(lo68) ? lo68 : r0, !std::isnan(hi68) ? hi68 : r0);
             if (do95_ && (!std::isnan(lo95) || !std::isnan(hi95))) {
                 rf.setRange("err95", !std::isnan(lo95) ? lo95 : r0, !std::isnan(hi95) ? hi95 : r0);
             }
 
             r.setVal(r0); r.setConstant(false);
         }
     }
 
     return ret;
 }

double FitterAlgoBase::findCrossing	(	CascadeMinimizer &	minim,
		RooAbsReal &	nll,
		RooRealVar &	r,
		double	level,
		double	rStart,
		double	rBound
	)

protected

Definition at line 179 of file FitterAlgoBase.cc.

References gather_cfg::cout, runtimedef::get(), CascadeMinimizer::improve(), max(), maxFailedSteps_, min, minimizerAlgoForMinos_, minimizerToleranceForMinos_, convertSQLiteXML::ok, CascadeMinimizer::save(), and stepSize_.

Referenced by MultiDimFit::doBox(), MultiDimFit::doContour2D(), and doFit().

                                                                                                                                        {
     ProfileLikelihood::MinimizerSentry minimizerConfig(minimizerAlgoForMinos_, minimizerToleranceForMinos_);
     if (verbose) std::cout << "Searching for crossing at nll = " << level << " in the interval " << rStart << ", " << rBound << std::endl; 
     double rInc = stepSize_*(rBound - rStart);
     r.setVal(rStart); 
     std::auto_ptr<RooFitResult> checkpoint;
     std::auto_ptr<RooArgSet>    allpars;
     bool ok = false;
     {
         CloseCoutSentry sentry(verbose < 3);    
         ok = minim.improve(verbose-1);
         checkpoint.reset(minim.save());
     }
     if (!ok) { std::cout << "Error: minimization failed at " << r.GetName() << " = " << rStart << std::endl; return NAN; }
     double here = nll.getVal();
     int nfail = 0;
     if (verbose > 0) { printf("      %s      lvl-here  lvl-there   stepping\n", r.GetName()); fflush(stdout); }
     do {
         rStart += rInc;
         if (rInc*(rStart - rBound) > 0) { // went beyond bounds
             rStart -= rInc;
             rInc    = 0.5*(rBound-rStart);
         }
         r.setVal(rStart);
         nll.clearEvalErrorLog(); nll.getVal();
         if (nll.numEvalErrors() > 0) {
             ok = false;
         } else {
             CloseCoutSentry sentry(verbose < 3);    
             ok = minim.improve(verbose-1);
         }
         if (!ok) { 
             nfail++;
             if (nfail >= maxFailedSteps_) {  std::cout << "Error: minimization failed at " << r.GetName() << " = " << rStart << std::endl; return NAN; }
             RooArgSet oldparams(checkpoint->floatParsFinal());
             if (allpars.get() == 0) allpars.reset(nll.getParameters((const RooArgSet *)0));
             *allpars = oldparams;
             rStart -= rInc; rInc *= 0.5; 
             continue;
         } else nfail = 0;
         double there = here;
         here = nll.getVal();
         if (verbose > 0) { printf("%f    %+.5f  %+.5f    %f\n", rStart, level-here, level-there, rInc); fflush(stdout); }
         if ( fabs(here - level) < 4*minimizerToleranceForMinos_ ) {
             // set to the right point with interpolation
             r.setVal(rStart + (level-here)*(level-there)/(here-there));
             return r.getVal();
         } else if (here > level) {
             // I'm above the level that I wanted, this means I stepped too long
             // First I take back all the step
             rStart -= rInc; 
             // Then I try to figure out a better step
             if (runtimedef::get("FITTER_DYN_STEP")) {
                 if (fabs(there - level) > 0.05) { // If started from far away, I still have to step carefully
                     double rIncFactor = std::max(0.2, std::min(0.7, 0.75*(level-there)/(here-there)));
                     //printf("\t\t\t\t\tCase A1: level-there = %f,  here-there = %f,   rInc(Old) = %f,  rInFactor = %f,  rInc(New) = %f\n", level-there, here-there, rInc, rIncFactor, rInc*rIncFactor);
                     rInc *= rIncFactor;
                 } else { // close enough to jump straight to target
                     double rIncFactor = std::max(0.05, std::min(0.95, 0.95*(level-there)/(here-there)));
                     //printf("\t\t\t\t\tCase A2: level-there = %f,  here-there = %f,   rInc(Old) = %f,  rInFactor = %f,  rInc(New) = %f\n", level-there, here-there, rInc, rIncFactor, rInc*rIncFactor);
                     rInc *= rIncFactor;
                 }
             } else {
                 rInc *= 0.3;
             }
             if (allpars.get() == 0) allpars.reset(nll.getParameters((const RooArgSet *)0));
             RooArgSet oldparams(checkpoint->floatParsFinal());
             *allpars = oldparams;
         } else if ((here-there)*(level-there) < 0 && // went wrong
                    fabs(here-there) > 0.1) {         // by more than roundoff
             if (allpars.get() == 0) allpars.reset(nll.getParameters((const RooArgSet *)0));
             RooArgSet oldparams(checkpoint->floatParsFinal());
             *allpars = oldparams;
             rStart -= rInc; rInc *= 0.5;
         } else {
             // I did this step, and I'm not there yet
             if (runtimedef::get("FITTER_DYN_STEP")) {
                 if (fabs(here - level) > 0.05) { // we still have to step carefully
                     if ((here-there)*(level-there) > 0) { // if we went in the right direction
                         // then optimize step size
                         double rIncFactor = std::max(0.2, std::min(2.0, 0.75*(level-there)/(here-there)));
                         //printf("\t\t\t\t\tCase B1: level-there = %f,  here-there = %f,   rInc(Old) = %f,  rInFactor = %f,  rInc(New) = %f\n", level-there, here-there, rInc, rIncFactor, rInc*rIncFactor);
                         rInc *= rIncFactor;
                     } //else printf("\t\t\t\t\tCase B3: level-there = %f,  here-there = %f,   rInc(Old) = %f\n", level-there, here-there, rInc);
                 } else { // close enough to jump straight to target
                     double rIncFactor = std::max(0.05, std::min(4.0, 0.95*(level-there)/(here-there)));
                     //printf("\t\t\t\t\tCase B2: level-there = %f,  here-there = %f,   rInc(Old) = %f,  rInFactor = %f,  rInc(New) = %f\n", level-there, here-there, rInc, rIncFactor, rInc*rIncFactor);
                     rInc *= rIncFactor;
                 }
             } else {
                 //nothing?
             }
             checkpoint.reset(minim.save());
         }
     } while (fabs(rInc) > minimizerToleranceForMinos_*stepSize_*std::max(1.0,rBound-rStart));
     if (fabs(here - level) > 0.01) {
         std::cout << "Error: closed range without finding crossing." << std::endl;
         return NAN;
     } else {
         return r.getVal();
     }
 }

bool FitterAlgoBase::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 79 of file FitterAlgoBase.cc.

References Combine::addBranch(), minimizerAlgo_, minimizerTolerance_, nllValue_, runSpecific(), and saveNLL_.

                                                                                                                                                                      { 
   ProfileLikelihood::MinimizerSentry minimizerConfig(minimizerAlgo_, minimizerTolerance_);
   CloseCoutSentry sentry(verbose < 0);
 
   static bool shouldCreateNLLBranch = saveNLL_;
   if (shouldCreateNLLBranch) { Combine::addBranch("nll", &nllValue_, "nll/F"); shouldCreateNLLBranch = false; }
 
   return runSpecific(w, mc_s, mc_b, data, limit, limitErr, hint);
 }

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