Reference Manual

Go to the documentation of this file.
00001 #include "../interface/ProfiledLikelihoodRatioTestStatExt.h"
00002 #include "../interface/CascadeMinimizer.h"
00003 #include "../interface/CloseCoutSentry.h"
00004 #include "../interface/utils.h"
00005 #include <stdexcept>
00006 #include <RooRealVar.h>
00007 #include <RooMinimizer.h>
00008 #include <RooFitResult.h>
00009 #include <RooSimultaneous.h>
00010 #include <RooCategory.h>
00011 #include <RooRandom.h>
00012 #include <Math/MinimizerOptions.h>
00013 #include <RooStats/RooStatsUtils.h>
00014 #include "../interface/ProfilingTools.h"
00015 
00016 //---- Uncomment this and run with --perfCounters to get statistics of successful and failed fits
00017 //#define DEBUG_FIT_STATUS
00018 #ifdef DEBUG_FIT_STATUS
00019 #define  COUNT_ONE(x) PerfCounter::add(x);
00020 #else
00021 #define  COUNT_ONE(X) 
00022 #endif
00023 
00024 //---- Uncomment this and set some of these to 1 to get more debugging
00025 #if 0
00026 #define DBG(X,Z) if (X) { Z; }
00027 #define DBGV(X,Z) if (X>1) { Z; }
00028 #define DBG_TestStat_params 0 // ProfiledLikelihoodRatioTestStatOpt::Evaluate; 1 = dump nlls; 2 = dump params at each eval
00029 #define DBG_TestStat_NOFIT  0 // FIXME HACK: if set, don't profile the likelihood, just evaluate it
00030 #define DBG_PLTestStat_ctor 2 // dump parameters in c-tor
00031 #define DBG_PLTestStat_pars 2 // dump parameters in eval
00032 #define DBG_PLTestStat_fit  2 // dump fit result
00033 #define DBG_PLTestStat_main 2 // limited debugging (final NLLs, failed fits)
00034 #else
00035 #define DBG(X,Z) 
00036 #define DBGV(X,Z) 
00037 #endif
00038 
00039 //============================================================================================================================================
00040 ProfiledLikelihoodRatioTestStatOpt::ProfiledLikelihoodRatioTestStatOpt(
00041     const RooArgSet & observables,
00042     RooAbsPdf &pdfNull, RooAbsPdf &pdfAlt,
00043     const RooArgSet *nuisances,
00044     const RooArgSet & paramsNull, const RooArgSet & paramsAlt,
00045     int verbosity)
00046 : 
00047     pdfNull_(&pdfNull), pdfAlt_(&pdfAlt),
00048     paramsNull_(pdfNull_->getVariables()), 
00049     paramsAlt_(pdfAlt_->getVariables()), 
00050     verbosity_(verbosity)
00051 {
00052     snapNull_.addClone(paramsNull);
00053     snapAlt_.addClone(paramsAlt);
00054     DBG(DBG_TestStat_params, (std::cout << "Null snapshot" << pdfNull_->GetName())) DBG(DBG_TestStat_params, (snapNull_.Print("V")))
00055     DBG(DBG_TestStat_params, (std::cout << "Alt  snapshot" << pdfAlt_->GetName()))  DBG(DBG_TestStat_params, (snapAlt_.Print("V")))
00056     if (nuisances) nuisances_.addClone(*nuisances);
00057 }
00058 
00059 Double_t ProfiledLikelihoodRatioTestStatOpt::Evaluate(RooAbsData& data, RooArgSet& nullPOI)
00060 {
00061     *paramsNull_ = nuisances_;
00062     *paramsNull_ = snapNull_;
00063     *paramsNull_ = nullPOI;
00064         
00065     DBGV(DBG_TestStat_params, (std::cout << "Parameters of null pdf (pre fit)" << pdfNull_->GetName()))
00066     DBGV(DBG_TestStat_params, (paramsNull_->Print("V")))
00067 
00068     bool canKeepNullNLL = createNLL(*pdfNull_, data, nllNull_);
00069     double nullNLL = minNLL(nllNull_);
00070     if (!canKeepNullNLL) nllNull_.reset();
00071     
00072     *paramsAlt_ = nuisances_;
00073     *paramsAlt_ = snapAlt_;
00074         
00075     DBGV(DBG_TestStat_params, (std::cout << "Parameters of alt pdf " << pdfAlt_->GetName()))
00076     DBGV(DBG_TestStat_params, (paramsAlt_->Print("V")))
00077     bool canKeepAltNLL = createNLL(*pdfAlt_, data, nllAlt_);
00078     double altNLL = minNLL(nllAlt_);
00079     if (!canKeepAltNLL) nllAlt_.reset();
00080     
00081     DBG(DBG_TestStat_params, (printf("Pln: null = %+8.4f, alt = %+8.4f\n", nullNLL, altNLL)))
00082     return nullNLL-altNLL;
00083 }
00084 
00085 bool ProfiledLikelihoodRatioTestStatOpt::createNLL(RooAbsPdf &pdf, RooAbsData &data, std::auto_ptr<RooAbsReal> &nll_) 
00086 {
00087     if (typeid(pdf) == typeid(RooSimultaneousOpt)) {
00088         if (nll_.get() == 0) nll_.reset(pdf.createNLL(data, RooFit::Constrain(nuisances_)));
00089         else ((cacheutils::CachingSimNLL&)(*nll_)).setData(data);
00090         return true;
00091     } else {
00092         nll_.reset(pdf.createNLL(data, RooFit::Constrain(nuisances_)));
00093         return false;
00094     }
00095 }
00096 
00097 
00098 double ProfiledLikelihoodRatioTestStatOpt::minNLL(std::auto_ptr<RooAbsReal> &nll_) 
00099 {
00100 #if defined(DBG_TestStat_NOFIT) && (DBG_TestStat_NOFIT > 0)
00101     if (verbosity_ > 0) std::cout << "Profiling likelihood for pdf " << pdf.GetName() << std::endl;
00102     std::auto_ptr<RooAbsReal> nll_(pdf.createNLL(data, RooFit::Constrain(nuisances_)));
00103     return nll_->getVal();
00104 #endif
00105     RooMinimizer minim(*nll_);
00106     minim.setStrategy(0);
00107     minim.setPrintLevel(verbosity_-2);
00108     minim.minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str());
00109     if (verbosity_ > 1) {
00110         std::auto_ptr<RooFitResult> res(minim.save());
00111         res->Print("V");
00112     }
00113     return nll_->getVal();
00114 }
00115 
00116 //============================================================================================================================================
00117 ProfiledLikelihoodTestStatOpt::ProfiledLikelihoodTestStatOpt(
00118     const RooArgSet & observables,
00119     RooAbsPdf &pdf, 
00120     const RooArgSet *nuisances, 
00121     const RooArgSet & params,
00122     const RooArgSet & poi,
00123     const RooArgList &gobsParams,
00124     const RooArgList &gobs,
00125     int verbosity, 
00126     OneSidedness oneSided)
00127 :
00128     pdf_(&pdf),
00129     gobsParams_(gobsParams),
00130     gobs_(gobs),
00131     verbosity_(verbosity),
00132     oneSided_(oneSided)
00133 {
00134     DBG(DBG_PLTestStat_main, (std::cout << "Created for " << pdf.GetName() << "." << std::endl))
00135 
00136     params.snapshot(snap_,false);
00137     ((RooRealVar*)snap_.find(params.first()->GetName()))->setConstant(false);
00138     if (nuisances) { nuisances_.add(*nuisances); snap_.addClone(*nuisances, /*silent=*/true); }
00139     params_.reset(pdf_->getParameters(observables));
00140     DBG(DBG_PLTestStat_ctor, (std::cout << "Observables: " << std::endl)) DBG(DBG_PLTestStat_ctor, (observables.Print("V")))
00141     DBG(DBG_PLTestStat_ctor, (std::cout << "All params: " << std::endl))  DBG(DBG_PLTestStat_ctor, (params_->Print("V")))
00142     DBG(DBG_PLTestStat_ctor, (std::cout << "Snapshot: " << std::endl))    DBG(DBG_PLTestStat_ctor, (snap_.Print("V")))
00143     DBG(DBG_PLTestStat_ctor, (std::cout << "POI: " << std::endl))         DBG(DBG_PLTestStat_ctor, (poi_.Print("V")))
00144     RooLinkedListIter it = poi.iterator();
00145     for (RooAbsArg *a = (RooAbsArg*) it.Next(); a != 0; a = (RooAbsArg*) it.Next()) {
00146         // search for this poi in the parameters and in the snapshot
00147         RooAbsArg *ps = snap_.find(a->GetName());   
00148         RooAbsArg *pp = params_->find(a->GetName());
00149         if (pp == 0) { std::cerr << "WARNING: NLL does not depend on POI " << a->GetName() << ", cannot profile" << std::endl; continue; }
00150         if (ps == 0) { std::cerr << "WARNING: no snapshot for POI " << a->GetName() << ", cannot profile"  << std::endl; continue; }
00151         poi_.add(*ps);
00152         poiParams_.add(*pp);
00153     }
00154 }
00155 
00156 
00157 Double_t ProfiledLikelihoodTestStatOpt::Evaluate(RooAbsData& data, RooArgSet& /*nullPOI*/)
00158 {
00159     bool do_debug = runtimedef::get("DEBUG_PLTSO");
00160     //static bool do_rescan = runtimedef::get("PLTSO_FULL_RESCAN");
00161     DBG(DBG_PLTestStat_main, (std::cout << "Being evaluated on " << data.GetName() << std::endl))
00162 
00163     // Take snapshot of initial state, to restore it at the end 
00164     RooArgSet initialState; params_->snapshot(initialState);
00165 
00166     DBG(DBG_PLTestStat_pars, std::cout << "Being evaluated on " << data.GetName() << ": params before snapshot are " << std::endl)
00167     DBG(DBG_PLTestStat_pars, params_->Print("V"))
00168     // Initialize parameters
00169     *params_ = snap_;
00170 
00171     //set value of "globalConstrained" nuisances to the value of the corresponding global observable and set them constant
00172     for (int i=0; i<gobsParams_.getSize(); ++i) {
00173       RooRealVar *nuis = (RooRealVar*)gobsParams_.at(i);
00174       RooRealVar *gobs = (RooRealVar*)gobs_.at(i);
00175       nuis->setVal(gobs->getVal());
00176       nuis->setConstant();
00177     }
00178         
00179     DBG(DBG_PLTestStat_pars, std::cout << "Being evaluated on " << data.GetName() << ": params after snapshot are " << std::endl)
00180     DBG(DBG_PLTestStat_pars, params_->Print("V"))
00181     // Initialize signal strength (or the first parameter)
00182     RooRealVar *rIn = (RooRealVar *) poi_.first();
00183     RooRealVar *r   = (RooRealVar *) params_->find(rIn->GetName());
00184     bool canKeepNLL = createNLL(*pdf_, data);
00185 
00186     double initialR = rIn->getVal();
00187 
00188     // Perform unconstrained minimization (denominator)
00189     if (poi_.getSize() == 1) {
00190         r->setMin(0); if (initialR == 0 || (oneSided_ != oneSidedDef)) r->removeMax(); else r->setMax(1.1*initialR); 
00191         r->setVal(initialR == 0 ? 0.5 : 0.5*initialR); //best guess
00192         r->setConstant(false);
00193     } else {
00194         utils::setAllConstant(poiParams_,false);
00195     }
00196     DBG(DBG_PLTestStat_pars, (std::cout << "r In: ")) DBG(DBG_PLTestStat_pars, (rIn->Print(""))) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
00197     DBG(DBG_PLTestStat_pars, std::cout << "r before the fit: ") DBG(DBG_PLTestStat_pars, r->Print("")) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
00198 
00199     //std::cout << "PERFORMING UNCONSTRAINED FIT " << r->GetName() << " [ " << r->getMin() << " - " << r->getMax() << " ] "<< std::endl;
00200     double nullNLL = minNLL(/*constrained=*/false, r);
00201     double bestFitR = r->getVal();
00202 
00203     DBG(DBG_PLTestStat_pars, (std::cout << "r after the fit: ")) DBG(DBG_PLTestStat_pars, (r->Print(""))) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
00204     DBG(DBG_PLTestStat_pars, std::cout << "Was evaluated on " << data.GetName() << ": params before snapshot are " << std::endl)
00205     DBG(DBG_PLTestStat_pars, params_->Print("V"))
00206 
00207     // Prepare for constrained minimization (numerator)
00208     if (poi_.getSize() == 1) {
00209         r->setVal(initialR); 
00210         r->setConstant(true);
00211     } else {
00212         poiParams_.assignValueOnly(poi_);
00213         utils::setAllConstant(poiParams_,true);
00214     }
00215     double thisNLL = nullNLL;
00216     if (initialR == 0 || oneSided_ != oneSidedDef || bestFitR < initialR) { 
00217         // must do constrained fit (if there's something to fit besides XS)
00218         //std::cout << "PERFORMING CONSTRAINED FIT " << r->GetName() << " == " << r->getVal() << std::endl;
00219         thisNLL = (nuisances_.getSize() > 0 ? minNLL(/*constrained=*/true, r) : nll_->getVal());
00220         if (thisNLL - nullNLL < -0.02) { 
00221             DBG(DBG_PLTestStat_main, (printf("  --> constrained fit is better... will repeat unconstrained fit\n")))
00222             utils::setAllConstant(poiParams_,false);
00223             nullNLL = minNLL(/*constrained=*/false, r);
00224             bestFitR = r->getVal();
00225             if (bestFitR > initialR && oneSided_ == oneSidedDef) {
00226                 DBG(DBG_PLTestStat_main, (printf("   after re-fit, signal %7.4f > %7.4f, test statistics will be zero.\n", bestFitR, initialR)))
00227                 thisNLL = nullNLL;
00228             }
00229         }
00230         /* This piece of debug code was added to see if we were finding a local minimum at zero instead of the global minimum.
00231            It doesn't seem to be the case, however  
00232         if (do_rescan && fabs(thisNLL - nullNLL) < 0.2 && initialR == 0) {
00233             if (do_debug) printf("Doing full rescan. best fit r = %.4f, -lnQ = %.5f\n", bestFitR, thisNLL-nullNLL);
00234             for (double rx = 0; rx <= 5; rx += 0.01) {
00235                 r->setVal(rx); double y = nll_->getVal();
00236                 if (y < nullNLL) {
00237                     printf("%4.2f\t%8.5f\t<<==== ALERT\n", rx, y - nullNLL);
00238                 } else {
00239                     printf("%4.2f\t%8.5f\n", rx, y - nullNLL);
00240                 }
00241             }
00242         } */
00243         if (bestFitR > initialR && oneSided_ == signFlipDef) {
00244             DBG(DBG_PLTestStat_main, (printf("   fitted signal %7.4f > %7.4f, test statistics will be negative.\n", bestFitR, initialR)))
00245             std::swap(thisNLL, nullNLL);
00246         }
00247     } else {
00248         DBG(DBG_PLTestStat_main, (printf("   signal fit %7.4f > %7.4f: don't need to compute numerator\n", bestFitR, initialR)))
00249     }
00250 
00251     //Restore initial state, to avoid issues with ToyMCSampler
00252     *params_ = initialState;
00253 
00254     if (!canKeepNLL) nll_.reset();
00255 
00256     DBG(DBG_PLTestStat_main, (printf("\nNLLs:  num % 10.4f, den % 10.4f (signal %7.4f), test stat % 10.4f\n", thisNLL, nullNLL, bestFitR, thisNLL-nullNLL)))
00257     if (do_debug) printf("\nNLLs:  num % 10.4f, den % 10.4f (signal %7.4f), test stat % 10.4f\n", thisNLL, nullNLL, bestFitR, thisNLL-nullNLL);
00258     //return std::min(thisNLL-nullNLL, 0.);
00259     return thisNLL-nullNLL;
00260 }
00261 
00262 std::vector<Double_t> ProfiledLikelihoodTestStatOpt::Evaluate(RooAbsData& data, RooArgSet& /*nullPOI*/, const std::vector<Double_t> &rVals)
00263 {
00264     static bool do_debug = runtimedef::get("DEBUG_PLTSO");
00265     std::vector<Double_t> ret(rVals.size(), 0);
00266 
00267     DBG(DBG_PLTestStat_main, (std::cout << "Being evaluated on " << data.GetName() << std::endl))
00268 
00269     // Take snapshot of initial state, to restore it at the end 
00270     RooArgSet initialState; params_->snapshot(initialState);
00271 
00272     DBG(DBG_PLTestStat_pars, std::cout << "Being evaluated on " << data.GetName() << ": params before snapshot are " << std::endl)
00273     DBG(DBG_PLTestStat_pars, params_->Print("V"))
00274     // Initialize parameters
00275     *params_ = snap_;
00276 
00277     //set value of "globalConstrained" nuisances to the value of the corresponding global observable and set them constant
00278     for (int i=0; i<gobsParams_.getSize(); ++i) {
00279       RooRealVar *nuis = (RooRealVar*)gobsParams_.at(i);
00280       RooRealVar *gobs = (RooRealVar*)gobs_.at(i);
00281       nuis->setVal(gobs->getVal());
00282       nuis->setConstant();
00283     }
00284         
00285     DBG(DBG_PLTestStat_pars, std::cout << "Being evaluated on " << data.GetName() << ": params after snapshot are " << std::endl)
00286     DBG(DBG_PLTestStat_pars, params_->Print("V"))
00287     // Initialize signal strength
00288     RooRealVar *rIn = (RooRealVar *) poi_.first();
00289     RooRealVar *r   = (RooRealVar *) params_->find(rIn->GetName());
00290     bool canKeepNLL = createNLL(*pdf_, data);
00291 
00292     double initialR = rVals.back();
00293 
00294     // Perform unconstrained minimization (denominator)
00295     r->setMin(0); if (initialR == 0 || (oneSided_ != oneSidedDef)) r->removeMax(); else r->setMax(1.1*initialR); 
00296     r->setVal(initialR == 0 ? 0.5 : 0.5*initialR); //best guess
00297     r->setConstant(false);
00298     DBG(DBG_PLTestStat_pars, (std::cout << "r In: ")) DBG(DBG_PLTestStat_pars, (rIn->Print(""))) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
00299     DBG(DBG_PLTestStat_pars, std::cout << "r before the fit: ") DBG(DBG_PLTestStat_pars, r->Print("")) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
00300 
00301     if (do_debug) std::cout << "PERFORMING UNCONSTRAINED FIT " << r->GetName() << " [ " << r->getMin() << " - " << r->getVal() << " - " << r->getMax() << " ] "<< std::endl;
00302     double nullNLL = minNLL(/*constrained=*/false, r);
00303     double bestFitR = r->getVal();
00304     // Take snapshot of initial state, to restore it at the end 
00305     RooArgSet bestFitState; params_->snapshot(bestFitState);
00306 
00307 
00308     DBG(DBG_PLTestStat_pars, (std::cout << "r after the fit: ")) DBG(DBG_PLTestStat_pars, (r->Print(""))) DBG(DBG_PLTestStat_pars, std::cout << std::endl)
00309     DBG(DBG_PLTestStat_pars, std::cout << "Was evaluated on " << data.GetName() << ": params before snapshot are " << std::endl)
00310     DBG(DBG_PLTestStat_pars, params_->Print("V"))
00311 
00312     double EPS = 0.25*ROOT::Math::MinimizerOptions::DefaultTolerance();
00313     for (int iR = 0, nR = rVals.size(); iR < nR; ++iR) {
00314         if (fabs(ret[iR]) > 10*EPS) continue; // don't bother re-update points which were too far from zero anyway.
00315         *params_ = bestFitState;
00316         initialR = rVals[iR];
00317         // Prepare for constrained minimization (numerator)
00318         r->setVal(initialR); 
00319         r->setConstant(true);
00320         double thisNLL = nullNLL, sign = +1.0;
00321         if (initialR == 0 || oneSided_ != oneSidedDef || bestFitR < initialR) { 
00322             // must do constrained fit (if there's something to fit besides XS)
00323             //std::cout << "PERFORMING CONSTRAINED FIT " << r->GetName() << " == " << r->getVal() << std::endl;
00324             thisNLL = (nuisances_.getSize() > 0 ? minNLL(/*constrained=*/true, r) : nll_->getVal());
00325             if (thisNLL - nullNLL < 0 && thisNLL - nullNLL >= -EPS) {
00326                 thisNLL = nullNLL;
00327             } else if (thisNLL - nullNLL < 0) {
00328                 DBG(DBG_PLTestStat_main, (printf("  --> constrained fit is better... will repeat unconstrained fit\n")))
00329                 r->setConstant(false);
00330                 r->setVal(bestFitR);
00331                 double oldNullNLL = nullNLL;
00332                 nullNLL = minNLL(/*constrained=*/false, r);
00333                 bestFitR = r->getVal();
00334                 bestFitState.removeAll(); params_->snapshot(bestFitState);
00335                 for (int iR2 = 0; iR2 < iR; ++iR2) {
00336                     ret[iR2] -= (nullNLL - oldNullNLL); // fixup already computed test statistics
00337                 }
00338                 iR = -1; continue; // restart over again, refitting those close to zero :-(
00339             }
00340             if (bestFitR > initialR && oneSided_ == signFlipDef) {
00341                 DBG(DBG_PLTestStat_main, (printf("   fitted signal %7.4f > %7.4f, test statistics will be negative.\n", bestFitR, initialR)))
00342                 sign = -1.0;
00343             }
00344         } else {
00345             DBG(DBG_PLTestStat_main, (printf("   signal fit %7.4f > %7.4f: don't need to compute numerator\n", bestFitR, initialR)))
00346         }
00347 
00348         ret[iR] = sign * (thisNLL-nullNLL);
00349         DBG(DBG_PLTestStat_main, (printf("\nNLLs for %7.4f:  num % 10.4f, den % 10.4f (signal %7.4f), test stat % 10.4f\n", initialR, thisNLL, nullNLL, bestFitR, ret[iR])))
00350         if (do_debug) { 
00351             printf("   Q(%.4f) = %.4f (best fit signal %.4f), from num %.4f, den %.4f\n", rVals[iR], ret[iR], bestFitR, thisNLL, nullNLL); fflush(stdout);
00352         }
00353     }
00354     //Restore initial state, to avoid issues with ToyMCSampler
00355     *params_ = initialState;
00356 
00357     if (!canKeepNLL) nll_.reset();
00358 
00359     //return std::min(thisNLL-nullNLL, 0.);
00360     return ret;
00361 }
00362 bool ProfiledLikelihoodTestStatOpt::createNLL(RooAbsPdf &pdf, RooAbsData &data) 
00363 {
00364     if (typeid(pdf) == typeid(RooSimultaneousOpt)) {
00365         if (nll_.get() == 0) nll_.reset(pdf.createNLL(data, RooFit::Constrain(nuisances_)));
00366         else ((cacheutils::CachingSimNLL&)(*nll_)).setData(data);
00367         return true;
00368     } else {
00369         nll_.reset(pdf.createNLL(data, RooFit::Constrain(nuisances_)));
00370         return false;
00371     }
00372 }
00373 
00374 double ProfiledLikelihoodTestStatOpt::minNLL(bool constrained, RooRealVar *r) 
00375 {
00376     CascadeMinimizer::Mode mode(constrained ? CascadeMinimizer::Constrained : CascadeMinimizer::Unconstrained);
00377     CascadeMinimizer minim(*nll_, mode, r);
00378     minim.minimize(verbosity_-2);
00379     return nll_->getVal();
00380 }
00381 
00382 //============================================================ProfiledLikelihoodRatioTestStatExt
00383 bool nllutils::robustMinimize(RooAbsReal &nll, RooMinimizerOpt &minim, int verbosity) 
00384 {
00385     static bool do_debug = (runtimedef::get("DEBUG_MINIM") || runtimedef::get("DEBUG_PLTSO") > 1);
00386     double initialNll = nll.getVal();
00387     std::auto_ptr<RooArgSet> pars;
00388     bool ret = false;
00389     for (int tries = 0, maxtries = 4; tries <= maxtries; ++tries) {
00390         int status = minim.minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str());
00391         //if (verbosity > 1) res->Print("V");
00392         if (status == 0 && nll.getVal() > initialNll + 0.02) {
00393             std::auto_ptr<RooFitResult> res(minim.save());
00394             //PerfCounter::add("Minimizer.save() called for false minimum"); 
00395             DBG(DBG_PLTestStat_main, (printf("\n  --> false minimum, status %d, cov. quality %d, edm %10.7f, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, res->covQual(), res->edm(), initialNll, nll.getVal(), initialNll - nll.getVal())))
00396             if (pars.get() == 0) pars.reset(nll.getParameters((const RooArgSet*)0));
00397             *pars = res->floatParsInit();
00398             if (tries == 0) {
00399                 COUNT_ONE("nllutils::robustMinimize: false minimum (first try)")
00400                 DBG(DBG_PLTestStat_main, (printf("    ----> Doing a re-scan and re-trying\n")))
00401                 minim.minimize("Minuit2","Scan");
00402             } else if (tries == 1) {
00403                 COUNT_ONE("nllutils::robustMinimize: false minimum (second try)")
00404                 DBG(DBG_PLTestStat_main, (printf("    ----> Re-trying with strategy = 1\n")))
00405                 minim.setStrategy(1);
00406             } else if (tries == 2) {
00407                 COUNT_ONE("nllutils::robustMinimize: false minimum (third try)")
00408                 DBG(DBG_PLTestStat_main, (printf("    ----> Re-trying with strategy = 2\n")))
00409                 minim.setStrategy(2);
00410             } else  {
00411                 COUNT_ONE("nllutils::robustMinimize: false minimum (third try)")
00412                 DBG(DBG_PLTestStat_main, (printf("    ----> Last attempt: simplex method \n")))
00413                 status = minim.minimize("Minuit2","Simplex");
00414                 if (nll.getVal() < initialNll + 0.02) {
00415                     DBG(DBG_PLTestStat_main, (printf("\n  --> success: status %d, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, initialNll, nll.getVal(), initialNll - nll.getVal())))
00416                     if (do_debug) printf("\n  --> success: status %d, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, initialNll, nll.getVal(), initialNll - nll.getVal());
00417                     COUNT_ONE("nllutils::robustMinimize: final success")
00418                     ret = true;
00419                     break;
00420                 } else {
00421                     COUNT_ONE("nllutils::robustMinimize: final fail")
00422                     DBG(DBG_PLTestStat_main, (printf("\n  --> final fail: status %d, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, initialNll, nll.getVal(), initialNll - nll.getVal())))
00423                     if (do_debug) printf("\n  --> final fail: status %d, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, initialNll, nll.getVal(), initialNll - nll.getVal());
00424                     return false;
00425                 }
00426             }
00427         } else if (status == 0) {  
00428             DBG(DBG_PLTestStat_main, (printf("\n  --> success: status %d, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, initialNll, nll.getVal(), initialNll - nll.getVal())))
00429             if (do_debug) printf("\n  --> success: status %d, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, initialNll, nll.getVal(), initialNll - nll.getVal());
00430             COUNT_ONE("nllutils::robustMinimize: final success")
00431             ret = true;
00432             break;
00433         } else if (tries != maxtries) {
00434             std::auto_ptr<RooFitResult> res(do_debug ? minim.save() : 0);
00435             //PerfCounter::add("Minimizer.save() called for failed minimization"); 
00436             if (tries > 0 && minim.edm() < 0.05*ROOT::Math::MinimizerOptions::DefaultTolerance()) {
00437                 DBG(DBG_PLTestStat_main, (printf("\n  --> acceptable: status %d, edm %10.7f, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, res->edm(), initialNll, nll.getVal(), initialNll - nll.getVal())))
00438                 if (do_debug) printf("\n  --> acceptable: status %d, edm %10.7f, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, res->edm(), initialNll, nll.getVal(), initialNll - nll.getVal());
00439                 COUNT_ONE("nllutils::robustMinimize: accepting fit with bad status but good EDM")
00440                 COUNT_ONE("nllutils::robustMinimize: final success")
00441                 ret = true;
00442                 break;
00443             }
00444             DBG(DBG_PLTestStat_main, (printf("\n  --> partial fail: status %d, cov. quality %d, edm %10.7f, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, res->covQual(), res->edm(), initialNll, nll.getVal(), initialNll - nll.getVal())))
00445             if (tries == 1) {
00446                 COUNT_ONE("nllutils::robustMinimize: failed first attempt")
00447                 DBG(DBG_PLTestStat_main, (printf("    ----> Doing a re-scan first, and switching to strategy 1\n")))
00448                 minim.minimize("Minuit2","Scan");
00449                 minim.setStrategy(1);
00450             }
00451             if (tries == 2) {
00452                 COUNT_ONE("nllutils::robustMinimize: failed second attempt")
00453                 DBG(DBG_PLTestStat_main, (printf("    ----> trying with strategy = 2\n")))
00454                 minim.minimize("Minuit2","Scan");
00455                 minim.setStrategy(2);
00456             }
00457         } else {
00458             std::auto_ptr<RooFitResult> res(do_debug ? minim.save() : 0);
00459             DBG(DBG_PLTestStat_main, (printf("\n  --> final fail: status %d, cov. quality %d, edm %10.7f, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, res->covQual(), res->edm(), initialNll, nll.getVal(), initialNll - nll.getVal())))
00460             if (do_debug) printf("\n  --> final fail: status %d, cov. quality %d, edm %10.7f, nll initial % 10.4f, nll final % 10.4f, change %10.5f\n", status, res->covQual(), res->edm(), initialNll, nll.getVal(), initialNll - nll.getVal());
00461             COUNT_ONE("nllutils::robustMinimize: final fail")
00462         }
00463     }
00464     return ret;
00465 }
00466
/data/refman/pasoursint/CMSSW_5_3_6/src/HiggsAnalysis/CombinedLimit/src/ProfiledLikelihoodRatioTestStatExt.cc
