diff --git a/[refs] b/[refs] index 98cc434895fbf..67a6426e5fb6e 100644 --- a/[refs] +++ b/[refs] @@ -1,2 +1,2 @@ --- -refs/heads/CMSSW_7_0_X: 8d707f24cff1d0e3fdb5bca850b657d54fc2111a +refs/heads/CMSSW_7_0_X: d2a66e3acfb93dac540352c1a9afb14c7c37895c diff --git a/trunk/HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h b/trunk/HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h index 2a2f4249854eb..2645cb036ede3 100644 --- a/trunk/HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h +++ b/trunk/HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h @@ -33,6 +33,8 @@ class CascadeMinimizer { int strategy_; RooRealVar * poi_; + bool improveOnce(int verbose); + /// options configured from command line static boost::program_options::options_description options_; /// compact information about an algorithm @@ -51,6 +53,8 @@ class CascadeMinimizer { static bool singleNuisFit_; /// do first a fit of only the POI static bool setZeroPoint_; + /// don't do old fallback using robustMinimize + static bool oldFallback_; }; #endif diff --git a/trunk/HiggsAnalysis/CombinedLimit/interface/SequentialMinimizer.h b/trunk/HiggsAnalysis/CombinedLimit/interface/SequentialMinimizer.h index b3b61090684b1..10d9303d60dc3 100644 --- a/trunk/HiggsAnalysis/CombinedLimit/interface/SequentialMinimizer.h +++ b/trunk/HiggsAnalysis/CombinedLimit/interface/SequentialMinimizer.h @@ -63,7 +63,7 @@ class SequentialMinimizer { bool minimize(double ytol=0, int bigsteps=0, int smallsteps=5); bool improve(double ytol=0, int bigsteps=0, int smallsteps=5); private: - enum State { Cleared, Ready, Active, Done }; + enum State { Cleared, Ready, Active, Done, Fixed }; struct Worker : public OneDimMinimizer { Worker() : OneDimMinimizer(), xtol(0), state(Cleared) {} Worker(RooAbsReal *nll, RooRealVar *var) : OneDimMinimizer(nll,var), xtol(0), state(Cleared) {} @@ -72,5 +72,6 @@ class SequentialMinimizer { }; RooAbsReal *nll_; std::vector workers_; + State state_; }; #endif diff --git a/trunk/HiggsAnalysis/CombinedLimit/src/CascadeMinimizer.cc b/trunk/HiggsAnalysis/CombinedLimit/src/CascadeMinimizer.cc index 5ef34f366dea4..eabf2171c11ef 100644 --- a/trunk/HiggsAnalysis/CombinedLimit/src/CascadeMinimizer.cc +++ b/trunk/HiggsAnalysis/CombinedLimit/src/CascadeMinimizer.cc @@ -13,6 +13,7 @@ bool CascadeMinimizer::preScan_; bool CascadeMinimizer::poiOnlyFit_; bool CascadeMinimizer::singleNuisFit_; bool CascadeMinimizer::setZeroPoint_; +bool CascadeMinimizer::oldFallback_ = true; CascadeMinimizer::CascadeMinimizer(RooAbsReal &nll, Mode mode, RooRealVar *poi, int initialStrategy) : nll_(nll), @@ -27,25 +28,17 @@ bool CascadeMinimizer::improve(int verbose, bool cascade) { if (setZeroPoint_) { cacheutils::CachingSimNLL *simnll = dynamic_cast(&nll_); - if (simnll) simnll->setZeroPoint(); + if (simnll) { + simnll->setZeroPoint(); + } } minimizer_.setPrintLevel(verbose-2); minimizer_.setStrategy(strategy_); - std::string nominalType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); - std::string nominalAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); - float nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); - bool outcome = false; - if (nominalType == "Sequential") { - if (seqmin_.get() == 0) { - seqmin_.reset(new SequentialMinimizer(& nll_, mode_ == Unconstrained ? poi_ : 0)); - seqmin_->minimize(nominalTol); - } else { - seqmin_->improve(nominalTol); - } - } else { - outcome = nllutils::robustMinimize(nll_, minimizer_, verbose); - } + bool outcome = improveOnce(verbose-2); if (cascade && !outcome && !fallbacks_.empty()) { + std::string nominalType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); + std::string nominalAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); + float nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); if (verbose > 0) std::cerr << "Failed minimization with " << nominalType << "," << nominalAlgo << " and tolerance " << nominalTol << std::endl; for (std::vector::const_iterator it = fallbacks_.begin(), ed = fallbacks_.end(); it != ed; ++it) { ProfileLikelihood::MinimizerSentry minimizerConfig(it->algo, it->tolerance != -1.f ? it->tolerance : nominalTol); @@ -53,7 +46,7 @@ bool CascadeMinimizer::improve(int verbose, bool cascade) nominalAlgo != ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo() || nominalTol != ROOT::Math::MinimizerOptions::DefaultTolerance()) { if (verbose > 0) std::cerr << "Will fallback to minimization using " << it->algo << " and tolerance " << it->tolerance << std::endl; - outcome = nllutils::robustMinimize(nll_, minimizer_, verbose); + outcome = improveOnce(verbose-2); if (outcome) break; } } @@ -65,6 +58,28 @@ bool CascadeMinimizer::improve(int verbose, bool cascade) return outcome; } +bool CascadeMinimizer::improveOnce(int verbose) +{ + std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); + std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); + float myTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); + bool outcome = false; + if (myType == "Sequential") { + if (seqmin_.get() == 0) { + seqmin_.reset(new SequentialMinimizer(& nll_, mode_ == Unconstrained ? poi_ : 0)); + outcome = seqmin_->minimize(myTol); + } else { + outcome = seqmin_->improve(myTol); + } + } else if (oldFallback_){ + outcome = nllutils::robustMinimize(nll_, minimizer_, verbose); + } else { + int status = minimizer_.minimize(myType.c_str(), myAlgo.c_str()); + outcome = (status == 0); + } + return outcome; +} + bool CascadeMinimizer::minimize(int verbose, bool cascade) { minimizer_.setPrintLevel(verbose-2); @@ -86,6 +101,7 @@ void CascadeMinimizer::initOptions() ("cminSingleNuisFit", "Do first a minimization of each nuisance parameter individually") ("cminFallbackAlgo", boost::program_options::value >(), "Fallback algorithms if the default minimizer fails (can use multiple ones). Syntax is algo[,subalgo][:tolerance]") ("cminSetZeroPoint", "Change the reference point of the NLL to be zero during minimization") + ("cminOldRobustMinimize", boost::program_options::value(&oldFallback_)->default_value(oldFallback_), "Use the old 'robustMinimize' logic in addition to the cascade") ; } diff --git a/trunk/HiggsAnalysis/CombinedLimit/src/ProfileLikelihood.cc b/trunk/HiggsAnalysis/CombinedLimit/src/ProfileLikelihood.cc index 4d9da87d3725c..d11d39e43d772 100644 --- a/trunk/HiggsAnalysis/CombinedLimit/src/ProfileLikelihood.cc +++ b/trunk/HiggsAnalysis/CombinedLimit/src/ProfileLikelihood.cc @@ -386,7 +386,7 @@ double ProfileLikelihood::significanceBruteForce(RooAbsPdf &pdf, RooAbsData &dat } poi.setConstant(true); CascadeMinimizer minim(*nll, CascadeMinimizer::Constrained); - if (!minim.improve(verbose-2)) { + if (!minim.minimize(verbose-2)) { std::cerr << "Initial minimization failed. Aborting." << std::endl; return -1; } else if (verbose > 0) { @@ -396,10 +396,14 @@ double ProfileLikelihood::significanceBruteForce(RooAbsPdf &pdf, RooAbsData &dat std::auto_ptr start(minim.save()); double minnll = nll->getVal(), thisnll = minnll, lastnll = thisnll; double rbest = poi.getVal(), rval = rbest; + TGraph *points = 0; if (verbose) { printf(" %-6s delta(NLL)\n", poi.GetName()); printf("%8.5f %8.5f\n", rval, 0.); fflush(stdout); + points = new TGraph(1); + points->SetName(Form("nll_scan_%g", mass_)); + points->SetPoint(0, rval, 0); } while (rval >= tolerance * poi.getMax()) { rval *= 0.8; @@ -412,7 +416,11 @@ double ProfileLikelihood::significanceBruteForce(RooAbsPdf &pdf, RooAbsData &dat std::cerr << "Minimization failed at " << poi.getVal() <<". exiting the loop" << std::endl; return -1; } - if (verbose) { printf("%8.5f %8.5f\n", rval, thisnll-minnll); fflush(stdout); } + if (verbose) { + printf("%8.5f %8.5f\n", rval, thisnll-minnll); fflush(stdout); + points->Set(points->GetN()+1); + points->SetPoint(points->GetN()-1, rval, thisnll - minnll); + } if (fabs(lastnll - thisnll) < 7*minimizerToleranceForBF_) { std::cout << "This is enough." << std::endl; if (thisnll < lastnll) { @@ -438,6 +446,7 @@ double ProfileLikelihood::significanceBruteForce(RooAbsPdf &pdf, RooAbsData &dat } #endif } + if (points) outputFile->WriteTObject(points); return (thisnll - minnll); } diff --git a/trunk/HiggsAnalysis/CombinedLimit/src/ProfiledLikelihoodRatioTestStatExt.cc b/trunk/HiggsAnalysis/CombinedLimit/src/ProfiledLikelihoodRatioTestStatExt.cc index 0fda7a8c80bcb..342cd05dbc08e 100644 --- a/trunk/HiggsAnalysis/CombinedLimit/src/ProfiledLikelihoodRatioTestStatExt.cc +++ b/trunk/HiggsAnalysis/CombinedLimit/src/ProfiledLikelihoodRatioTestStatExt.cc @@ -14,7 +14,7 @@ #include "../interface/ProfilingTools.h" //---- Uncomment this and run with --perfCounters to get statistics of successful and failed fits -//#define DEBUG_FIT_STATUS +#define DEBUG_FIT_STATUS #ifdef DEBUG_FIT_STATUS #define COUNT_ONE(x) PerfCounter::add(x); #else @@ -147,6 +147,7 @@ ProfiledLikelihoodTestStatOpt::ProfiledLikelihoodTestStatOpt( Double_t ProfiledLikelihoodTestStatOpt::Evaluate(RooAbsData& data, RooArgSet& /*nullPOI*/) { bool do_debug = runtimedef::get("DEBUG_PLTSO"); + //static bool do_rescan = runtimedef::get("PLTSO_FULL_RESCAN"); DBG(DBG_PLTestStat_main, (std::cout << "Being evaluated on " << data.GetName() << std::endl)) // Take snapshot of initial state, to restore it at the end @@ -207,6 +208,19 @@ Double_t ProfiledLikelihoodTestStatOpt::Evaluate(RooAbsData& data, RooArgSet& /* thisNLL = nullNLL; } } + /* This piece of debug code was added to see if we were finding a local minimum at zero instead of the global minimum. + It doesn't seem to be the case, however + if (do_rescan && fabs(thisNLL - nullNLL) < 0.2 && initialR == 0) { + if (do_debug) printf("Doing full rescan. best fit r = %.4f, -lnQ = %.5f\n", bestFitR, thisNLL-nullNLL); + for (double rx = 0; rx <= 5; rx += 0.01) { + r->setVal(rx); double y = nll_->getVal(); + if (y < nullNLL) { + printf("%4.2f\t%8.5f\t<<==== ALERT\n", rx, y - nullNLL); + } else { + printf("%4.2f\t%8.5f\n", rx, y - nullNLL); + } + } + } */ if (bestFitR > initialR && oneSided_ == signFlipDef) { DBG(DBG_PLTestStat_main, (printf(" fitted signal %7.4f > %7.4f, test statistics will be negative.\n", bestFitR, initialR))) std::swap(thisNLL, nullNLL); diff --git a/trunk/HiggsAnalysis/CombinedLimit/src/SequentialMinimizer.cc b/trunk/HiggsAnalysis/CombinedLimit/src/SequentialMinimizer.cc index 4ba9e73bf315a..9b2ddb50aac95 100644 --- a/trunk/HiggsAnalysis/CombinedLimit/src/SequentialMinimizer.cc +++ b/trunk/HiggsAnalysis/CombinedLimit/src/SequentialMinimizer.cc @@ -9,13 +9,22 @@ #include "RooLinkedListIter.h" #include #include +#include "../interface/ProfilingTools.h" #define foreach BOOST_FOREACH +//#define DEBUG_SMD +#ifdef DEBUG_SMD +#define printf if (RT_DEBUG_SMD) printf +#else #define printf if (0) printf +#endif namespace { const double GOLD_R1 = 0.61803399 ; const double GOLD_R2 = 1-0.61803399 ; const double XTOL = 10*std::sqrt(std::numeric_limits::epsilon()); +#ifdef DEBUG_SMD + bool RT_DEBUG_SMD = 0; +#endif } bool OneDimMinimizer::minimize(int steps, double ytol, double xtol) @@ -32,8 +41,16 @@ bool OneDimMinimizer::minimize(int steps, double ytol, double xtol) OneDimMinimizer::ImproveRet OneDimMinimizer::improve(int steps, double ytol, double xtol, bool force) { - double x0 = var_->getVal(), y0 = nll_->getVal(); - xi_[1] = x0; + double x0 = var_->getVal(); + if (x0 < xi_[0] || x0 > xi_[2]) { + // could happen if somebody outside this loop modifies some parameters + printf("ODM: ALERT: variable %s outside bounds x = [%.4f, %.4f, %.4f], x0 = %.4f\n", var_->GetName(), xi_[0], xi_[1], xi_[2], x0); + x0 = xi_[1]; + var_->setVal(x0); + } else { + xi_[1] = x0; + } + double y0 = nll_->getVal(); yi_[1] = y0; yi_[0] = eval(xi_[0]); yi_[2] = eval(xi_[2]); @@ -54,10 +71,24 @@ OneDimMinimizer::ImproveRet OneDimMinimizer::improve(int steps, double ytol, dou reseek(); } + //post-condition: always a sorted interval + assert(xi_[0] < xi_[2]); + assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]); + // if midpoint is not not one of the extremes, it's not higher than that extreme + assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); + assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]); bool done = doloop(steps,ytol,xtol); parabolaStep(); - + + //post-condition: always a sorted interval + assert(xi_[0] < xi_[2]); + assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]); + // if midpoint is not not one of the extremes, it's not higher than that extreme + assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); + assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]); + + if (ytol > 0 && fabs(yi_[1] - y0) < ytol) { printf("ODM: final ytol for %s: ymin(old) %.8f, ymin(new) %.8f, diff %.8f\n", var_->GetName(), y0, yi_[1], y0 -yi_[1]); if (!force) var_->setVal(x0); @@ -69,7 +100,7 @@ OneDimMinimizer::ImproveRet OneDimMinimizer::improve(int steps, double ytol, dou return Unchanged; } printf("ODM: doloop for %s is %s\n", var_->GetName(), done ? "Done" : "NotDone"); - printf("ODM: start of improve %s x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", var_->GetName(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]); + printf("ODM: end of improve %s x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", var_->GetName(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]); var_->setVal(xi_[1]); var_->setError(xi_[2] - xi_[0]); return done ? Done : NotDone; @@ -116,6 +147,7 @@ void OneDimMinimizer::reseek() xi_[2] = std::min(xmax, xi_[1]+xerr); yi_[2] = (xi_[2] == xi_[1] ? yi_[1] : eval(xi_[2])); assert(xi_[0] < xi_[2]); + assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]); for (;;) { //printf("ODM: bracketing %s in x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", var_->GetName(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]); @@ -136,16 +168,37 @@ void OneDimMinimizer::reseek() xerr *= 2; } //printf("ODM: bracketed minimum of %s in [%.4f, %.4f, %.4f]\n", var_->GetName(), xi_[0], xi_[1], xi_[2]); + //post-condition: always a sorted interval assert(xi_[0] < xi_[2]); + assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]); + // if midpoint is not not one of the extremes, it's not higher than that extreme + assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); + assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]); var_->setError(xerr); } void OneDimMinimizer::goldenBisection() { + //pre-condition: always a sorted interval assert(xi_[0] < xi_[2]); + assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]); + // if midpoint is not not one of the extremes, it's not higher than that extreme + assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); + assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]); if (xi_[0] == xi_[1] || xi_[1] == xi_[2]) { + int isame = (xi_[0] == xi_[1] ? 0 : 2); + /// pre-condition: the endpoint equal to x1 is not the highest + assert(yi_[isame] <= yi_[2-isame]); xi_[1] = 0.5*(xi_[0]+xi_[2]); yi_[1] = eval(xi_[1]); + if (yi_[1] < yi_[isame]) { + // maximum is in the interval- + // leave as is, next bisection steps will find it + } else { + // maximum remains on the boundary, leave both points there + assign(2-isame, 1); + assign(1, isame); + } } else { int inear = 2, ifar = 0; if (xi_[2]-xi_[1] > xi_[1] - xi_[0]) { @@ -164,13 +217,19 @@ void OneDimMinimizer::goldenBisection() xi_[ifar] = xc; yi_[ifar] = yc; } } + //post-condition: always a sorted interval assert(xi_[0] < xi_[2]); + assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]); + // if midpoint is not not one of the extremes, it's not higher than that extreme + assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); + assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]); + } double OneDimMinimizer::parabolaFit() { if (xi_[0] == xi_[1] || xi_[1] == xi_[2]) { - return 0.5*(xi_[1]+xi_[2]); + return xi_[1]; } double dx0 = xi_[1] - xi_[0], dx2 = xi_[1] - xi_[2]; double dy0 = yi_[1] - yi_[0], dy2 = yi_[1] - yi_[2]; @@ -189,9 +248,15 @@ bool OneDimMinimizer::parabolaStep() { double xc = parabolaFit(); if (xc != xi_[1]) { double yc = eval(xc); - if (yc < xi_[1]) { + if (yc < yi_[1]) { xi_[1] = xc; yi_[1] = yc; + //post-condition: always a sorted interval + assert(xi_[0] < xi_[2]); + assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]); + // if midpoint is not not one of the extremes, it's not higher than that extreme + assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); + assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]); return true; } } @@ -200,13 +265,20 @@ bool OneDimMinimizer::parabolaStep() { double OneDimMinimizer::eval(double x) { + double x0 = var_->getVal(); var_->setVal(x); - return nll_->getVal(); + double y = nll_->getVal(); + var_->setVal(x0); + return y; } SequentialMinimizer::SequentialMinimizer(RooAbsReal *nll, RooRealVar *poi) : - nll_(nll) + nll_(nll), + state_(Cleared) { +#ifdef DEBUG_SMD + RT_DEBUG_SMD = runtimedef::get("DEBUG_SMD"); +#endif std::auto_ptr args(nll->getParameters((const RooArgSet *)0)); workers_.reserve(args->getSize()); if (poi != 0) workers_.push_back(Worker(nll,poi)); @@ -222,27 +294,42 @@ SequentialMinimizer::SequentialMinimizer(RooAbsReal *nll, RooRealVar *poi) : bool SequentialMinimizer::minimize(double ytol, int bigsteps, int smallsteps) { foreach(Worker &w, workers_) { - w.minimize(1); - w.state = Ready; + if (w.var().isConstant()) { + w.state = Fixed; + } else { + w.minimize(1); + w.state = Ready; + } } + state_ = Ready; return improve(ytol, bigsteps, smallsteps); } bool SequentialMinimizer::improve(double ytol, int bigsteps, int smallsteps) { + // catch improve before minimize case + if (state_ == Cleared) return minimize(ytol,bigsteps,smallsteps); + + // setup default tolerances and steps if (ytol == 0) ytol = ROOT::Math::MinimizerOptions::DefaultTolerance()/sqrt(workers_.size()); if (bigsteps == 0) bigsteps = 100 * (workers_.size()); - State state = Active; + + // list of done workers (latest-done on top) std::list doneWorkers; + + // start with active workers, for all except constants foreach(Worker &w, workers_) { - if (w.state == Done) doneWorkers.push_back(&w); + if (w.var().isConstant()) w.state = Fixed; + else w.state = Active; } + + state_ = Active; for (int i = 0; i < bigsteps; ++i) { - printf("Start of loop. State is %s\n",(state == Done ? "DONE" : "ACTIVE")); + printf("Start of loop. State is %s\n",(state_ == Done ? "DONE" : "ACTIVE")); State newstate = Done; foreach(Worker &w, workers_) { OneDimMinimizer::ImproveRet iret = OneDimMinimizer::Unchanged; - if (w.state == Done) continue; + if (w.state == Done || w.state == Fixed) continue; iret = w.improve(smallsteps,ytol); if (iret == OneDimMinimizer::Unchanged) { printf("\tMinimized %s: Unchanged. NLL = %.8f\n", w.var().GetName(), nll_->getVal()); @@ -251,7 +338,7 @@ bool SequentialMinimizer::improve(double ytol, int bigsteps, int smallsteps) } else { printf("\tMinimized %s: Changed. NLL = %.8f\n", w.var().GetName(), nll_->getVal()); w.state = Active; - newstate = Active; + newstate = Active; } } if (newstate == Done) { @@ -272,13 +359,13 @@ bool SequentialMinimizer::improve(double ytol, int bigsteps, int smallsteps) } } printf("End of loop. New state is %s\n",(newstate == Done ? "DONE" : "ACTIVE")); - if (state == Done && newstate == Done) { - std::cout << "Converged after " << i << " big steps" << std::endl; + if (state_ == Done && newstate == Done) { + //std::cout << "Converged after " << i << " big steps" << std::endl; return true; } - state = newstate; + state_ = newstate; } - std::cout << "Did not converge after " << bigsteps << " big steps" << std::endl; + //std::cout << "Did not converge after " << bigsteps << " big steps" << std::endl; return false; } diff --git a/trunk/HiggsAnalysis/CombinedLimit/test/unit/testSeqMinim.cxx b/trunk/HiggsAnalysis/CombinedLimit/test/unit/testSeqMinim.cxx index 7d0ddb3472126..952ae39e396c5 100644 --- a/trunk/HiggsAnalysis/CombinedLimit/test/unit/testSeqMinim.cxx +++ b/trunk/HiggsAnalysis/CombinedLimit/test/unit/testSeqMinim.cxx @@ -18,6 +18,8 @@ #include "HiggsAnalysis/CombinedLimit/interface/RooSimultaneousOpt.h" #include "HiggsAnalysis/CombinedLimit/interface/SequentialMinimizer.h" #include "HiggsAnalysis/CombinedLimit/interface/ProfiledLikelihoodRatioTestStatExt.h" +#include "HiggsAnalysis/CombinedLimit/interface/ProfilingTools.h" + RooWorkspace *w; @@ -69,11 +71,12 @@ double testProfile(Func f, RooAbsReal *nll, RooRealVar *poi, double &xmin, doubl -void runExternal(const char *file, const char *algo, double tol, const char *wsp, const char *datan, const char *mcn) { +void runExternal(const char *file, double mh, const char *algo, double tol, const char *wsp, const char *datan, const char *mcn) { + runtimedef::set("DEBUG_SMD",1); ProfileLikelihood::MinimizerSentry cfg(algo,tol); TFile *f = TFile::Open(file); if (f == 0) return; w = (RooWorkspace *) f->Get(wsp); if (w == 0) return; - if (w->var("MH")) w->var("MH")->setVal(115.); + if (w->var("MH")) w->var("MH")->setVal(mh); //w->Print("V"); RooAbsData *data = w->data(datan); if (data == 0) return; RooStats::ModelConfig *mc = (RooStats::ModelConfig *) w->genobj(mcn); if (mc == 0) return; @@ -93,9 +96,12 @@ void runExternal(const char *file, const char *algo, double tol, const char *wsp std::cout << "dnR = " << dnR << ", dnS = " << dnS << ", diff = " << dnR - dnS << std::endl; #else - double xminR = 1.0, xminS = 1.0, yminR, yminS, y1R, y1S, q0R, q0S; + + double xminR = 1.0, xminS = 1.0, yminR = 0, yminS, y1R = 0, y1S, q0R = 0, q0S; +#ifdef CMP nuis = clean; q0R = testProfile(runRooMin, &nll, poi, xminR, 0.0, yminR, y1R); +#endif nuis = clean; q0S = testProfile(runSeqMin, &nll, poi, xminS, 0.0, yminS, y1S); printf("RooFit minimizer: minimum at r = %.6f, NLL(min) = %+16.8f, NLL(0) = %+16.8f, q0 = %8.5f\n", xminR, yminR, y1R, q0R); @@ -120,8 +126,9 @@ int main(int argc, char **argv) { if (argc >= 2) { if (strstr(argv[1],"root")) { runExternal(argv[1], - argc >= 3 ? argv[2] : "Minuit2", - argc >= 4 ? atof(argv[2]) : 0.0001, + argc >= 3 ? atof(argv[2]) : 120., + argc >= 4 ? argv[3] : "Minuit2", + argc >= 5 ? atof(argv[4]) : 0.0001, "w", "data_obs", "ModelConfig");