From 45fc551da052678ebff9624645d248ce0025cf6c Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 3 Dec 2021 16:52:18 -0500 Subject: [PATCH 01/37] Add ability to try rand start pts for prof POIs --- src/MultiDimFit.cc | 150 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 132 insertions(+), 18 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index ebf2d536633..7bad1ee0096 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -186,6 +186,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 " < 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; @@ -626,8 +630,12 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } if (n == 1) { + if (verbose > 1) std::cout << "The nll0 from initial fit: " << nll0 << std::endl; unsigned int points = pointsPerPoi.size() == 0 ? points_ : pointsPerPoi[0]; + // Set seed for random points + srand(1); + double xspacing = (pmax[0]-pmin[0]) / points; double xspacingOffset = 0.5; if (alignEdges_) { @@ -665,42 +673,148 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } //if (verbose > 1) std::cout << "Point " << i << "/" << points << " " << poiVars_[0]->GetName() << " = " << x << std::endl; - std::cout << "Point " << i << "/" << points << " " << poiVars_[0]->GetName() << " = " << x << std::endl; + std::cout << "Point " << i << "/" << points << " " << poiVars_[0]->GetName() << " = " << x << std::endl; + *params = snap; + poiVals_[0] = x; + poiVars_[0]->setVal(x); + + + /////////////////////////////////////////////////////////////////////////////////////////////////// + /////////////////// Loop over rand points for each profiled POI to get best nll /////////////////// + + bool ok; + int n_rand_start_pts_to_try = 0; + int n_prof_params = specifiedVars_.size(); + std::vector nll_at_alt_start_pts; + + // Get vector of points to try + float prof_start_pt_range_max = 20.0; // Is this what we want? + std::vector> wc_vals_vec_of_vec = {}; + + // Append the defualt start pt to the list of points to try + 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 vecotr of points to try + for (int pt_idx=0; pt_idx wc_vals_vec; + for (int prof_param_idx=0; prof_param_idx 1) { + std::cout << "List of points to try:" << std::endl; + for (auto vals_vec: wc_vals_vec_of_vec) { + std::cout << "\tThe vals at this point: " << std::endl; + for (auto val: vals_vec) { + std::cout << "\t\tPoint val: " << val << std::endl; + } + } + } + + // Loop over starting points to try for the prof WCs + for (unsigned int start_pt_idx=0; start_pt_idxsetVal(x); + + // Loop over prof POIs and set their values + if (verbose > 1) std::cout << "\n\tStart pt idx: " << start_pt_idx << std::endl; + for (unsigned int var_idx=0; var_idx 1) std::cout << "\t\tThe var name: " << specifiedVars_[var_idx]->GetName() << std::endl; + if (strcmp(specifiedVars_[var_idx]->GetName(),"r")==0) { + // Don't bother setting r to anything + if (verbose > 1) std::cout << "\t\t\tSkipping var " << specifiedVars_[var_idx]->GetName() << std::endl; + continue; + } + if (verbose > 1) std::cout << "\t\t\tRange before: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; + if (verbose > 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]->removeRange(); // Do not impose a range + specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(start_pt_idx).at(var_idx)); + if (verbose > 1) std::cout << "\t\t\tRange after: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; + if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " after setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; + } + + // 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; + } + ok = fastScan_ || (hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_) || utils::countFloating(*params)==0 ? + true : + minim.minimize(verbose-1); + + if (verbose > 1) std::cout << "\t\tThe nll.getVal() is: " << nll.getVal() << std::endl; + nll_at_alt_start_pts.push_back(nll.getVal()); + + } + + // Rerun the fit for the point with the min nll + int min_nll_idx; + min_nll_idx = std::min_element(nll_at_alt_start_pts.begin(),nll_at_alt_start_pts.end()) - nll_at_alt_start_pts.begin(); + if (verbose > 1) std::cout << "\tMin nll at idx: " << min_nll_idx << " " << nll_at_alt_start_pts.at(min_nll_idx) << std::endl; + if (verbose > 1) std::cout << "\tResetting for rerunning at the point that gives best nll:" << std::endl; *params = snap; poiVals_[0] = x; poiVars_[0]->setVal(x); + for (unsigned int var_idx=0; var_idx 1) std::cout << "\t\tSetting " << specifiedVars_[var_idx]->GetName() << " to " << wc_vals_vec_of_vec.at(min_nll_idx).at(var_idx) << std::endl; + specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(min_nll_idx).at(var_idx)); + } // 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(); - } + 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 ? + 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(); - } + for(unsigned int j=0; jgetVal(); + } + for(unsigned int j=0; jgetVal(); + } + for(unsigned int j=0; jgetIndex(); + } Combine::commitPoint(true, /*quantile=*/prob); } + + std::cout << "" << std::endl; } } else if (n == 2) { RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); From 3326195e5b203aaef8babe46f8d0273acb6ee255 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 3 Dec 2021 17:52:37 -0500 Subject: [PATCH 02/37] Option for specifying number of random start points --- interface/MultiDimFit.h | 2 ++ src/MultiDimFit.cc | 12 +++++++----- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index cf9e77a8353..6253e01b339 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -65,6 +65,8 @@ class MultiDimFit : public FitterAlgoBase { static std::string robustHesseLoad_; static std::string robustHesseSave_; + static int pointsRandProf_; + static std::string saveSpecifiedFuncs_; static std::string saveSpecifiedNuis_; static std::string saveSpecifiedIndex_; diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 7bad1ee0096..9e9eee53e79 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -60,6 +60,7 @@ float MultiDimFit::centeredRange_ = -1.0; bool MultiDimFit::robustHesse_ = false; std::string MultiDimFit::robustHesseLoad_ = ""; std::string MultiDimFit::robustHesseSave_ = ""; +int MultiDimFit::pointsRandProf_ = 0; std::string MultiDimFit::saveSpecifiedFuncs_; @@ -110,6 +111,7 @@ MultiDimFit::MultiDimFit() : ("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") ; } @@ -630,7 +632,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } if (n == 1) { - if (verbose > 1) std::cout << "The nll0 from initial fit: " << nll0 << std::endl; + 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 @@ -683,7 +685,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) /////////////////// Loop over rand points for each profiled POI to get best nll /////////////////// bool ok; - int n_rand_start_pts_to_try = 0; int n_prof_params = specifiedVars_.size(); std::vector nll_at_alt_start_pts; @@ -699,7 +700,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) wc_vals_vec_of_vec.push_back(default_start_pt_vec); // Append the random points to the vecotr of points to try - for (int pt_idx=0; pt_idx wc_vals_vec; for (int prof_param_idx=0; prof_param_idx 1) std::cout << "\t\t\tRange before: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; if (verbose > 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]->removeRange(); // Do not impose a range + if (pointsRandProf_ > 0) { + specifiedVars_[var_idx]->removeRange(); // Do not impose a range if we're trying multiple random start points + } specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(start_pt_idx).at(var_idx)); if (verbose > 1) std::cout << "\t\t\tRange after: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " after setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; @@ -814,7 +817,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) Combine::commitPoint(true, /*quantile=*/prob); } - std::cout << "" << std::endl; } } else if (n == 2) { RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); From cf60fa46fb2c9dd28e134fce81ecec2f82d3ab7a Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 3 Dec 2021 18:11:54 -0500 Subject: [PATCH 03/37] If using pointsRandProf_, turn on saveInactivePOI_ --- src/MultiDimFit.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 9e9eee53e79..0a1ceb40591 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -141,6 +141,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); From 4c10a7c065136dda5b909593b42b826a7b94646b Mon Sep 17 00:00:00 2001 From: nckw Date: Tue, 31 May 2022 16:22:50 +0200 Subject: [PATCH 04/37] Update nonstandard.md --- docs/part3/nonstandard.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/docs/part3/nonstandard.md b/docs/part3/nonstandard.md index 7284dfb888d..01a267d3bc7 100644 --- a/docs/part3/nonstandard.md +++ b/docs/part3/nonstandard.md @@ -1041,3 +1041,16 @@ The example given here is extremely basic and it should be noted that additiona !!! danger If trying to implement parametric uncertainties in this setup (eg on transfer factors) which are correlated with other channels and implemented separately, you ***MUST*** normalise the uncertainty effect so that the datacard line can read `param name X 1`. That is the uncertainty on this parameter must be 1. Without this, there will be inconsistency with other nuisances of the same name in other channels implemented as **shape** or **lnN**. + + +## Look-elsewhere effect for one parameter + +In case you see an excess somewhere in your analysis, you can evaluate the look-elsewhere effect (LEE) of that excess. For an explanation of the LEE, take a look at the CMS Statistics Committee Twiki [here](https://twiki.cern.ch/twiki/bin/viewauth/CMS/LookElsewhereEffect). + +To calculate the look-elsewhere effect for a single parameter (in this case the mass of the resonance), you can follow the instructions. Note that these instructions assume you have a workspace which is parametric in your resonance mass $m$, otherwise you need to fit each background toy with separate workspaces. 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 - eg 1 GeV. `for i in $(seq 10 35); do combine ws.root -M Significance --redefineSignalPOI r --freezeParameters MH --setParameter MH=$i -n $i`. 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) From 018fe0e76757e1a6c0c714a5d558ebb378173b77 Mon Sep 17 00:00:00 2001 From: nckw Date: Tue, 31 May 2022 16:25:21 +0200 Subject: [PATCH 05/37] Update nonstandard.md --- docs/part3/nonstandard.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/part3/nonstandard.md b/docs/part3/nonstandard.md index 01a267d3bc7..8eb629a502e 100644 --- a/docs/part3/nonstandard.md +++ b/docs/part3/nonstandard.md @@ -1051,6 +1051,6 @@ To calculate the look-elsewhere effect for a single parameter (in this case the * 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 - eg 1 GeV. `for i in $(seq 10 35); do combine ws.root -M Significance --redefineSignalPOI r --freezeParameters MH --setParameter MH=$i -n $i`. 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. + * 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) From eff8006b3dc29afaa67917c057feb5de999bb31d Mon Sep 17 00:00:00 2001 From: nckw Date: Fri, 10 Jun 2022 12:01:55 +0100 Subject: [PATCH 06/37] Update commonstatsmethods.md --- docs/part3/commonstatsmethods.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/part3/commonstatsmethods.md b/docs/part3/commonstatsmethods.md index 2842735ab73..94215be353e 100644 --- a/docs/part3/commonstatsmethods.md +++ b/docs/part3/commonstatsmethods.md @@ -332,8 +332,8 @@ The choice of test-statistic can be made via the option `--testStat` and differe * **LHC-style**: `--LHCmode LHC-limits` , which is the shortcut for `--testStat LHC --generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1` - * The test statistic is defined using the ratio of likelihoods $q_{r} = -2\ln[\mathcal{L}(\mathrm{data}|r,\hat{\theta}_{r})/\mathcal{L}(\mathrm{data}|r=\hat{r},\hat{\theta}])$ , in which the nuisance parameters are profiled separately for $r=\hat{r}$ and $r$. - * The value of $q_{r}$ set to 0 when $\hat{r}>r$ giving a one sided limit. Furthermore, the constraint $r>0$ is enforced in the fit. This means that if the unconstrained value of $\hat{r}$ would be negative, the test statistic $q_{r}$ is evaluated as $-2\ln[\mathcal{L}(\mathrm{data}|r,\hat{\theta}_{r})/\mathcal{L}(\mathrm{data}|r=0,\hat{\theta}_{0}])$ + * The test statistic is defined using the ratio of likelihoods $q_{r} = -2\ln[\mathcal{L}(\mathrm{data}|r,\hat{\theta}_{r})/\mathcal{L}(\mathrm{data}|r=\hat{r},\hat{\theta})]$ , in which the nuisance parameters are profiled separately for $r=\hat{r}$ and $r$. + * The value of $q_{r}$ set to 0 when $\hat{r}>r$ giving a one sided limit. Furthermore, the constraint $r>0$ is enforced in the fit. This means that if the unconstrained value of $\hat{r}$ would be negative, the test statistic $q_{r}$ is evaluated as $-2\ln[\mathcal{L}(\mathrm{data}|r,\hat{\theta}_{r})/\mathcal{L}(\mathrm{data}|r=0,\hat{\theta}_{0})]$ * For the purposes of toy generation, the nuisance parameters are fixed to their **post-fit** values from the data (conditionally on the value of **r**), while the constraint terms are randomized in the evaluation of the likelihood. !!! warning From 032dc174be6b0159c78f7f40ceb9139d64e8c809 Mon Sep 17 00:00:00 2001 From: nckw Date: Fri, 10 Jun 2022 14:47:24 +0100 Subject: [PATCH 07/37] Update index.md Update standalone to 8.2.0 --- docs/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/index.md b/docs/index.md index e86c240cb61..a6b11bb196f 100644 --- a/docs/index.md +++ b/docs/index.md @@ -74,7 +74,7 @@ Access to `/cvmfs/cms.cern.ch/` can be obtained from lxplus machines or via `Cer git clone https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit.git HiggsAnalysis/CombinedLimit cd HiggsAnalysis/CombinedLimit/ git fetch origin -git checkout v8.1.0 +git checkout v8.2.0 . env_standalone.sh make ``` From 6ed9c9a9a4b066eb0efb4615d208b8b044939a08 Mon Sep 17 00:00:00 2001 From: nckw Date: Tue, 14 Jun 2022 14:53:28 +0100 Subject: [PATCH 08/37] Update commonstatsmethods.md --- docs/part3/commonstatsmethods.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/part3/commonstatsmethods.md b/docs/part3/commonstatsmethods.md index 94215be353e..e203b6098c2 100644 --- a/docs/part3/commonstatsmethods.md +++ b/docs/part3/commonstatsmethods.md @@ -7,7 +7,9 @@ This section will assume that you are using the default model unless otherwise s ## Asymptotic Frequentist Limits The `AsymptoticLimits` method allows to compute quickly an estimate of the observed and expected limits, which is fairly accurate when the event yields are not too small and the systematic uncertainties don't play a major role in the result. -The limit calculation relies on an asymptotic approximation of the distributions of the **LHC** test-statistic, which is based on a profile likelihood ratio, under signal and background hypotheses to compute two p-values $p_{\mu}, p_{b}$ and therefore $CL_s=p_{\mu}/(1-p_{b})$ (see the (see the [FAQ](http://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part4/usefullinks/#faq) section for a description of these) - i.e it is the asymptotic approximation of computing limits with frequentist toys. +The limit calculation relies on an asymptotic approximation of the distributions of the **LHC** test-statistic, which is based on a profile likelihood ratio, under signal and background hypotheses to compute two p-values $p_{\mu}, p_{b}$ and therefore $CL_s=p_{\mu}/(1-p_{b})$ (see the (see the [FAQ](http://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part4/usefullinks/#faq) section for a description of these) - i.e it is the asymptotic approximation of computing limits with frequentist toys using the LHC test-statistic for limits: + + * The test statistic is defined using the ratio of likelihoods $q_{r} = -2\ln[\mathcal{L}(\mathrm{data}|r,\hat{\theta}_{r})/\mathcal{L}(\mathrm{data}|r=\hat{r},\hat{\theta})]$ , in which the nuisance parameters are profiled separately for $r=\hat{r}$ and $r$. The value of $q_{r}$ set to 0 when $\hat{r}>r$ giving a one sided limit. Furthermore, the constraint $r>0$ is enforced in the fit. This means that if the unconstrained value of $\hat{r}$ would be negative, the test statistic $q_{r}$ is evaluated as $-2\ln[\mathcal{L}(\mathrm{data}|r,\hat{\theta}_{r})/\mathcal{L}(\mathrm{data}|r=0,\hat{\theta}_{0})]$ This method is so commonly used that it is the default method (i.e not specifying `-M` will run `AsymptoticLimits`) From 8feba2d0ad6fe154f93da46fe5170b3539279798 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Tue, 5 Jul 2022 00:11:01 -0400 Subject: [PATCH 09/37] Do not save nll if fit fails, also add ability to pass ranges for rand start points --- interface/MultiDimFit.h | 4 ++++ src/MultiDimFit.cc | 39 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index 6253e01b339..43475d1c2c9 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -66,6 +66,7 @@ class MultiDimFit : public FitterAlgoBase { static std::string robustHesseSave_; static int pointsRandProf_; + static std::string randPointsRanges_; static std::string saveSpecifiedFuncs_; static std::string saveSpecifiedNuis_; @@ -95,6 +96,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/src/MultiDimFit.cc b/src/MultiDimFit.cc index 0a1ceb40591..6f99cabae0e 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -61,6 +61,7 @@ bool MultiDimFit::robustHesse_ = false; std::string MultiDimFit::robustHesseLoad_ = ""; std::string MultiDimFit::robustHesseSave_ = ""; int MultiDimFit::pointsRandProf_ = 0; +std::string MultiDimFit::randPointsRanges_; std::string MultiDimFit::saveSpecifiedFuncs_; @@ -112,6 +113,7 @@ MultiDimFit::MultiDimFit() : ("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") + ("randPointsRanges", boost::program_options::value(&randPointsRanges_)->default_value(""), "Range from which to draw random start points for the profiled POIs") ; } @@ -694,7 +696,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) std::vector nll_at_alt_start_pts; // Get vector of points to try - float prof_start_pt_range_max = 20.0; // Is this what we want? std::vector> wc_vals_vec_of_vec = {}; // Append the defualt start pt to the list of points to try @@ -705,9 +706,19 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) wc_vals_vec_of_vec.push_back(default_start_pt_vec); // Append the random points to the vecotr 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 (randPointsRanges_ != "") { + rand_ranges_dict = getRangesDictFromInString(randPointsRanges_); + } for (int pt_idx=0; pt_idx wc_vals_vec; for (int prof_param_idx=0; prof_param_idxGetName()][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)); + } // 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); @@ -770,7 +781,11 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) minim.minimize(verbose-1); if (verbose > 1) std::cout << "\t\tThe nll.getVal() is: " << nll.getVal() << std::endl; - nll_at_alt_start_pts.push_back(nll.getVal()); + if (ok) { + nll_at_alt_start_pts.push_back(nll.getVal()); + } else { + nll_at_alt_start_pts.push_back(999999.9); // Placeholder for now + } } @@ -1354,3 +1369,23 @@ 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; +} From 2f19e028fb8cdf11b3578d093127c03b90816ed2 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Tue, 5 Jul 2022 01:08:53 -0400 Subject: [PATCH 10/37] Add option for skipping default start point --- interface/MultiDimFit.h | 1 + src/MultiDimFit.cc | 12 ++++++++---- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index 43475d1c2c9..65dded8e250 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -84,6 +84,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) ; diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 6f99cabae0e..36056222552 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -81,6 +81,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") @@ -104,6 +105,7 @@ 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") @@ -699,11 +701,13 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) std::vector> wc_vals_vec_of_vec = {}; // Append the defualt start pt to the list of points to try - std::vector default_start_pt_vec; - for (int prof_param_idx=0; prof_param_idxgetVal()); + 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); } - wc_vals_vec_of_vec.push_back(default_start_pt_vec); // Append the random points to the vecotr of points to try float prof_start_pt_range_max = 20.0; // Default to 20 if we're not asking for custom ranges From e1f2b668024272bef06cd132d43a58b996c5cbaf Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Tue, 5 Jul 2022 23:07:08 -0400 Subject: [PATCH 11/37] Save each point if best nll so far to avoid rerunning fit at end --- src/MultiDimFit.cc | 81 ++++++++++++++-------------------------------- 1 file changed, 25 insertions(+), 56 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 36056222552..f068dc531a4 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -742,6 +742,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } // Loop over starting points to try for the prof WCs + float current_best_nll; for (unsigned int start_pt_idx=0; start_pt_idx maxDeltaNLLForProf_) || utils::countFloating(*params)==0 ? - true : + ok = fastScan_ || (hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_) || utils::countFloating(*params)==0 ? + true : minim.minimize(verbose-1); - if (verbose > 1) std::cout << "\t\tThe nll.getVal() is: " << nll.getVal() << std::endl; if (ok) { - nll_at_alt_start_pts.push_back(nll.getVal()); - } else { - nll_at_alt_start_pts.push_back(999999.9); // Placeholder for now - } - - } - - // Rerun the fit for the point with the min nll - int min_nll_idx; - min_nll_idx = std::min_element(nll_at_alt_start_pts.begin(),nll_at_alt_start_pts.end()) - nll_at_alt_start_pts.begin(); - if (verbose > 1) std::cout << "\tMin nll at idx: " << min_nll_idx << " " << nll_at_alt_start_pts.at(min_nll_idx) << std::endl; - if (verbose > 1) std::cout << "\tResetting for rerunning at the point that gives best nll:" << std::endl; - *params = snap; - poiVals_[0] = x; - poiVars_[0]->setVal(x); - for (unsigned int var_idx=0; var_idx 1) std::cout << "\t\tSetting " << specifiedVars_[var_idx]->GetName() << " to " << wc_vals_vec_of_vec.at(min_nll_idx).at(var_idx) << std::endl; - specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(min_nll_idx).at(var_idx)); - } - // 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(); + 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 (start_pt_idx==0) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } else if (nll.getVal() < current_best_nll) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } } - Combine::commitPoint(true, /*quantile=*/0); - continue; - } - 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); - } + } // End loop over random start points - } + } // End of the loop over scan points } else if (n == 2) { RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); CloseCoutSentry sentry(verbose < 2); From 317704142071fee4bb97ce9d27f7617b375d1c15 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Wed, 6 Jul 2022 16:15:34 -0400 Subject: [PATCH 12/37] Pass rand seed for rand start points --- interface/MultiDimFit.h | 1 + src/MultiDimFit.cc | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index 65dded8e250..eb398032bc0 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -67,6 +67,7 @@ class MultiDimFit : public FitterAlgoBase { static int pointsRandProf_; static std::string randPointsRanges_; + static int randPointsSeed_; static std::string saveSpecifiedFuncs_; static std::string saveSpecifiedNuis_; diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index f068dc531a4..ddba3f93b4a 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -61,6 +61,7 @@ bool MultiDimFit::robustHesse_ = false; std::string MultiDimFit::robustHesseLoad_ = ""; std::string MultiDimFit::robustHesseSave_ = ""; int MultiDimFit::pointsRandProf_ = 0; +int MultiDimFit::randPointsSeed_ = 0; std::string MultiDimFit::randPointsRanges_; @@ -115,6 +116,7 @@ MultiDimFit::MultiDimFit() : ("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") ("randPointsRanges", boost::program_options::value(&randPointsRanges_)->default_value(""), "Range from which to draw random start points for the profiled POIs") ; } @@ -645,7 +647,9 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) unsigned int points = pointsPerPoi.size() == 0 ? points_ : pointsPerPoi[0]; // Set seed for random points - srand(1); + if (pointsRandProf_ > 0) { + srand(randPointsSeed_); + } double xspacing = (pmax[0]-pmin[0]) / points; double xspacingOffset = 0.5; From 2c1e05b19f45daddacb2002d239087430690999c Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 8 Jul 2022 15:08:54 -0400 Subject: [PATCH 13/37] Small cleanups --- src/MultiDimFit.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index ddba3f93b4a..14595ef538e 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -699,7 +699,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) bool ok; int n_prof_params = specifiedVars_.size(); - std::vector nll_at_alt_start_pts; // Get vector of points to try std::vector> wc_vals_vec_of_vec = {}; @@ -746,7 +745,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } // Loop over starting points to try for the prof WCs - float current_best_nll; + float current_best_nll = 0; // This variable will be properly initialized within the for loop, just initialize to a dummy value here to avoid compiling warnings for (unsigned int start_pt_idx=0; start_pt_idx Date: Thu, 21 Jul 2022 17:53:43 +0100 Subject: [PATCH 14/37] Update systematicsAnalyzer.py Include default for "noData" option (avoid complaint when running with --t2w) --- test/systematicsAnalyzer.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/systematicsAnalyzer.py b/test/systematicsAnalyzer.py index 07f8f87f8b5..1990a78fd5d 100755 --- a/test/systematicsAnalyzer.py +++ b/test/systematicsAnalyzer.py @@ -44,6 +44,7 @@ options.optimizeTemplateBins=False options.forceNonSimPdf = False options.removeMultiPdf = False +options.noData = False options.physModel = "HiggsAnalysis.CombinedLimit.PhysicsModel:defaultModel" # import ROOT with a fix to get batch mode (http://root.cern.ch/phpBB3/viewtopic.php?t=3198) import sys From a22bb471eeae992f7e1d02f909f28c79aed01730 Mon Sep 17 00:00:00 2001 From: nckw Date: Tue, 15 Nov 2022 17:14:04 +0000 Subject: [PATCH 15/37] Include ATLAS ref for Impacts --- docs/part3/nonstandard.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/part3/nonstandard.md b/docs/part3/nonstandard.md index 8eb629a502e..81701bc9c38 100644 --- a/docs/part3/nonstandard.md +++ b/docs/part3/nonstandard.md @@ -227,7 +227,7 @@ A similar option, `--X-nuisance-group-function`, can be used to scale whole grou ## Nuisance parameter impacts -The impact of a nuisance parameter (NP) θ on a parameter of interest (POI) μ is defined as the shift Δμ that is induced as θ is fixed and brought to its +1σ or −1σ post-fit values, with all other parameters profiled as normal. +The impact of a nuisance parameter (NP) θ on a parameter of interest (POI) μ is defined as the shift Δμ that is induced as θ is fixed and brought to its +1σ or −1σ post-fit values, with all other parameters profiled as normal (see [JHEP 01 (2015) 069](https://link.springer.com/article/10.1007/JHEP01(2015)069) for a description of this method). This is effectively a measure of the correlation between the NP and the POI, and is useful for determining which NPs have the largest effect on the POI uncertainty. From a17df940fe21bc25f26737d666c76be0c5d6b39d Mon Sep 17 00:00:00 2001 From: nckw Date: Mon, 28 Nov 2022 12:47:28 +0000 Subject: [PATCH 16/37] Update commonstatsmethods.md --- docs/part3/commonstatsmethods.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/part3/commonstatsmethods.md b/docs/part3/commonstatsmethods.md index e203b6098c2..a61f2275a40 100644 --- a/docs/part3/commonstatsmethods.md +++ b/docs/part3/commonstatsmethods.md @@ -297,7 +297,9 @@ For more heavy methods (eg the `MarkovChainMC`) you'll probably want to split th The `MarkovChainMC` method allows the user to produce the posterior pdf as a function of (in principle) any number of parameter of interest. In order to do so, you first need to create a workspace with more than one parameter, as explained in the [physics models](http://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part2/physicsmodels/) section. -For example, lets use the toy datacard [test/multiDim/toy-hgg-125.txt](https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/81x-root606/test/multiDim/toy-hgg-125.txt) (counting experiment which vaguely resembles the H→γγ analysis at 125 GeV) and convert the datacard into a workspace with 2 parameters, ggH and qqH cross sections using `text2workspace` with the option `-P HiggsAnalysis.CombinedLimit.PhysicsModel:floatingXSHiggs --PO modes=ggH,qqH`. +For example, lets use the toy datacard [test/multiDim/toy-hgg-125.txt](https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/81x-root606/test/multiDim/toy-hgg-125.txt) (counting experiment which vaguely resembles the H→γγ analysis at 125 GeV) and convert the datacard into a workspace with 2 parameters, ggH and qqH cross sections using `text2workspace`. + + text2workspace.py test/multiDim/toy-hgg-125.txt -P HiggsAnalysis.CombinedLimit.PhysicsModel:floatingXSHiggs --PO modes=ggH,qqH -o workspace.root Now we just run one (or more) MCMC chain(s) and save them in the output tree.By default, the nuisance parameters will be marginalized (integrated) over their pdfs. You can ignore the complaints about not being able to compute an upper limit (since for more than 1D, this isn't well defined), From 42c3c821f6b8a256234fd17ecb088d27bf0f4f65 Mon Sep 17 00:00:00 2001 From: nckw Date: Mon, 28 Nov 2022 12:48:09 +0000 Subject: [PATCH 17/37] Update commonstatsmethods.md --- docs/part3/commonstatsmethods.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/part3/commonstatsmethods.md b/docs/part3/commonstatsmethods.md index a61f2275a40..7f579489e59 100644 --- a/docs/part3/commonstatsmethods.md +++ b/docs/part3/commonstatsmethods.md @@ -303,7 +303,7 @@ For example, lets use the toy datacard [test/multiDim/toy-hgg-125.txt](https://g Now we just run one (or more) MCMC chain(s) and save them in the output tree.By default, the nuisance parameters will be marginalized (integrated) over their pdfs. You can ignore the complaints about not being able to compute an upper limit (since for more than 1D, this isn't well defined), - combine -M MarkovChainMC workspace.root --tries 1 --saveChain -i 1000000 -m 125 -s seed --noDefaultPrior=0 + combine -M MarkovChainMC workspace.root --tries 1 --saveChain -i 1000000 -m 125 -s 12345 --noDefaultPrior=0 The output of the markov chain is again a RooDataSet of weighted events distributed according to the posterior pdf (after you cut out the burn in part), so it can be used to make histograms or other distributions of the posterior pdf. See as an example [bayesPosterior2D.cxx](https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/81x-root606/test/multiDim/bayesPosterior2D.cxx). From 1889f090cc5b475f51e3cf37592a3d229a33bb70 Mon Sep 17 00:00:00 2001 From: nckw Date: Wed, 7 Dec 2022 11:47:27 +0000 Subject: [PATCH 18/37] Update index.md --- docs/index.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/index.md b/docs/index.md index a6b11bb196f..b4d7b8f8a45 100644 --- a/docs/index.md +++ b/docs/index.md @@ -66,9 +66,9 @@ scramv1 b clean; scramv1 b # always make a clean build ## Standalone version -The standalone version can be easily compiled using the \verb@cvmfs@ as it relies on dependencies which are already installed at [/cvmfs/cms.cern.ch/](/cvmfs/cms.cern.ch/). +The standalone version can be easily compiled using the `cvmfs` as it relies on dependencies which are already installed at [/cvmfs/cms.cern.ch/](/cvmfs/cms.cern.ch/). -Access to `/cvmfs/cms.cern.ch/` can be obtained from lxplus machines or via `CernVM`, by adding the `CMS` group to the CVMFS Configuration. A minimal `CernVM` working context setup can be found in the CernVM Marketplace under `Experimental/HiggsCombine` or at [https://cernvm-online.cern.ch/context/view/9ee5960ce4b143f5829e72bbbb26d382](https://cernvm-online.cern.ch/context/view/9ee5960ce4b143f5829e72bbbb26d382). At least 2GB of disk space should be reserved on the virtual machine for Combine to work properly. In case you do not want to use the `cvmfs` area, you will need to adapt the location of the dependencies listed in both the `Makefile` and `env_standalone.sh` files. +Access to `/cvmfs/cms.cern.ch/` can be obtained from lxplus machines or via `CernVM`, by adding the `CMS` group to the CVMFS Configuration. See instructions for [standalone version](http://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/#standalone-version-of-combine). In case you do not want to use the `cvmfs` area, you will need to adapt the location of the dependencies listed in both the `Makefile` and `env_standalone.sh` files. ``` git clone https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit.git HiggsAnalysis/CombinedLimit From 56c447151e4c978362444b62ee8b3e403f5088df Mon Sep 17 00:00:00 2001 From: nckw Date: Wed, 7 Dec 2022 15:22:39 +0000 Subject: [PATCH 19/37] Update runningthetool.md --- docs/part3/runningthetool.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/part3/runningthetool.md b/docs/part3/runningthetool.md index 96ffe2cb5d3..8e18e7cf774 100644 --- a/docs/part3/runningthetool.md +++ b/docs/part3/runningthetool.md @@ -15,7 +15,7 @@ The option `-M` allows to chose the method used. There are several groups of sta - `BayesianSimple`: performing a classical numerical integration (for simple models only) - `MarkovChainMC`: performing Markov Chain integration, for arbitrarily complex models. - **Frequentist** or hybrid bayesian-frequentist methods: - - `HybridNew`: compute modified frequentist limits according to several possible prescriptions + - `HybridNew`: compute modified frequentist limits or significance/p-values according to several possible prescriptions with toys. - **Fitting** - `FitDiagnostics`: performs maximum likelihood fits to extract the signal yield and provide diagnostic tools such as pre and post-fit models and correlations - `MultiDimFit`: perform maximum likelihood fits in multiple parameters and likelihood scans From 1491724f8cfbc7f27da983980e34d8dcc40127e7 Mon Sep 17 00:00:00 2001 From: nckw Date: Wed, 7 Dec 2022 16:05:10 +0000 Subject: [PATCH 20/37] Update runningthetool.md --- docs/part3/runningthetool.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/part3/runningthetool.md b/docs/part3/runningthetool.md index 8e18e7cf774..c2880bf4664 100644 --- a/docs/part3/runningthetool.md +++ b/docs/part3/runningthetool.md @@ -15,7 +15,7 @@ The option `-M` allows to chose the method used. There are several groups of sta - `BayesianSimple`: performing a classical numerical integration (for simple models only) - `MarkovChainMC`: performing Markov Chain integration, for arbitrarily complex models. - **Frequentist** or hybrid bayesian-frequentist methods: - - `HybridNew`: compute modified frequentist limits or significance/p-values according to several possible prescriptions with toys. + - `HybridNew`: compute modified frequentist limits, significance/p-values and confidence intervals according to several possible prescriptions with toys. - **Fitting** - `FitDiagnostics`: performs maximum likelihood fits to extract the signal yield and provide diagnostic tools such as pre and post-fit models and correlations - `MultiDimFit`: perform maximum likelihood fits in multiple parameters and likelihood scans From ce5c9d2b88d915f1566e46d1258ba591d07e86f4 Mon Sep 17 00:00:00 2001 From: nckw Date: Tue, 13 Dec 2022 11:13:33 +0000 Subject: [PATCH 21/37] Update runningthetool.md --- docs/part3/runningthetool.md | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/part3/runningthetool.md b/docs/part3/runningthetool.md index c2880bf4664..0c168db872b 100644 --- a/docs/part3/runningthetool.md +++ b/docs/part3/runningthetool.md @@ -145,11 +145,15 @@ The option `-t` is used to specify to combine to first generate a toy dataset(s) * `-t -1` will produce an Asimov dataset in which statistical fluctuations are suppressed. The procedure to generate this Asimov dataset depends on which type of analysis you are using, see below for details. -The output file will contain the toys (as `RooDataSets` for the observables, including global observables) in the **toys** directory if the option `--saveToys` is provided. If you include this option, the `limit` TTree in the output will have an entry corresponding to the state of the POI used for the generation of the toy, with the value of **`quantileExpected`** set to **-2**. You can add additional branches using the `--trackParameters` option as described in the [common command line options](#common-command-line-options) section above. - !!! warning The default values of the nuisance parameters (or any parameter) are used to generate the toy. This means that if, for example, you are using parametric shapes and the parameters inside the workspace are set to arbitrary values, *those* arbitrary values will be used to generate the toy. This behaviour can be modified through the use of the option `--setParameters x=value_x,y=value_y...` which will set the values of the parameters (`x` and `y`) before toy generation. You can also load a snap-shot from a previous fit to set the nuisances to their *post-fit* values (see below). +The output file will contain the toys (as `RooDataSets` for the observables, including global observables) in the **toys** directory if the option `--saveToys` is provided. If you include this option, the `limit` TTree in the output will have an entry corresponding to the state of the POI used for the generation of the toy, with the value of **`quantileExpected`** set to **-2**. + +!!! info + The branches that are created by methods like `MultiDimFit` *will not* show the values used to generate the toy. If you also want the TTree to show the values of the POIs used to generate to toy, you should add additional branches using the `--trackParameters` option as described in the [common command line options](#common-command-line-options) section above. These branches will behave as expected when adding the option `--saveToys`. + + #### Asimov datasets If you are using wither `-t -1` or using `AsymptoticLimits`, combine will calculate results based on an Asimov dataset. From b34142c98b6c12b597a7afa5d95966cf4cd2b701 Mon Sep 17 00:00:00 2001 From: nckw Date: Thu, 22 Dec 2022 10:46:23 +0000 Subject: [PATCH 22/37] Update commonstatsmethods.md --- docs/part3/commonstatsmethods.md | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/docs/part3/commonstatsmethods.md b/docs/part3/commonstatsmethods.md index 7f579489e59..0b0faa22011 100644 --- a/docs/part3/commonstatsmethods.md +++ b/docs/part3/commonstatsmethods.md @@ -184,10 +184,11 @@ Done in 0.14 min (cpu), 0.15 min (real) Again, the resulting limit tree will contain the result. You can also save the chains using the option `--saveChain` which will then also be included in the output file. -Exclusion regions can be made from the posterior once an ordering principle is defined to decide how to grow the contour (there's infinite possible regions that contain 68% of the posterior pdf...) -Below is a simple example script which can be used to plot the posterior distribution from these chains and calculate the *smallest* such region. Note that in this example we are ignoring the burn-in (but you can add it by just editing `for i in range(mychain.numEntries()):` to `for i in range(200,mychain.numEntries()):` eg for a burn-in of 200. +Exclusion regions can be made from the posterior once an ordering principle is defined to decide how to grow the contour (there's infinite possible regions that contain 68% of the posterior pdf). Below is a simple example script which can be used to plot the posterior distribution from these chains and calculate the *smallest* such region. Note that in this example we are ignoring the burn-in (but you can add it by just editing `for i in range(mychain.numEntries()):` to `for i in range(200,mychain.numEntries()):` eg for a burn-in of 200. -```python +
+Show example script +

 import ROOT
 
 rmin = 0
