diff --git a/docs/part3/nonstandard.md b/docs/part3/nonstandard.md index fcc66ebe89e..ff76ae297af 100644 --- a/docs/part3/nonstandard.md +++ b/docs/part3/nonstandard.md @@ -1052,7 +1052,9 @@ In case you see an excess somewhere in your analysis, you can evaluate the look- To calculate the look-elsewhere effect for a single parameter (in this case the mass of the resonance), you can follow the instructions below. Note that these instructions assume you have a workspace that is parametric in your resonance mass $m$, otherwise you need to fit each background toy with separate workspaces. We will assume the local significance for your excess is $\sigma$. * Generate background-only toys `combine ws.root -M GenerateOnly --toysFrequentist -m 16.5 -t 100 --saveToys --expectSignal=0`. The output will be something like `higgsCombineTest.GenerateOnly.mH16.5.123456.root`. + * For each toy, calculate the significance for a predefined range (e.g $m\in [10,35]$ GeV) in steps suitable to the resolution (e.g. 1 GeV). For `toy_1` the procedure would be: `for i in $(seq 10 35); do combine ws.root -M Significance --redefineSignalPOI r --freezeParameters MH --setParameter MH=$i -n $i -D higgsCombineTest.GenerateOnly.mH16.5.123456.root:toys/toy_1`. Calculate the maximum significance over all of these mass points - call this $\sigma_{max}$. + * Count how many toys have a maximum significance larger than the local one for your observed excess. This fraction of toys with $\sigma_{max}>\sigma$ is the global p-value. You can find more tutorials on the LEE [here](https://indico.cern.ch/event/456547/contributions/1126036/attachments/1188691/1724680/20151117_comb_tutorial_Lee.pdf) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index 5fb28aab025..91da1be399d 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -23,6 +23,8 @@ class MultiDimFit : public FitterAlgoBase { } virtual void applyOptions(const boost::program_options::variables_map &vm) ; + friend class RandStartPt; + protected: virtual bool runSpecific(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint); @@ -67,6 +69,10 @@ class MultiDimFit : public FitterAlgoBase { static std::string robustHesseLoad_; static std::string robustHesseSave_; + static int pointsRandProf_; + static std::string setParameterRandomInitialValueRanges_; + static int randPointsSeed_; + static std::string saveSpecifiedFuncs_; static std::string saveSpecifiedNuis_; static std::string saveSpecifiedIndex_; @@ -83,6 +89,7 @@ class MultiDimFit : public FitterAlgoBase { static std::vector specifiedVals_; static RooArgList specifiedList_; static bool saveInactivePOI_; + static bool skipDefaultStart_; // initialize variables void initOnce(RooWorkspace *w, RooStats::ModelConfig *mc_s) ; @@ -95,6 +102,9 @@ class MultiDimFit : public FitterAlgoBase { void doStitch2D(RooWorkspace *w, RooAbsReal &nll) ; void doImpact(RooFitResult &res, RooAbsReal &nll) ; + std::map> getRangesDictFromInString(std::string) ; + + // utilities /// for each RooRealVar, set a range 'box' from the PL profiling all other parameters void doBox(RooAbsReal &nll, double cl, const char *name="box", bool commitPoints=true) ; diff --git a/interface/RandStartPt.h b/interface/RandStartPt.h new file mode 100644 index 00000000000..a8e60518ffc --- /dev/null +++ b/interface/RandStartPt.h @@ -0,0 +1,45 @@ +#ifndef HiggsAnalysis_CombinedLimit_RandStartPt_h +#define HiggsAnalaysis_CombinedLimit_RandStartPt_h + +#include "HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h" +#include "HiggsAnalysis/CombinedLimit/interface/MultiDimFit.h" + +#include +#include +#include "RooRealVar.h" +#include "RooAbsReal.h" +#include "RooArgSet.h" +#include "RooCategory.h" + +class RandStartPt { + private: + RooAbsReal& nll_; + std::vector &specifiedvars_; + std::vector &specifiedvals_; + bool skipdefaultstart_; + std::string parameterRandInitialValranges_; + int numrandpts_; + int verbosity_; + bool fastscan_; + bool hasmaxdeltaNLLforprof_; + float maxdeltaNLLforprof_; + std::vector &specifiednuis_; + std::vector &specifiedfuncnames_; + std::vector &specifiedfunc_; + std::vector &specifiedfuncvals_; + std::vector &specifiedcatnames_; + std::vector &specifiedcat_; + std::vector &specifiedcatvals_; + unsigned int nOtherFloatingPOI_; + public: + RandStartPt(RooAbsReal& nll, std::vector &specifiedvars, std::vector &specifiedvals, bool skipdefaultstart, std::string parameterRandInitialValranges, int numrandpts, int verbose, bool fastscan, bool hasmaxdeltaNLLforprof, float maxdeltaNLLforprof, std::vector &specifiednuis, std::vector &specifiedfuncnames, std::vector &specifiedfunc, std::vector &specifiedfuncvals, std::vector &specifiedcatnames, std::vector &specifiedcat, std::vector &specifiedcatvals, unsigned int nOtherFloatingPOI); + std::map> getRangesDictFromInString(std::string params_ranges_string_in); + std::vector> vectorOfPointsToTry (); + void commitBestNLLVal(unsigned int idx, float &nllVal, double &probVal); + void setProfPOIvalues(unsigned int startptIdx, std::vector> &nested_vector_of_wc_vals); + void setValSpecifiedObjs(); + void doRandomStartPt1DGridScan(double &xval, unsigned int poiSize, std::vector &poival, std::vector &poivars, std::unique_ptr ¶m, RooArgSet &snap, float &deltaNLL, double &nll_init, CascadeMinimizer &minimObj); + void doRandomStartPt2DGridScan(double &xval, double &yval, unsigned int poiSize, std::vector &poival, std::vector &poivars, std::unique_ptr ¶m, RooArgSet &snap, float &deltaNLL, double &nll_init, MultiDimFit::GridType gridType, double deltaX, double deltaY, CascadeMinimizer &minimObj); + +}; +#endif diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 683d9ff480e..eaa83498ffe 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -19,6 +19,7 @@ #include "../interface/utils.h" #include "../interface/RobustHesse.h" #include "../interface/ProfilingTools.h" +#include "../interface/RandStartPt.h" #include "../interface/CombineLogger.h" #include @@ -61,6 +62,9 @@ float MultiDimFit::centeredRange_ = -1.0; bool MultiDimFit::robustHesse_ = false; std::string MultiDimFit::robustHesseLoad_ = ""; std::string MultiDimFit::robustHesseSave_ = ""; +int MultiDimFit::pointsRandProf_ = 0; +int MultiDimFit::randPointsSeed_ = 0; +std::string MultiDimFit::setParameterRandomInitialValueRanges_; std::string MultiDimFit::saveSpecifiedFuncs_; @@ -80,6 +84,7 @@ std::vector MultiDimFit::specifiedVars_; std::vector MultiDimFit::specifiedVals_; RooArgList MultiDimFit::specifiedList_; bool MultiDimFit::saveInactivePOI_= false; +bool MultiDimFit::skipDefaultStart_ = false; MultiDimFit::MultiDimFit() : FitterAlgoBase("MultiDimFit specific options") @@ -103,15 +108,19 @@ MultiDimFit::MultiDimFit() : ("saveSpecifiedFunc", boost::program_options::value(&saveSpecifiedFuncs_)->default_value(""), "Save specified function values (default = none)") ("saveSpecifiedIndex", boost::program_options::value(&saveSpecifiedIndex_)->default_value(""), "Save specified indexes/discretes (default = none)") ("saveInactivePOI", boost::program_options::value(&saveInactivePOI_)->default_value(saveInactivePOI_), "Save inactive POIs in output (1) or not (0, default)") + ("skipDefaultStart", boost::program_options::value(&skipDefaultStart_)->default_value(skipDefaultStart_), "Do not include the default start point in list of points to fit") ("startFromPreFit", boost::program_options::value(&startFromPreFit_)->default_value(startFromPreFit_), "Start each point of the likelihood scan from the pre-fit values") - ("alignEdges", boost::program_options::value(&alignEdges_)->default_value(alignEdges_), "Align the grid points such that the endpoints of the ranges are included") - ("setParametersForGrid", boost::program_options::value(&setParametersForGrid_)->default_value(""), "Set the values of relevant physics model parameters. Give a comma separated list of parameter value assignments. Example: CV=1.0,CF=1.0") - ("saveFitResult", "Save RooFitResult to multidimfit.root") - ("out", boost::program_options::value(&out_)->default_value(out_), "Directory to put the diagnostics output file in") - ("robustHesse", boost::program_options::value(&robustHesse_)->default_value(robustHesse_), "Use a more robust calculation of the hessian/covariance matrix") - ("robustHesseLoad", boost::program_options::value(&robustHesseLoad_)->default_value(robustHesseLoad_), "Load the pre-calculated Hessian") - ("robustHesseSave", boost::program_options::value(&robustHesseSave_)->default_value(robustHesseSave_), "Save the calculated Hessian") - ; + ("alignEdges", boost::program_options::value(&alignEdges_)->default_value(alignEdges_), "Align the grid points such that the endpoints of the ranges are included") + ("setParametersForGrid", boost::program_options::value(&setParametersForGrid_)->default_value(""), "Set the values of relevant physics model parameters. Give a comma separated list of parameter value assignments. Example: CV=1.0,CF=1.0") + ("saveFitResult", "Save RooFitResult to multidimfit.root") + ("out", boost::program_options::value(&out_)->default_value(out_), "Directory to put the diagnostics output file in") + ("robustHesse", boost::program_options::value(&robustHesse_)->default_value(robustHesse_), "Use a more robust calculation of the hessian/covariance matrix") + ("robustHesseLoad", boost::program_options::value(&robustHesseLoad_)->default_value(robustHesseLoad_), "Load the pre-calculated Hessian") + ("robustHesseSave", boost::program_options::value(&robustHesseSave_)->default_value(robustHesseSave_), "Save the calculated Hessian") + ("pointsRandProf", boost::program_options::value(&pointsRandProf_)->default_value(pointsRandProf_), "Number of random start points to try for the profiled POIs") + ("randPointsSeed", boost::program_options::value(&randPointsSeed_)->default_value(randPointsSeed_), "Seed to use when generating random start points to try for the profiled POIs") + ("setParameterRandomInitialValueRanges", boost::program_options::value(&setParameterRandomInitialValueRanges_)->default_value(""), "Range from which to draw random start points for the profiled POIs. This range should be equal to or smaller than the max and min values for the profiled POIs. Does not override max/min ranges for the given POIs. E.g. usage: c1=-5,5:c2=-1,1") + ; } void MultiDimFit::applyOptions(const boost::program_options::variables_map &vm) @@ -140,6 +149,11 @@ void MultiDimFit::applyOptions(const boost::program_options::variables_map &vm) if (vm["floatOtherPOIs"].defaulted()) floatOtherPOIs_ = true; if (vm["saveInactivePOI"].defaulted()) saveInactivePOI_ = true; } else throw std::invalid_argument(std::string("Unknown algorithm: "+algo)); + if (pointsRandProf_ > 0) { + // Probably not the best way of doing this + // But right now the pointsRandProf_ option relies on saveInactivePOI_ being true to include the profiled POIs in specifiedVars_ + saveInactivePOI_ = true; + } fastScan_ = (vm.count("fastScan") > 0); squareDistPoiStep_ = (vm.count("squareDistPoiStep") > 0); skipInitialFit_ = (vm.count("skipInitialFit") > 0); @@ -186,6 +200,7 @@ bool MultiDimFit::runSpecific(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooS if (verbose <= 3) RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); bool doHesse = (algo_ == Singles || algo_ == Impact) || (saveFitResult_) ; if ( !skipInitialFit_){ + std::cout << "Doing initial fit: " << std::endl; res.reset(doFit(pdf, data, (doHesse ? poiList_ : RooArgList()), constrainCmdArg, (saveFitResult_ && !robustHesse_), 1, true, false)); if (!res.get()) { std::cout << "\n " < 2) throw std::logic_error("Don't know how to do a grid with more than 2 POIs."); double nll0 = nll.getVal(); - if (setParametersForGrid_ != "") { RooArgSet allParams(w->allVars()); allParams.add(w->allCats()); @@ -604,7 +618,10 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) //minim.setStrategy(minimizerStrategy_); std::unique_ptr params(nll.getParameters((const RooArgSet *)0)); RooArgSet snap; params->snapshot(snap); - //snap.Print("V"); + if (verbose > 1) { + std::cout << "Print snap: " << std::endl; + snap.Print("V"); + } // check if gridPoints are defined std::vector pointsPerPoi; @@ -620,8 +637,23 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } if (n == 1) { + if (verbose > 1){ + std::cout << "\nStarting n==1. The nll0 from initial fit: " << nll0 << std::endl; + } unsigned int points = pointsPerPoi.size() == 0 ? points_ : pointsPerPoi[0]; - + // Set seed for random points + if (pointsRandProf_ > 0) { + std::cout<<"\n************************************************************************************"<::max()) { lastPoint_ = points - 1; } + for (unsigned int i = 0; i < points; ++i) { if (i < firstPoint_) continue; if (i > lastPoint_) break; @@ -665,42 +698,33 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) *params = snap; poiVals_[0] = x; poiVars_[0]->setVal(x); - // now we minimize - nll.clearEvalErrorLog(); - deltaNLL_ = nll.getVal() - nll0; - if (nll.numEvalErrors() > 0) { - deltaNLL_ = 9990; - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetVal(); - } - Combine::commitPoint(true, /*quantile=*/0); - continue; - } - bool ok = fastScan_ || (hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_) || utils::countFloating(*params)==0 ? - true : - minim.minimize(verbose-1); - if (ok) { - deltaNLL_ = nll.getVal() - nll0; - double qN = 2*(deltaNLL_); - double prob = ROOT::Math::chisquared_cdf_c(qN, n+nOtherFloatingPoi_); - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetIndex(); - } - Combine::commitPoint(true, /*quantile=*/prob); - } - } + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + ////////////////// Rand starting points for each profiled POI to get best nll //////////////////////////////////////////////////////////////////////////// + ///////////////// The default behavior (i.e. no random start point) is incorporated within the function below //////////////////////////////////////////// + ///////////////// To retrieve only default start point usage, set _ to 0 or leave it unspecified. ///////////////////////////////////////// + RandStartPt randStartPt(nll, specifiedVars_, specifiedVals_, skipDefaultStart_, setParameterRandomInitialValueRanges_, _, verbose, fastScan_, hasMaxDeltaNLLForProf_, maxDeltaNLLForProf_, specifiedNuis_, specifiedFuncNames_, specifiedFunc_, specifiedFuncVals_, specifiedCatNames_, specifiedCat_, specifiedCatVals_, nOtherFloatingPoi_); + randStartPt.doRandomStartPt1DGridScan(x, n, poiVals_, poiVars_, params, snap, deltaNLL_, nll0, minim); + + } // End of the loop over scan points } else if (n == 2) { + if (verbose > 1){ + std::cout << "\nStarting n==2. The nll0 from initial fit: " << nll0 << std::endl; + } RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); - CloseCoutSentry sentry(verbose < 2); + //Set seed for random points + if (pointsRandProf_ > 0) { + std::cout<<"\n************************************************************************************"<GetName(), x, poiVars_[1]->GetName(), y); - } + //if (verbose && (ipoint % nprint == 0)) { + //fprintf(sentry.trueStdOut(), "Point %d/%d, (i,j) = (%d,%d), %s = %f, %s = %f\n", + // fprintf("Point %d/%d, (i,j) = (%d,%d), %s = %f, %s = %f\n", + // ipoint,nTotal, i,j, poiVars_[0]->GetName(), x, poiVars_[1]->GetName(), y); + //} + //Explicitly printing this out to allow users to monitor the progress + std::cout << "Point " << ipoint << "/" << nTotal << " " <<"(i,j)= "<<"("<GetName() << " = " << x <<" "<GetName() << " = " <setVal(x); poiVars_[1]->setVal(y); - nll.clearEvalErrorLog(); nll.getVal(); - if (nll.numEvalErrors() > 0) { - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetIndex(); - } - deltaNLL_ = 9999; Combine::commitPoint(true, /*quantile=*/0); - if (gridType_ == G3x3) { - for (int i2 = -1; i2 <= +1; ++i2) { - for (int j2 = -1; j2 <= +1; ++j2) { - if (i2 == 0 && j2 == 0) continue; - poiVals_[0] = x + 0.33333333*i2*deltaX; - poiVals_[1] = y + 0.33333333*j2*deltaY; - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetIndex(); - } - deltaNLL_ = 9999; Combine::commitPoint(true, /*quantile=*/0); - } - } - } - continue; - } - // now we minimize - bool skipme = hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_; - bool ok = fastScan_ || skipme ? true : minim.minimize(verbose-1); - if (ok) { - deltaNLL_ = nll.getVal() - nll0; - double qN = 2*(deltaNLL_); - double prob = ROOT::Math::chisquared_cdf_c(qN, n+nOtherFloatingPoi_); - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetIndex(); - } - Combine::commitPoint(true, /*quantile=*/prob); - } - if (gridType_ == G3x3) { - bool forceProfile = !fastScan_ && std::min(fabs(deltaNLL_ - 1.15), fabs(deltaNLL_ - 2.995)) < 0.5; - utils::CheapValueSnapshot center(*params); - double x0 = x, y0 = y; - for (int i2 = -1; i2 <= +1; ++i2) { - for (int j2 = -1; j2 <= +1; ++j2) { - if (i2 == 0 && j2 == 0) continue; - center.writeTo(*params); - x = x0 + 0.33333333*i2*deltaX; - y = y0 + 0.33333333*j2*deltaY; - poiVals_[0] = x; poiVars_[0]->setVal(x); - poiVals_[1] = y; poiVars_[1]->setVal(y); - nll.clearEvalErrorLog(); nll.getVal(); - if (nll.numEvalErrors() > 0) { - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetIndex(); - } - deltaNLL_ = 9999; Combine::commitPoint(true, /*quantile=*/0); - continue; - } - deltaNLL_ = nll.getVal() - nll0; - if (forceProfile || (!fastScan_ && std::min(fabs(deltaNLL_ - 1.15), fabs(deltaNLL_ - 2.995)) < 0.5)) { - minim.minimize(verbose-1); - deltaNLL_ = nll.getVal() - nll0; - } - double qN = 2*(deltaNLL_); - double prob = ROOT::Math::chisquared_cdf_c(qN, n+nOtherFloatingPoi_); - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetVal(); - } - for(unsigned int j=0; jgetIndex(); - } - Combine::commitPoint(true, /*quantile=*/prob); - } - } - } - } - } + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////// Rand starting points for each profiled POI to get best nll ///////////////////////////////////////////////////////////////////////// + ////////////////// The default behavior (i.e. no random start pt) is incorporated within the function below //////////////////////////////////////////// + ///////////////// To retrieve only default start point usage, set pointsRandProf_ to 0 or leave it unspecified. ///////////////////////////////// + + RandStartPt randStartPt(nll, specifiedVars_, specifiedVals_, skipDefaultStart_, setParameterRandomInitialValueRanges_, pointsRandProf_, verbose, fastScan_, hasMaxDeltaNLLForProf_, maxDeltaNLLForProf_, specifiedNuis_, specifiedFuncNames_, specifiedFunc_, specifiedFuncVals_, specifiedCatNames_, specifiedCat_, specifiedCatVals_, nOtherFloatingPoi_); + randStartPt.doRandomStartPt2DGridScan(x, y, n, poiVals_, poiVars_, params, snap, deltaNLL_, nll0, gridType_, deltaX, deltaY, minim); + } //End of loop over y scan points + } //End of loop over x scan points } else { // Use utils routine if n > 2 RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); @@ -1229,3 +1167,24 @@ void MultiDimFit::splitGridPoints(const std::string& s, std::vector> MultiDimFit::getRangesDictFromInString(std::string params_ranges_string_in) { + std::map> out_range_dict; + std::vector params_ranges_string_lst; + boost::split(params_ranges_string_lst, params_ranges_string_in, boost::is_any_of(":")); + for (UInt_t p = 0; p < params_ranges_string_lst.size(); ++p) { + std::vector params_ranges_string; + boost::split(params_ranges_string, params_ranges_string_lst[p], boost::is_any_of("=,")); + if (params_ranges_string.size() != 3) { + std::cout << "Error parsing expression : " << params_ranges_string_lst[p] << std::endl; + } + std::string wc_name =params_ranges_string[0]; + float lim_lo = atof(params_ranges_string[1].c_str()); + float lim_hi = atof(params_ranges_string[2].c_str()); + out_range_dict.insert({wc_name,{lim_lo,lim_hi}}); + } + return out_range_dict; +} + diff --git a/src/RandStartPt.cc b/src/RandStartPt.cc new file mode 100644 index 00000000000..621ab07f882 --- /dev/null +++ b/src/RandStartPt.cc @@ -0,0 +1,265 @@ +#include "HiggsAnalysis/CombinedLimit/interface/RandStartPt.h" +#include "HiggsAnalysis/CombinedLimit/interface/Combine.h" +#include "HiggsAnalysis/CombinedLimit/interface/utils.h" +#include "HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h" + +#include +#include +#include + +#include "TMath.h" +#include "TFile.h" +#include "RooArgSet.h" +#include "RooArgList.h" +#include "RooRealVar.h" +#include "RooMinimizer.h" +#include + +#include +#include +#include +#include + +RandStartPt::RandStartPt(RooAbsReal& nll, std::vector &specifiedvars, std::vector &specifiedvals, bool skipdefaultstart, std::string parameterRandInitialValranges, int numrandpts, int verbose, bool fastscan, bool hasmaxdeltaNLLforprof, float maxdeltaNLLforprof, std::vector &specifiednuis, std::vector &specifiedfuncnames, std::vector &specifiedfunc, std::vector &specifiedfuncvals, std::vector &specifiedcatnames, std::vector &specifiedcat, std::vector &specifiedcatvals, unsigned int nOtherFloatingPOI) : + nll_(nll), + specifiedvars_(specifiedvars), + specifiedvals_(specifiedvals), + skipdefaultstart_(skipdefaultstart), + parameterRandInitialValranges_(parameterRandInitialValranges), + numrandpts_(numrandpts), + verbosity_(verbose), + fastscan_(fastscan), + hasmaxdeltaNLLforprof_(hasmaxdeltaNLLforprof), + maxdeltaNLLforprof_(maxdeltaNLLforprof), + specifiednuis_(specifiednuis), + specifiedfuncnames_(specifiedfuncnames), + specifiedfunc_(specifiedfunc), + specifiedfuncvals_(specifiedfuncvals), + specifiedcatnames_(specifiedcatnames), + specifiedcat_(specifiedcat), + specifiedcatvals_(specifiedcatvals), + nOtherFloatingPOI_(nOtherFloatingPOI) + + {} + +std::vector> RandStartPt::vectorOfPointsToTry (){ + std::vector> wc_vals_vec_of_vec = {}; + int n_prof_params = specifiedvars_.size(); + + if(!skipdefaultstart_) { + std::vector default_start_pt_vec; + for (int prof_param_idx = 0; prof_param_idxgetVal()); + } + wc_vals_vec_of_vec.push_back(default_start_pt_vec); + } + + // Append the random points to the vector of points to try + float prof_start_pt_range_max = 20.0; // Default to 20 if we're not asking for custom ranges + std::map> rand_ranges_dict; + if (parameterRandInitialValranges_ != "") { + rand_ranges_dict = RandStartPt::getRangesDictFromInString(parameterRandInitialValranges_); + } + + for (int pt_idx = 0; pt_idx wc_vals_vec; + for (int prof_param_idx=0; prof_param_idxGetName()) != rand_ranges_dict.end()){ //if the random starting point range for this floating POI was supplied during runtime + float rand_range_lo = rand_ranges_dict[specifiedvars_[prof_param_idx]->GetName()][0]; + float rand_range_hi = rand_ranges_dict[specifiedvars_[prof_param_idx]->GetName()][1]; + prof_start_pt_range_max = std::max(abs(rand_range_lo),abs(rand_range_hi)); + } + else { //if the random starting point range for this floating POI was not supplied during runtime, set the default low to -20 and high to +20 + rand_ranges_dict.insert({specifiedvars_[prof_param_idx]->GetName(),{-20.0,20.0}}); + } + } + //Get a random number in the range [-prof_start_pt_range_max,prof_start_pt_range_max] + float rand_num = (rand()*2.0*prof_start_pt_range_max)/RAND_MAX - prof_start_pt_range_max; + wc_vals_vec.push_back(rand_num); + } + wc_vals_vec_of_vec.push_back(wc_vals_vec); + + } + + //Print vector of points to try + if (verbosity_ > 1) { + std::cout<<"List of points to try for : "<GetName() <<" = "<< val << std::endl; + index++; + } + } + } + return wc_vals_vec_of_vec; +} + +// Extract the ranges map from the input string +// // Assumes the string is formatted with colons like "poi_name1=lo_lim,hi_lim:poi_name2=lo_lim,hi_lim" +std::map> RandStartPt::getRangesDictFromInString(std::string params_ranges_string_in) { + std::map> out_range_dict; + std::vector params_ranges_string_lst; + boost::split(params_ranges_string_lst, params_ranges_string_in, boost::is_any_of(":")); + for (UInt_t p = 0; p < params_ranges_string_lst.size(); ++p) { + std::vector params_ranges_string; + boost::split(params_ranges_string, params_ranges_string_lst[p], boost::is_any_of("=,")); + if (params_ranges_string.size() != 3) { + std::cout << "Error parsing expression : " << params_ranges_string_lst[p] << std::endl; + } + std::string wc_name =params_ranges_string[0]; + float lim_lo = atof(params_ranges_string[1].c_str()); + float lim_hi = atof(params_ranges_string[2].c_str()); + out_range_dict.insert({wc_name,{lim_lo,lim_hi}}); + } + return out_range_dict; +} + +void RandStartPt::commitBestNLLVal(unsigned int idx, float &nllVal, double &probVal){//, RooAbsReal& nll_){ + if (idx==0){ + Combine::commitPoint(true, /*quantile=*/probVal); + nllVal = nll_.getVal(); + } else if (nll_.getVal() < nllVal){ + Combine::commitPoint(true, /*quantile=*/probVal); + nllVal = nll_.getVal(); + } +} + +void RandStartPt::setProfPOIvalues(unsigned int startptIdx, std::vector> &nested_vector_of_wc_vals){ + if (verbosity_ > 1) std::cout << "\n\tStart pt idx: " << startptIdx << std::endl; + for (unsigned int var_idx = 0; var_idx 1) std::cout << "\t\tThe var name: " << specifiedvars_[var_idx]->GetName() << std::endl; + if (verbosity_ > 1) std::cout << "\t\t\tRange before: " << specifiedvars_[var_idx]->getMin() << " " << specifiedvars_[var_idx]->getMax() << std::endl; + if (verbosity_ > 1) std::cout << "\t\t\t" << specifiedvars_[var_idx]->GetName() << " before setting: " << specifiedvars_[var_idx]->getVal() << " += " << specifiedvars_[var_idx]->getError() << std::endl; + specifiedvars_[var_idx]->setVal(nested_vector_of_wc_vals.at(startptIdx).at(var_idx)); + if (verbosity_ > 1) std::cout << "\t\t\tRange after: " << specifiedvars_[var_idx]->getMin() << " " << specifiedvars_[var_idx]->getMax() << std::endl; + if (verbosity_ > 1) std::cout << "\t\t\t" << specifiedvars_[var_idx]->GetName() << " after setting: " << specifiedvars_[var_idx]->getVal() << " += " << specifiedvars_[var_idx]->getError() << std::endl; + } +} + +void RandStartPt::setValSpecifiedObjs(){ + for(unsigned int j=0; jgetVal(); + } + for(unsigned int j=0; jgetVal(); + } + for(unsigned int j=0; jgetIndex(); + } +} + +void RandStartPt::doRandomStartPt1DGridScan(double &xval, unsigned int poiSize, std::vector &poival, std::vector &poivars, std::unique_ptr ¶m, RooArgSet &snap, float &deltaNLL, double &nll_init, CascadeMinimizer &minimObj){ + float current_best_nll = 0; + //the nested vector to hold random starting points to try + std::vector> nested_vector_of_wc_vals = vectorOfPointsToTry (); + for (unsigned int start_pt_idx = 0; start_pt_idxsetVal(xval); + + //Loop over prof POIs and set their values + setProfPOIvalues(start_pt_idx, nested_vector_of_wc_vals); + + //now we minimize + nll_.clearEvalErrorLog(); + deltaNLL = nll_.getVal() - nll_init; + if (nll_.numEvalErrors() > 0){ + deltaNLL = 9990; + setValSpecifiedObjs(); + Combine::commitPoint(true, /*quantile=*/0); + continue; + } + bool ok = fastscan_ || (hasmaxdeltaNLLforprof_ && (nll_.getVal() - nll_init) > maxdeltaNLLforprof_) || utils::countFloating(*param)==0 ? + true : + minimObj.minimize(verbosity_-1); + if (ok) { + deltaNLL = nll_.getVal() - nll_init; + double qN = 2*(deltaNLL); + double prob = ROOT::Math::chisquared_cdf_c(qN, poiSize + nOtherFloatingPOI_); + setValSpecifiedObjs(); + //finally, commit best NLL value + commitBestNLLVal(start_pt_idx, current_best_nll, prob); + } + } +} + +void RandStartPt::doRandomStartPt2DGridScan(double &xval, double &yval, unsigned int poiSize, std::vector &poival, std::vector &poivars, std::unique_ptr ¶m, RooArgSet &snap, float &deltaNLL, double &nll_init, MultiDimFit::GridType gridType, double deltaX, double deltaY, CascadeMinimizer &minimObj){ + float current_best_nll = 0; + //the nested vector to hold random starting points to try + std::vector> nested_vector_of_wc_vals = vectorOfPointsToTry (); + for (unsigned int start_pt_idx = 0; start_pt_idxsetVal(xval); + poivars[1]->setVal(yval); + + //Loop over prof POIs and set their values + setProfPOIvalues(start_pt_idx, nested_vector_of_wc_vals); + + //now we minimize + nll_.clearEvalErrorLog(); + nll_.getVal(); + deltaNLL = nll_.getVal() - nll_init; + if (nll_.numEvalErrors() > 0) { + setValSpecifiedObjs(); + deltaNLL = 9999; + Combine::commitPoint(true, /*quantile=*/0); + if (gridType == MultiDimFit::G3x3){ + for (int i2 = -1; i2 <= +1; ++i2){ + for (int j2 = -1; j2 <= +1; ++j2) { + if (i2 == 0 && j2 == 0) continue; + poival[0] = xval + 0.33333333*i2*deltaX; + poival[1] = yval + 0.33333333*j2*deltaY; + setValSpecifiedObjs(); + deltaNLL = 9999; Combine::commitPoint(true, /*quantile=*/0); + } + } + } + continue; + } + bool ok = fastscan_ || (hasmaxdeltaNLLforprof_ && (nll_.getVal() - nll_init) > maxdeltaNLLforprof_) ? + true : + minimObj.minimize(verbosity_-1); + if (ok) { + deltaNLL = nll_.getVal() - nll_init; + double qN = 2*(deltaNLL); + double prob = ROOT::Math::chisquared_cdf_c(qN, poiSize + nOtherFloatingPOI_); + setValSpecifiedObjs(); + commitBestNLLVal(start_pt_idx, current_best_nll, prob); + } + if (gridType == MultiDimFit::G3x3){ + bool forceProfile = !fastscan_ && std::min(fabs(deltaNLL - 1.15), fabs(deltaNLL - 2.995)) < 0.5; + utils::CheapValueSnapshot center(*param); + double x0 = xval, y0 = yval; + for (int i2 = -1; i2 <= +1; ++i2){ + for (int j2 = -1; j2 <= +1; ++j2){ + if (i2 == 0 && j2 == 0) continue; + center.writeTo(*param); + xval = x0 + 0.33333333*i2*deltaX; + yval = y0 + 0.33333333*j2*deltaY; + poival[0] = xval; poivars[0]->setVal(xval); + poival[1] = yval; poivars[1]->setVal(yval); + nll_.clearEvalErrorLog(); nll_.getVal(); + if (nll_.numEvalErrors() > 0){ + setValSpecifiedObjs(); + deltaNLL = 9999; Combine::commitPoint(true, /*quantile*/0); + continue; + } + deltaNLL = nll_.getVal() - nll_init; + if (forceProfile || (fastscan_ && std::min(fabs(deltaNLL - 1.15), fabs(deltaNLL - 2.995)) < 0.5)) { + minimObj.minimize(verbosity_-1); + deltaNLL = nll_.getVal() - nll_init; + } + double qN = 2*(deltaNLL); + double prob = ROOT::Math::chisquared_cdf_c(qN, poiSize + nOtherFloatingPOI_); + setValSpecifiedObjs(); + commitBestNLLVal(start_pt_idx, current_best_nll, prob); + } + } + } + } +} diff --git a/test/systematicsAnalyzer.py b/test/systematicsAnalyzer.py index 5f966d30dab..6e7326ebecc 100755 --- a/test/systematicsAnalyzer.py +++ b/test/systematicsAnalyzer.py @@ -114,6 +114,7 @@ options.optimizeTemplateBins = False options.forceNonSimPdf = False options.removeMultiPdf = False +options.noData = False options.physModel = "HiggsAnalysis.CombinedLimit.PhysicsModel:defaultModel" sys.argv = ["-b-"]