@@ -248,7 +249,8 @@ ll.Draw()
 lu.Draw()
 
 print " %g %% (%g %%) interval (target)  = %g < r < %g "%(trueCL,CL,vl,vu)
-```
+
+
Running the script on the output file produced for the same datacard (including the `--saveChain` option) will produce the following output @@ -317,7 +319,6 @@ bayesPosterior2D("bayes2D","Posterior PDF") ![](images/bayes2D.png) - ## Computing Limits with toys The `HybridNew` method is used to compute either the hybrid bayesian-frequentist limits popularly known as "CLs of LEP or Tevatron type" or the fully frequentist limits which are the current recommended method by the LHC Higgs Combination Group. Note that these methods can be resource intensive for complex models. From abb232a67eb00126a54e18bc9f5a98ea9b6f959c Mon Sep 17 00:00:00 2001 From: Dominic Stafford Date: Thu, 5 Jan 2023 17:39:56 +0100 Subject: [PATCH 23/37] Allow ChannelCompatibilityCheck to skip channels with no signal templates --- src/ChannelCompatibilityCheck.cc | 4 +++- src/FitterAlgoBase.cc | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/ChannelCompatibilityCheck.cc b/src/ChannelCompatibilityCheck.cc index 35dc35bea7c..f8cdcfd36cc 100644 --- a/src/ChannelCompatibilityCheck.cc +++ b/src/ChannelCompatibilityCheck.cc @@ -113,7 +113,9 @@ bool ChannelCompatibilityCheck::runSpecific(RooWorkspace *w, RooStats::ModelConf } for (std::map::const_iterator it = rs.begin(), ed = rs.end(); it != ed; ++it) { RooRealVar *ri = (RooRealVar*) result_freeform->floatParsFinal().find(it->second.c_str()); - if (runMinos_ && do95_) { + if (ri == NULL){ + printf("Parameter %s not found in channel %s. Does this region contain signal templates?\n", r->GetName(), it->first.c_str()); + } else if (runMinos_ && do95_) { printf("Alternate fit: %s = %7.4f %+6.4f/%+6.4f (68%% CL) in channel %s\n", r->GetName(), ri->getVal(), ri->getAsymErrorLo(), ri->getAsymErrorHi(), it->first.c_str()); printf(" %s = %7.4f %+6.4f/%+6.4f (95%% CL) in channel %s\n", r->GetName(), ri->getVal(), ri->getMin("err95")-ri->getVal(), ri->getMax("err95")-ri->getVal(), it->first.c_str()); } else if (runMinos_) { diff --git a/src/FitterAlgoBase.cc b/src/FitterAlgoBase.cc index bd0617abe75..8c0d95147f2 100644 --- a/src/FitterAlgoBase.cc +++ b/src/FitterAlgoBase.cc @@ -317,6 +317,7 @@ RooFitResult *FitterAlgoBase::doFit(RooAbsPdf &pdf, RooAbsData &data, const RooA rfloat = ret->constPars().find(r.GetName()); if (!rfloat) { fprintf(sentry.trueStdOut(), "Skipping %s. Parameter not found in the RooFitResult.\n",r.GetName()); + continue ; } } //rfloat->Print("V"); From 09637694d6bf6515b2b699b7aae2c934d575acb3 Mon Sep 17 00:00:00 2001 From: Aashwin Basnet Date: Fri, 3 Mar 2023 17:36:02 +0100 Subject: [PATCH 24/37] random starting point implementation for 2d grid scan --- src/MultiDimFit.cc | 298 ++++++++++++++++++++++++++++++--------------- 1 file changed, 201 insertions(+), 97 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 14595ef538e..ad57c069e40 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -712,7 +712,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) wc_vals_vec_of_vec.push_back(default_start_pt_vec); } - // Append the random points to the vecotr of points to try + // 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 (randPointsRanges_ != "") { @@ -815,9 +815,15 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } // 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) { + srand(randPointsSeed_); + } + // get number of points per axis unsigned int nX, nY; if (pointsPerPoi.size() == 0) { @@ -871,106 +877,204 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) poiVals_[1] = y; poiVars_[0]->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; + + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////// Loop over rand points for each profiled POI to get best nll ///////////////////// + + bool ok; + int n_prof_params = specifiedVars_.size(); + + //Get vector of points to try + std::vector> wc_vals_vec_of_vec = {}; + + //Append the default start pt to the list of points to try + 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); } - // 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); + + //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 (randPointsRanges_ != "") { + rand_ranges_dict = getRangesDictFromInString(randPointsRanges_); } - 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); - } - } + for (int pt_idx=0; pt_idx wc_vals_vec; + for (int prof_param_idx=0; prof_param_idxGetName()][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)); + } + // 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 (verbose > 1){ + std::cout << "List of points to try:" << std::endl; + for (auto vals_vec: wc_vals_vec_of_vec) { + std::cout << "\tThe vals at this point: " << std::endl; + for (auto val: vals_vec) { + std::cout << "\t\tPoint val: " << val << std::endl; + } + } } - } - } + // Loop over starting points to try for the prof WCs + float current_best_nll = 0; // This variable will be properly initialized within the for loop, just initialize to a dummy value here to avoid compiling warnings + for (unsigned int start_pt_idx=0; start_pt_idxsetVal(x); + poiVars_[1]->setVal(y); + + //Loop over prof POIs and set their values + if (verbose > 1) std::cout<<"\n\tStart pt idx: "<< start_pt_idx << std::endl; + for (unsigned int var_idx=0; var_idx 1) std::cout<<"\t\tThe var name: "<GetName() << std::endl; + if(strcmp(specifiedVars_[var_idx]->GetName(),"r")==0){ + //Don't bother setting r to anything + if (verbose > 1) std::cout<<"\t\t\tSkipping var "<< specifiedVars_[var_idx]->GetName() << std::endl; + continue; + } + if (verbose > 1) std::cout << "\t\t\tRange before: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; + if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " before setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; + if (pointsRandProf_ > 0) { + specifiedVars_[var_idx]->removeRange(); // Do not impose a range if we're trying multiple random start points + } + specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(start_pt_idx).at(var_idx)); + if (verbose > 1) std::cout << "\t\t\tRange after: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; + if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " after setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; + } + //now we minimize + nll.clearEvalErrorLog(); + nll.getVal(); + deltaNLL_ = nll.getVal() - nll0; + 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; + } + bool skipme = hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_; + 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 (start_pt_idx==0) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } else if (nll.getVal() < current_best_nll) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } + } + 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); + if (start_pt_idx==0) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } else if (nll.getVal() < current_best_nll) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } + + } + } + } + } // End loop over random start points + } // End of the loop over y dim scan points + } // End of the loop over x dim scan points } else { // Use utils routine if n > 2 RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); CloseCoutSentry sentry(verbose < 2); From 58bbff0d8e363bf6243ae617dcb4c2f92537b8df Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 3 Dec 2021 16:52:18 -0500 Subject: [PATCH 25/37] Add ability to try rand start pts for prof POIs --- src/MultiDimFit.cc | 150 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 132 insertions(+), 18 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 25c52ca4741..6942d9a3f04 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -186,6 +186,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 " < 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; @@ -627,8 +631,12 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } if (n == 1) { + if (verbose > 1) std::cout << "The nll0 from initial fit: " << nll0 << std::endl; unsigned int points = pointsPerPoi.size() == 0 ? points_ : pointsPerPoi[0]; + // Set seed for random points + srand(1); + double xspacing = (pmax[0]-pmin[0]) / points; double xspacingOffset = 0.5; if (alignEdges_) { @@ -666,42 +674,148 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } //if (verbose > 1) std::cout << "Point " << i << "/" << points << " " << poiVars_[0]->GetName() << " = " << x << std::endl; - std::cout << "Point " << i << "/" << points << " " << poiVars_[0]->GetName() << " = " << x << std::endl; + std::cout << "Point " << i << "/" << points << " " << poiVars_[0]->GetName() << " = " << x << std::endl; + *params = snap; + poiVals_[0] = x; + poiVars_[0]->setVal(x); + + + /////////////////////////////////////////////////////////////////////////////////////////////////// + /////////////////// Loop over rand points for each profiled POI to get best nll /////////////////// + + bool ok; + int n_rand_start_pts_to_try = 0; + int n_prof_params = specifiedVars_.size(); + std::vector nll_at_alt_start_pts; + + // Get vector of points to try + float prof_start_pt_range_max = 20.0; // Is this what we want? + std::vector> wc_vals_vec_of_vec = {}; + + // Append the defualt start pt to the list of points to try + 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 vecotr of points to try + for (int pt_idx=0; pt_idx wc_vals_vec; + for (int prof_param_idx=0; prof_param_idx 1) { + std::cout << "List of points to try:" << std::endl; + for (auto vals_vec: wc_vals_vec_of_vec) { + std::cout << "\tThe vals at this point: " << std::endl; + for (auto val: vals_vec) { + std::cout << "\t\tPoint val: " << val << std::endl; + } + } + } + + // Loop over starting points to try for the prof WCs + for (unsigned int start_pt_idx=0; start_pt_idxsetVal(x); + + // Loop over prof POIs and set their values + if (verbose > 1) std::cout << "\n\tStart pt idx: " << start_pt_idx << std::endl; + for (unsigned int var_idx=0; var_idx 1) std::cout << "\t\tThe var name: " << specifiedVars_[var_idx]->GetName() << std::endl; + if (strcmp(specifiedVars_[var_idx]->GetName(),"r")==0) { + // Don't bother setting r to anything + if (verbose > 1) std::cout << "\t\t\tSkipping var " << specifiedVars_[var_idx]->GetName() << std::endl; + continue; + } + if (verbose > 1) std::cout << "\t\t\tRange before: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; + if (verbose > 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]->removeRange(); // Do not impose a range + specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(start_pt_idx).at(var_idx)); + if (verbose > 1) std::cout << "\t\t\tRange after: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; + if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " after setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; + } + + // 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; + } + ok = fastScan_ || (hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_) || utils::countFloating(*params)==0 ? + true : + minim.minimize(verbose-1); + + if (verbose > 1) std::cout << "\t\tThe nll.getVal() is: " << nll.getVal() << std::endl; + nll_at_alt_start_pts.push_back(nll.getVal()); + + } + + // Rerun the fit for the point with the min nll + int min_nll_idx; + min_nll_idx = std::min_element(nll_at_alt_start_pts.begin(),nll_at_alt_start_pts.end()) - nll_at_alt_start_pts.begin(); + if (verbose > 1) std::cout << "\tMin nll at idx: " << min_nll_idx << " " << nll_at_alt_start_pts.at(min_nll_idx) << std::endl; + if (verbose > 1) std::cout << "\tResetting for rerunning at the point that gives best nll:" << std::endl; *params = snap; poiVals_[0] = x; poiVars_[0]->setVal(x); + for (unsigned int var_idx=0; var_idx 1) std::cout << "\t\tSetting " << specifiedVars_[var_idx]->GetName() << " to " << wc_vals_vec_of_vec.at(min_nll_idx).at(var_idx) << std::endl; + specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(min_nll_idx).at(var_idx)); + } // 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(); - } + 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 ? + 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(); - } + for(unsigned int j=0; jgetVal(); + } + for(unsigned int j=0; jgetVal(); + } + for(unsigned int j=0; jgetIndex(); + } Combine::commitPoint(true, /*quantile=*/prob); } + + std::cout << "" << std::endl; } } else if (n == 2) { RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); From 0af5422002beb57fe47ee1fa8d7bbf73210c47a4 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 3 Dec 2021 17:52:37 -0500 Subject: [PATCH 26/37] Option for specifying number of random start points --- interface/MultiDimFit.h | 2 ++ src/MultiDimFit.cc | 12 +++++++----- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index cf9e77a8353..6253e01b339 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -65,6 +65,8 @@ class MultiDimFit : public FitterAlgoBase { static std::string robustHesseLoad_; static std::string robustHesseSave_; + static int pointsRandProf_; + static std::string saveSpecifiedFuncs_; static std::string saveSpecifiedNuis_; static std::string saveSpecifiedIndex_; diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 6942d9a3f04..b3c7e22fbe2 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -60,6 +60,7 @@ float MultiDimFit::centeredRange_ = -1.0; bool MultiDimFit::robustHesse_ = false; std::string MultiDimFit::robustHesseLoad_ = ""; std::string MultiDimFit::robustHesseSave_ = ""; +int MultiDimFit::pointsRandProf_ = 0; std::string MultiDimFit::saveSpecifiedFuncs_; @@ -110,6 +111,7 @@ MultiDimFit::MultiDimFit() : ("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") ; } @@ -631,7 +633,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } if (n == 1) { - if (verbose > 1) std::cout << "The nll0 from initial fit: " << nll0 << std::endl; + 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 @@ -684,7 +686,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) /////////////////// Loop over rand points for each profiled POI to get best nll /////////////////// bool ok; - int n_rand_start_pts_to_try = 0; int n_prof_params = specifiedVars_.size(); std::vector nll_at_alt_start_pts; @@ -700,7 +701,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) wc_vals_vec_of_vec.push_back(default_start_pt_vec); // Append the random points to the vecotr of points to try - for (int pt_idx=0; pt_idx wc_vals_vec; for (int prof_param_idx=0; prof_param_idx 1) std::cout << "\t\t\tRange before: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; if (verbose > 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]->removeRange(); // Do not impose a range + if (pointsRandProf_ > 0) { + specifiedVars_[var_idx]->removeRange(); // Do not impose a range if we're trying multiple random start points + } specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(start_pt_idx).at(var_idx)); if (verbose > 1) std::cout << "\t\t\tRange after: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " after setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; @@ -815,7 +818,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) Combine::commitPoint(true, /*quantile=*/prob); } - std::cout << "" << std::endl; } } else if (n == 2) { RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); From 09ad136fa96b431b61ac92f4c401bf7001be211e Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 3 Dec 2021 18:11:54 -0500 Subject: [PATCH 27/37] If using pointsRandProf_, turn on saveInactivePOI_ --- src/MultiDimFit.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index b3c7e22fbe2..164a2655237 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -141,6 +141,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); From 3360d4c41e0f86b94633adce6d6bb5e6dc7d08d1 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Tue, 5 Jul 2022 00:11:01 -0400 Subject: [PATCH 28/37] Do not save nll if fit fails, also add ability to pass ranges for rand start points --- interface/MultiDimFit.h | 4 ++++ src/MultiDimFit.cc | 39 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index 6253e01b339..43475d1c2c9 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -66,6 +66,7 @@ class MultiDimFit : public FitterAlgoBase { static std::string robustHesseSave_; static int pointsRandProf_; + static std::string randPointsRanges_; static std::string saveSpecifiedFuncs_; static std::string saveSpecifiedNuis_; @@ -95,6 +96,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/src/MultiDimFit.cc b/src/MultiDimFit.cc index 164a2655237..1b846bf90d2 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -61,6 +61,7 @@ bool MultiDimFit::robustHesse_ = false; std::string MultiDimFit::robustHesseLoad_ = ""; std::string MultiDimFit::robustHesseSave_ = ""; int MultiDimFit::pointsRandProf_ = 0; +std::string MultiDimFit::randPointsRanges_; std::string MultiDimFit::saveSpecifiedFuncs_; @@ -112,6 +113,7 @@ MultiDimFit::MultiDimFit() : ("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") + ("randPointsRanges", boost::program_options::value(&randPointsRanges_)->default_value(""), "Range from which to draw random start points for the profiled POIs") ; } @@ -695,7 +697,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) std::vector nll_at_alt_start_pts; // Get vector of points to try - float prof_start_pt_range_max = 20.0; // Is this what we want? std::vector> wc_vals_vec_of_vec = {}; // Append the defualt start pt to the list of points to try @@ -706,9 +707,19 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) wc_vals_vec_of_vec.push_back(default_start_pt_vec); // Append the random points to the vecotr 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 (randPointsRanges_ != "") { + rand_ranges_dict = getRangesDictFromInString(randPointsRanges_); + } for (int pt_idx=0; pt_idx wc_vals_vec; for (int prof_param_idx=0; prof_param_idxGetName()][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)); + } // 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); @@ -771,7 +782,11 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) minim.minimize(verbose-1); if (verbose > 1) std::cout << "\t\tThe nll.getVal() is: " << nll.getVal() << std::endl; - nll_at_alt_start_pts.push_back(nll.getVal()); + if (ok) { + nll_at_alt_start_pts.push_back(nll.getVal()); + } else { + nll_at_alt_start_pts.push_back(999999.9); // Placeholder for now + } } @@ -1355,3 +1370,23 @@ 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; +} From 709c63757a098c64eb25d0ff474b63427189fffe Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Tue, 5 Jul 2022 01:08:53 -0400 Subject: [PATCH 29/37] Add option for skipping default start point --- interface/MultiDimFit.h | 1 + src/MultiDimFit.cc | 12 ++++++++---- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index 43475d1c2c9..65dded8e250 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -84,6 +84,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) ; diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 1b846bf90d2..2779206cad8 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -81,6 +81,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") @@ -104,6 +105,7 @@ 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") @@ -700,11 +702,13 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) std::vector> wc_vals_vec_of_vec = {}; // Append the defualt start pt to the list of points to try - std::vector default_start_pt_vec; - for (int prof_param_idx=0; prof_param_idxgetVal()); + 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); } - wc_vals_vec_of_vec.push_back(default_start_pt_vec); // Append the random points to the vecotr of points to try float prof_start_pt_range_max = 20.0; // Default to 20 if we're not asking for custom ranges From 8e7e17d1072a6a8539c810d1c8e0bd58d38a14a3 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Tue, 5 Jul 2022 23:07:08 -0400 Subject: [PATCH 30/37] Save each point if best nll so far to avoid rerunning fit at end --- src/MultiDimFit.cc | 81 ++++++++++++++-------------------------------- 1 file changed, 25 insertions(+), 56 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 2779206cad8..000492751b9 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -743,6 +743,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } // Loop over starting points to try for the prof WCs + float current_best_nll; for (unsigned int start_pt_idx=0; start_pt_idx maxDeltaNLLForProf_) || utils::countFloating(*params)==0 ? - true : + ok = fastScan_ || (hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_) || utils::countFloating(*params)==0 ? + true : minim.minimize(verbose-1); - if (verbose > 1) std::cout << "\t\tThe nll.getVal() is: " << nll.getVal() << std::endl; if (ok) { - nll_at_alt_start_pts.push_back(nll.getVal()); - } else { - nll_at_alt_start_pts.push_back(999999.9); // Placeholder for now - } - - } - - // Rerun the fit for the point with the min nll - int min_nll_idx; - min_nll_idx = std::min_element(nll_at_alt_start_pts.begin(),nll_at_alt_start_pts.end()) - nll_at_alt_start_pts.begin(); - if (verbose > 1) std::cout << "\tMin nll at idx: " << min_nll_idx << " " << nll_at_alt_start_pts.at(min_nll_idx) << std::endl; - if (verbose > 1) std::cout << "\tResetting for rerunning at the point that gives best nll:" << std::endl; - *params = snap; - poiVals_[0] = x; - poiVars_[0]->setVal(x); - for (unsigned int var_idx=0; var_idx 1) std::cout << "\t\tSetting " << specifiedVars_[var_idx]->GetName() << " to " << wc_vals_vec_of_vec.at(min_nll_idx).at(var_idx) << std::endl; - specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(min_nll_idx).at(var_idx)); - } - // 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(); + 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 (start_pt_idx==0) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } else if (nll.getVal() < current_best_nll) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } } - Combine::commitPoint(true, /*quantile=*/0); - continue; - } - 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); - } + } // End loop over random start points - } + } // End of the loop over scan points } else if (n == 2) { RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); CloseCoutSentry sentry(verbose < 2); From e2fda6ecbfd4e25e4af5c396f072be010bebaa0b Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Wed, 6 Jul 2022 16:15:34 -0400 Subject: [PATCH 31/37] Pass rand seed for rand start points --- interface/MultiDimFit.h | 1 + src/MultiDimFit.cc | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index 65dded8e250..eb398032bc0 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -67,6 +67,7 @@ class MultiDimFit : public FitterAlgoBase { static int pointsRandProf_; static std::string randPointsRanges_; + static int randPointsSeed_; static std::string saveSpecifiedFuncs_; static std::string saveSpecifiedNuis_; diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 000492751b9..0bfc1ffea30 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -61,6 +61,7 @@ bool MultiDimFit::robustHesse_ = false; std::string MultiDimFit::robustHesseLoad_ = ""; std::string MultiDimFit::robustHesseSave_ = ""; int MultiDimFit::pointsRandProf_ = 0; +int MultiDimFit::randPointsSeed_ = 0; std::string MultiDimFit::randPointsRanges_; @@ -115,6 +116,7 @@ MultiDimFit::MultiDimFit() : ("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") ("randPointsRanges", boost::program_options::value(&randPointsRanges_)->default_value(""), "Range from which to draw random start points for the profiled POIs") ; } @@ -646,7 +648,9 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) unsigned int points = pointsPerPoi.size() == 0 ? points_ : pointsPerPoi[0]; // Set seed for random points - srand(1); + if (pointsRandProf_ > 0) { + srand(randPointsSeed_); + } double xspacing = (pmax[0]-pmin[0]) / points; double xspacingOffset = 0.5; From c0fcac95564d4d191f4496169655f2a6339d985c Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 8 Jul 2022 15:08:54 -0400 Subject: [PATCH 32/37] Small cleanups --- src/MultiDimFit.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 0bfc1ffea30..334b6f9a5ff 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -700,7 +700,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) bool ok; int n_prof_params = specifiedVars_.size(); - std::vector nll_at_alt_start_pts; // Get vector of points to try std::vector> wc_vals_vec_of_vec = {}; @@ -747,7 +746,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } // Loop over starting points to try for the prof WCs - float current_best_nll; + float current_best_nll = 0; // This variable will be properly initialized within the for loop, just initialize to a dummy value here to avoid compiling warnings for (unsigned int start_pt_idx=0; start_pt_idx Date: Fri, 3 Mar 2023 17:36:02 +0100 Subject: [PATCH 33/37] random starting point implementation for 2d grid scan --- src/MultiDimFit.cc | 298 ++++++++++++++++++++++++++++++--------------- 1 file changed, 201 insertions(+), 97 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 334b6f9a5ff..7c24f4ef39e 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -713,7 +713,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) wc_vals_vec_of_vec.push_back(default_start_pt_vec); } - // Append the random points to the vecotr of points to try + // 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 (randPointsRanges_ != "") { @@ -816,9 +816,15 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) } // 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) { + srand(randPointsSeed_); + } + // get number of points per axis unsigned int nX, nY; if (pointsPerPoi.size() == 0) { @@ -872,106 +878,204 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) poiVals_[1] = y; poiVars_[0]->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; + + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////// Loop over rand points for each profiled POI to get best nll ///////////////////// + + bool ok; + int n_prof_params = specifiedVars_.size(); + + //Get vector of points to try + std::vector> wc_vals_vec_of_vec = {}; + + //Append the default start pt to the list of points to try + 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); } - // 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); + + //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 (randPointsRanges_ != "") { + rand_ranges_dict = getRangesDictFromInString(randPointsRanges_); } - 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); - } - } + for (int pt_idx=0; pt_idx wc_vals_vec; + for (int prof_param_idx=0; prof_param_idxGetName()][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)); + } + // 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 (verbose > 1){ + std::cout << "List of points to try:" << std::endl; + for (auto vals_vec: wc_vals_vec_of_vec) { + std::cout << "\tThe vals at this point: " << std::endl; + for (auto val: vals_vec) { + std::cout << "\t\tPoint val: " << val << std::endl; + } + } } - } - } + // Loop over starting points to try for the prof WCs + float current_best_nll = 0; // This variable will be properly initialized within the for loop, just initialize to a dummy value here to avoid compiling warnings + for (unsigned int start_pt_idx=0; start_pt_idxsetVal(x); + poiVars_[1]->setVal(y); + + //Loop over prof POIs and set their values + if (verbose > 1) std::cout<<"\n\tStart pt idx: "<< start_pt_idx << std::endl; + for (unsigned int var_idx=0; var_idx 1) std::cout<<"\t\tThe var name: "<GetName() << std::endl; + if(strcmp(specifiedVars_[var_idx]->GetName(),"r")==0){ + //Don't bother setting r to anything + if (verbose > 1) std::cout<<"\t\t\tSkipping var "<< specifiedVars_[var_idx]->GetName() << std::endl; + continue; + } + if (verbose > 1) std::cout << "\t\t\tRange before: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; + if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " before setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; + if (pointsRandProf_ > 0) { + specifiedVars_[var_idx]->removeRange(); // Do not impose a range if we're trying multiple random start points + } + specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(start_pt_idx).at(var_idx)); + if (verbose > 1) std::cout << "\t\t\tRange after: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; + if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " after setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; + } + //now we minimize + nll.clearEvalErrorLog(); + nll.getVal(); + deltaNLL_ = nll.getVal() - nll0; + 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; + } + bool skipme = hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_; + 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 (start_pt_idx==0) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } else if (nll.getVal() < current_best_nll) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } + } + 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); + if (start_pt_idx==0) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } else if (nll.getVal() < current_best_nll) { + Combine::commitPoint(true, /*quantile=*/prob); + current_best_nll = nll.getVal(); + } + + } + } + } + } // End loop over random start points + } // End of the loop over y dim scan points + } // End of the loop over x dim scan points } else { // Use utils routine if n > 2 RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); CloseCoutSentry sentry(verbose < 2); From 1f7b676ac7c376ee05245d659a0608a51e5ba094 Mon Sep 17 00:00:00 2001 From: Aashwin Basnet Date: Mon, 10 Apr 2023 11:44:57 +0200 Subject: [PATCH 34/37] created a separate class for random starting point implementation --- interface/MultiDimFit.h | 2 + interface/RandStartPt.h | 45 ++++++ src/MultiDimFit.cc | 349 ++++------------------------------------ src/RandStartPt.cc | 260 ++++++++++++++++++++++++++++++ 4 files changed, 341 insertions(+), 315 deletions(-) create mode 100644 interface/RandStartPt.h create mode 100644 src/RandStartPt.cc diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index eb398032bc0..a7c1e0ab773 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -21,6 +21,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); diff --git a/interface/RandStartPt.h b/interface/RandStartPt.h new file mode 100644 index 00000000000..09db515680f --- /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 randptsranges_; + 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 randptsranges, 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::auto_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::auto_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 7c24f4ef39e..d6971518ae9 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -19,6 +19,7 @@ #include "HiggsAnalysis/CombinedLimit/interface/utils.h" #include "HiggsAnalysis/CombinedLimit/interface/RobustHesse.h" #include "HiggsAnalysis/CombinedLimit/interface/ProfilingTools.h" +#include "HiggsAnalysis/CombinedLimit/interface/RandStartPt.h" #include #include @@ -596,7 +597,6 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) unsigned int n = poi_.size(); //if (poi_.size() > 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()); @@ -649,8 +649,16 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) // 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; @@ -694,125 +703,13 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) poiVals_[0] = x; poiVars_[0]->setVal(x); + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + ////////////////// 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 pointsRandProf_ to 0 or leave it unspecified. ///////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////////// Loop over rand points for each profiled POI to get best nll /////////////////// - - bool ok; - int n_prof_params = specifiedVars_.size(); - - // Get vector of points to try - std::vector> wc_vals_vec_of_vec = {}; - - // Append the defualt start pt to the list of points to try - 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 (randPointsRanges_ != "") { - rand_ranges_dict = getRangesDictFromInString(randPointsRanges_); - } - for (int pt_idx=0; pt_idx wc_vals_vec; - for (int prof_param_idx=0; prof_param_idxGetName()][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)); - } - // 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 (verbose > 1) { - std::cout << "List of points to try:" << std::endl; - for (auto vals_vec: wc_vals_vec_of_vec) { - std::cout << "\tThe vals at this point: " << std::endl; - for (auto val: vals_vec) { - std::cout << "\t\tPoint val: " << val << std::endl; - } - } - } - - // Loop over starting points to try for the prof WCs - float current_best_nll = 0; // This variable will be properly initialized within the for loop, just initialize to a dummy value here to avoid compiling warnings - for (unsigned int start_pt_idx=0; start_pt_idxsetVal(x); - - // Loop over prof POIs and set their values - if (verbose > 1) std::cout << "\n\tStart pt idx: " << start_pt_idx << std::endl; - for (unsigned int var_idx=0; var_idx 1) std::cout << "\t\tThe var name: " << specifiedVars_[var_idx]->GetName() << std::endl; - if (strcmp(specifiedVars_[var_idx]->GetName(),"r")==0) { - // Don't bother setting r to anything - if (verbose > 1) std::cout << "\t\t\tSkipping var " << specifiedVars_[var_idx]->GetName() << std::endl; - continue; - } - if (verbose > 1) std::cout << "\t\t\tRange before: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; - if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " before setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; - if (pointsRandProf_ > 0) { - specifiedVars_[var_idx]->removeRange(); // Do not impose a range if we're trying multiple random start points - } - specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(start_pt_idx).at(var_idx)); - if (verbose > 1) std::cout << "\t\t\tRange after: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; - if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " after setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; - } - - // 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; - } - 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); - if (start_pt_idx==0) { - Combine::commitPoint(true, /*quantile=*/prob); - current_best_nll = nll.getVal(); - } else if (nll.getVal() < current_best_nll) { - Combine::commitPoint(true, /*quantile=*/prob); - current_best_nll = nll.getVal(); - } - } - - } // End loop over random start points + RandStartPt randStartPt(nll, specifiedVars_, specifiedVals_, skipDefaultStart_, randPointsRanges_, pointsRandProf_, 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) { @@ -821,9 +718,17 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) CloseCoutSentry sentry(verbose < 2); //Set seed for random points - if(pointsRandProf_ > 0) { + if (pointsRandProf_ > 0) { + std::cout<<"\n************************************************************************************"<setVal(x); poiVars_[1]->setVal(y); - ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //////////////////////////////////// Loop over rand points for each profiled POI to get best nll ///////////////////// - - bool ok; - int n_prof_params = specifiedVars_.size(); - - //Get vector of points to try - std::vector> wc_vals_vec_of_vec = {}; - - //Append the default start pt to the list of points to try - 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 (randPointsRanges_ != "") { - rand_ranges_dict = getRangesDictFromInString(randPointsRanges_); - } - for (int pt_idx=0; pt_idx wc_vals_vec; - for (int prof_param_idx=0; prof_param_idxGetName()][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)); - } - // 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); - } + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////// 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. ///////////////////////////////// - // Print vector of points to try - if (verbose > 1){ - std::cout << "List of points to try:" << std::endl; - for (auto vals_vec: wc_vals_vec_of_vec) { - std::cout << "\tThe vals at this point: " << std::endl; - for (auto val: vals_vec) { - std::cout << "\t\tPoint val: " << val << std::endl; - } - } - } + RandStartPt randStartPt(nll, specifiedVars_, specifiedVals_, skipDefaultStart_, randPointsRanges_, 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 - // Loop over starting points to try for the prof WCs - float current_best_nll = 0; // This variable will be properly initialized within the for loop, just initialize to a dummy value here to avoid compiling warnings - for (unsigned int start_pt_idx=0; start_pt_idxsetVal(x); - poiVars_[1]->setVal(y); - - //Loop over prof POIs and set their values - if (verbose > 1) std::cout<<"\n\tStart pt idx: "<< start_pt_idx << std::endl; - for (unsigned int var_idx=0; var_idx 1) std::cout<<"\t\tThe var name: "<GetName() << std::endl; - if(strcmp(specifiedVars_[var_idx]->GetName(),"r")==0){ - //Don't bother setting r to anything - if (verbose > 1) std::cout<<"\t\t\tSkipping var "<< specifiedVars_[var_idx]->GetName() << std::endl; - continue; - } - if (verbose > 1) std::cout << "\t\t\tRange before: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; - if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " before setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; - if (pointsRandProf_ > 0) { - specifiedVars_[var_idx]->removeRange(); // Do not impose a range if we're trying multiple random start points - } - specifiedVars_[var_idx]->setVal(wc_vals_vec_of_vec.at(start_pt_idx).at(var_idx)); - if (verbose > 1) std::cout << "\t\t\tRange after: " << specifiedVars_[var_idx]->getMin() << " " << specifiedVars_[var_idx]->getMax() << std::endl; - if (verbose > 1) std::cout << "\t\t\t" << specifiedVars_[var_idx]->GetName() << " after setting: " << specifiedVars_[var_idx]->getVal() << " += " << specifiedVars_[var_idx]->getError() << std::endl; - } - //now we minimize - nll.clearEvalErrorLog(); - nll.getVal(); - deltaNLL_ = nll.getVal() - nll0; - 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; - } - bool skipme = hasMaxDeltaNLLForProf_ && (nll.getVal() - nll0) > maxDeltaNLLForProf_; - 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 (start_pt_idx==0) { - Combine::commitPoint(true, /*quantile=*/prob); - current_best_nll = nll.getVal(); - } else if (nll.getVal() < current_best_nll) { - Combine::commitPoint(true, /*quantile=*/prob); - current_best_nll = nll.getVal(); - } - } - 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); - if (start_pt_idx==0) { - Combine::commitPoint(true, /*quantile=*/prob); - current_best_nll = nll.getVal(); - } else if (nll.getVal() < current_best_nll) { - Combine::commitPoint(true, /*quantile=*/prob); - current_best_nll = nll.getVal(); - } - - } - } - } - } // End loop over random start points - } // End of the loop over y dim scan points - } // End of the loop over x dim scan points } else { // Use utils routine if n > 2 RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); CloseCoutSentry sentry(verbose < 2); @@ -1470,3 +1188,4 @@ std::map> MultiDimFit::getRangesDictFromInString } return out_range_dict; } + diff --git a/src/RandStartPt.cc b/src/RandStartPt.cc new file mode 100644 index 00000000000..554c855ba19 --- /dev/null +++ b/src/RandStartPt.cc @@ -0,0 +1,260 @@ +#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 randptsranges, 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), + randptsranges_(randptsranges), + 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 (randptsranges_ != "") { + rand_ranges_dict = RandStartPt::getRangesDictFromInString(randptsranges_); + } + for (int pt_idx = 0; pt_idx wc_vals_vec; + for (int prof_param_idx=0; prof_param_idxGetName()][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)); + } + //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: "<> 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; + if (numrandpts_ > 0){ + specifiedvars_[var_idx]->removeRange(); //Do not impose a range if we are trying multiple random start points + } + 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::auto_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::auto_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); + } + } + } + } +} From 0e06eb359ec3610f23ab24802cce339d6633f2af Mon Sep 17 00:00:00 2001 From: Aashwin Basnet Date: Tue, 2 Jan 2024 16:49:54 +0100 Subject: [PATCH 35/37] Making changes to be able to compile with CMSSW_11_X plus some other code cleaning --- interface/RandStartPt.h | 4 ++-- src/MultiDimFit.cc | 49 ++++++++++++++++++++++------------------- src/RandStartPt.cc | 4 ++-- 3 files changed, 30 insertions(+), 27 deletions(-) diff --git a/interface/RandStartPt.h b/interface/RandStartPt.h index 09db515680f..8d41e823f20 100644 --- a/interface/RandStartPt.h +++ b/interface/RandStartPt.h @@ -38,8 +38,8 @@ class RandStartPt { 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::auto_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::auto_ptr ¶m, RooArgSet &snap, float &deltaNLL, double &nll_init, MultiDimFit::GridType gridType, double deltaX, double deltaY, CascadeMinimizer &minimObj); + 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 d4a4d08ba6f..14b12afb593 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -109,17 +109,17 @@ MultiDimFit::MultiDimFit() : ("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") - ("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") - ("randPointsRanges", boost::program_options::value(&randPointsRanges_)->default_value(""), "Range from which to draw random start points for the profiled POIs") - ; + ("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") + ("randPointsRanges", boost::program_options::value(&randPointsRanges_)->default_value(""), "Range from which to draw random start points for the profiled POIs") + ; } void MultiDimFit::applyOptions(const boost::program_options::variables_map &vm) @@ -644,14 +644,15 @@ 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; + 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************************************************************************************"< 1) std::cout << "\nStarting n==2. The nll0 from initial fit: " << nll0 << std::endl; + 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); diff --git a/src/RandStartPt.cc b/src/RandStartPt.cc index 554c855ba19..1adf021facf 100644 --- a/src/RandStartPt.cc +++ b/src/RandStartPt.cc @@ -146,7 +146,7 @@ void RandStartPt::setValSpecifiedObjs(){ } } -void RandStartPt::doRandomStartPt1DGridScan(double &xval, unsigned int poiSize, std::vector &poival, std::vector &poivars, std::auto_ptr ¶m, RooArgSet &snap, float &deltaNLL, double &nll_init, CascadeMinimizer &minimObj){ +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 (); @@ -181,7 +181,7 @@ void RandStartPt::doRandomStartPt1DGridScan(double &xval, unsigned int poiSize, } } -void RandStartPt::doRandomStartPt2DGridScan(double &xval, double &yval, unsigned int poiSize, std::vector &poival, std::vector &poivars, std::auto_ptr ¶m, RooArgSet &snap, float &deltaNLL, double &nll_init, MultiDimFit::GridType gridType, double deltaX, double deltaY, CascadeMinimizer &minimObj){ +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 (); From 2561169bcae4097a11c9cab4b8292fe6661f32fa Mon Sep 17 00:00:00 2001 From: Aashwin Basnet Date: Wed, 24 Jan 2024 21:11:16 +0100 Subject: [PATCH 36/37] edited behavior of random starting point range to not override POI ranges --- interface/MultiDimFit.h | 2 +- interface/RandStartPt.h | 4 ++-- src/MultiDimFit.cc | 9 ++++----- src/RandStartPt.cc | 31 ++++++++++++++++++------------- 4 files changed, 25 insertions(+), 21 deletions(-) diff --git a/interface/MultiDimFit.h b/interface/MultiDimFit.h index e963d02ba3b..91da1be399d 100644 --- a/interface/MultiDimFit.h +++ b/interface/MultiDimFit.h @@ -70,7 +70,7 @@ class MultiDimFit : public FitterAlgoBase { static std::string robustHesseSave_; static int pointsRandProf_; - static std::string randPointsRanges_; + static std::string setParameterRandomInitialValueRanges_; static int randPointsSeed_; static std::string saveSpecifiedFuncs_; diff --git a/interface/RandStartPt.h b/interface/RandStartPt.h index 8d41e823f20..a8e60518ffc 100644 --- a/interface/RandStartPt.h +++ b/interface/RandStartPt.h @@ -17,7 +17,7 @@ class RandStartPt { std::vector &specifiedvars_; std::vector &specifiedvals_; bool skipdefaultstart_; - std::string randptsranges_; + std::string parameterRandInitialValranges_; int numrandpts_; int verbosity_; bool fastscan_; @@ -32,7 +32,7 @@ class RandStartPt { std::vector &specifiedcatvals_; unsigned int nOtherFloatingPOI_; public: - RandStartPt(RooAbsReal& nll, std::vector &specifiedvars, std::vector &specifiedvals, bool skipdefaultstart, std::string randptsranges, 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); + 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); diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index d70e845481e..f527f27d194 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -64,7 +64,7 @@ std::string MultiDimFit::robustHesseLoad_ = ""; std::string MultiDimFit::robustHesseSave_ = ""; int MultiDimFit::pointsRandProf_ = 0; int MultiDimFit::randPointsSeed_ = 0; -std::string MultiDimFit::randPointsRanges_; +std::string MultiDimFit::setParameterRandomInitialValueRanges_; std::string MultiDimFit::saveSpecifiedFuncs_; @@ -119,7 +119,7 @@ MultiDimFit::MultiDimFit() : ("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") - ("randPointsRanges", boost::program_options::value(&randPointsRanges_)->default_value(""), "Range from which to draw random start points 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") ; } @@ -702,8 +702,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) ////////////////// 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 pointsRandProf_ to 0 or leave it unspecified. ///////////////////////////////////////// - - RandStartPt randStartPt(nll, specifiedVars_, specifiedVals_, skipDefaultStart_, randPointsRanges_, pointsRandProf_, verbose, fastScan_, hasMaxDeltaNLLForProf_, maxDeltaNLLForProf_, specifiedNuis_, specifiedFuncNames_, specifiedFunc_, specifiedFuncVals_, specifiedCatNames_, specifiedCat_, specifiedCatVals_, nOtherFloatingPoi_); + RandStartPt randStartPt(nll, specifiedVars_, specifiedVals_, skipDefaultStart_, setParameterRandomInitialValueRanges_, pointsRandProf_, 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 @@ -787,7 +786,7 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &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_, randPointsRanges_, pointsRandProf_, verbose, fastScan_, hasMaxDeltaNLLForProf_, maxDeltaNLLForProf_, specifiedNuis_, specifiedFuncNames_, specifiedFunc_, specifiedFuncVals_, specifiedCatNames_, specifiedCat_, specifiedCatVals_, nOtherFloatingPoi_); + 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 diff --git a/src/RandStartPt.cc b/src/RandStartPt.cc index 1adf021facf..621ab07f882 100644 --- a/src/RandStartPt.cc +++ b/src/RandStartPt.cc @@ -20,12 +20,12 @@ #include #include -RandStartPt::RandStartPt(RooAbsReal& nll, std::vector &specifiedvars, std::vector &specifiedvals, bool skipdefaultstart, std::string randptsranges, 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) : +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), - randptsranges_(randptsranges), + parameterRandInitialValranges_(parameterRandInitialValranges), numrandpts_(numrandpts), verbosity_(verbose), fastscan_(fastscan), @@ -57,16 +57,22 @@ std::vector> RandStartPt::vectorOfPointsToTry (){ // 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 (randptsranges_ != "") { - rand_ranges_dict = RandStartPt::getRangesDictFromInString(randptsranges_); + 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()][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)); + if (parameterRandInitialValranges_ != "") { + if (rand_ranges_dict.find(specifiedvars_[prof_param_idx]->GetName()) != 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; @@ -78,11 +84,13 @@ std::vector> RandStartPt::vectorOfPointsToTry (){ //Print vector of points to try if (verbosity_ > 1) { - std::cout<<"List of points to try: "<GetName() <<" = "<< val << std::endl; + index++; } } } @@ -125,9 +133,6 @@ void RandStartPt::setProfPOIvalues(unsigned int startptIdx, std::vector 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; - if (numrandpts_ > 0){ - specifiedvars_[var_idx]->removeRange(); //Do not impose a range if we are trying multiple random start points - } 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; From 25db367b02a2e093b421687927888027ffdc1985 Mon Sep 17 00:00:00 2001 From: nckw Date: Thu, 25 Jan 2024 14:16:08 +0000 Subject: [PATCH 37/37] Update MultiDimFit.cc Remove unecessary COUT and add logger command --- src/MultiDimFit.cc | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index f527f27d194..eaa83498ffe 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -646,13 +646,14 @@ void MultiDimFit::doGrid(RooWorkspace *w, RooAbsReal &nll) std::cout<<"\n************************************************************************************"<