Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
xRooFit.cxx
Go to the documentation of this file.
1/*
2 * Project: xRooFit
3 * Author:
4 * Will Buttinger, RAL 2022
5 *
6 * Copyright (c) 2022, CERN
7 *
8 * Redistribution and use in source and binary forms,
9 * with or without modification, are permitted according to the terms
10 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
11 */
12
13#include "RVersion.h"
14
15// #define private public
16// #include "Minuit2/Minuit2Minimizer.h"
17// #undef private
18// #include "Minuit2/FunctionMinimum.h"
19
20#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
21#define protected public
22#endif
23#include "RooFitResult.h"
24#ifdef protected
25#undef protected
26#endif
27
28#include "xRooFit/xRooFit.h"
29
30#include "RooDataSet.h"
31#include "RooSimultaneous.h"
32#include "RooArgSet.h"
33#include "RooRandom.h"
34#include "RooAbsPdf.h"
35#include "TUUID.h"
36#include "RooProdPdf.h"
37#include "RooGamma.h"
38#include "RooPoisson.h"
39#include "RooGaussian.h"
40#include "RooBifurGauss.h"
41#include "RooLognormal.h"
42#include "RooBinning.h"
43#include "RooUniformBinning.h"
44
46#include "Math/GenAlgoOptions.h"
47#include "Math/Minimizer.h"
48#include "RooMinimizer.h"
49#include "coutCapture.h"
50
51#include "TCanvas.h"
52#include "TGraphErrors.h"
53#include "TLegend.h"
54#include "TKey.h"
55#include "RooAbsTestStatistic.h"
56#include "TPRegexp.h"
57#include "RooStringVar.h"
58
59#include "RooRealProxy.h"
60#include "RooSuperCategory.h"
61
62#include "xRooFitVersion.h"
63
64#include <csignal>
65#include "TROOT.h"
66#include "TBrowser.h"
67
69
70std::shared_ptr<RooLinkedList> xRooFit::sDefaultNLLOptions = nullptr;
71std::shared_ptr<ROOT::Fit::FitConfig> xRooFit::sDefaultFitConfig = nullptr;
72
73const char *xRooFit::GetVersion()
74{
75 return GIT_COMMIT_HASH;
76}
77const char *xRooFit::GetVersionDate()
78{
79 return GIT_COMMIT_DATE;
80}
81
82RooCmdArg xRooFit::ReuseNLL(bool flag)
83{
84 return RooCmdArg("ReuseNLL", flag, 0, 0, 0, nullptr, nullptr, nullptr, nullptr);
85}
86
87RooCmdArg xRooFit::Tolerance(double val)
88{
89 return RooCmdArg("Tolerance", 0, 0, val);
90}
91
92RooCmdArg xRooFit::StrategySequence(const char *val)
93{
94 return RooCmdArg("StrategySequence", 0, 0, 0, 0, val);
95}
96
97xRooNLLVar xRooFit::createNLL(const std::shared_ptr<RooAbsPdf> pdf, const std::shared_ptr<RooAbsData> data,
98 const RooLinkedList &nllOpts)
99{
100 return xRooNLLVar(pdf, data, nllOpts);
101}
102
103xRooNLLVar xRooFit::createNLL(RooAbsPdf &pdf, RooAbsData *data, const RooLinkedList &nllOpts)
104{
105 return createNLL(std::shared_ptr<RooAbsPdf>(&pdf, [](RooAbsPdf *) {}),
106 std::shared_ptr<RooAbsData>(data, [](RooAbsData *) {}), nllOpts);
107}
108
109xRooNLLVar xRooFit::createNLL(RooAbsPdf &pdf, RooAbsData *data, const RooCmdArg &arg1, const RooCmdArg &arg2,
110 const RooCmdArg &arg3, const RooCmdArg &arg4, const RooCmdArg &arg5,
111 const RooCmdArg &arg6, const RooCmdArg &arg7, const RooCmdArg &arg8)
112{
113
115 l.Add((TObject *)&arg1);
116 l.Add((TObject *)&arg2);
117 l.Add((TObject *)&arg3);
118 l.Add((TObject *)&arg4);
119 l.Add((TObject *)&arg5);
120 l.Add((TObject *)&arg6);
121 l.Add((TObject *)&arg7);
122 l.Add((TObject *)&arg8);
123 return createNLL(pdf, data, l);
124}
125
126std::shared_ptr<const RooFitResult>
127xRooFit::fitTo(RooAbsPdf &pdf,
128 const std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> &data,
129 const RooLinkedList &nllOpts, const ROOT::Fit::FitConfig &fitConf)
130{
131 return xRooNLLVar(std::shared_ptr<RooAbsPdf>(&pdf, [](RooAbsPdf *) {}), data, nllOpts)
132 .minimize(std::shared_ptr<ROOT::Fit::FitConfig>(const_cast<ROOT::Fit::FitConfig *>(&fitConf),
133 [](ROOT::Fit::FitConfig *) {}));
134}
135
136std::shared_ptr<const RooFitResult> xRooFit::fitTo(RooAbsPdf &pdf,
137 const std::pair<RooAbsData *, const RooAbsCollection *> &data,
138 const RooLinkedList &nllOpts, const ROOT::Fit::FitConfig &fitConf)
139{
140 return xRooNLLVar(pdf, data, nllOpts)
141 .minimize(std::shared_ptr<ROOT::Fit::FitConfig>(const_cast<ROOT::Fit::FitConfig *>(&fitConf),
142 [](ROOT::Fit::FitConfig *) {}));
143}
144
145std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>>
146xRooFit::generateFrom(RooAbsPdf &pdf, const RooFitResult &_fr, bool expected, int seed)
147{
148
149 std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> out;
150
151 auto fr = &_fr;
152 if (!fr)
153 return out;
154
155 auto _allVars = std::unique_ptr<RooAbsCollection>(pdf.getVariables());
156 auto _snap = std::unique_ptr<RooAbsCollection>(_allVars->snapshot());
157 *_allVars = fr->constPars();
158 *_allVars = fr->floatParsFinal();
159
160 // determine globs from fr constPars
161 auto _globs = std::unique_ptr<RooAbsCollection>(fr->constPars().selectByAttrib("global", true));
162
163 // bool doBinned = false;
164 // RooAbsPdf::GenSpec** gs = nullptr;
165
166 if (seed == 0)
167 seed = RooRandom::randomGenerator()->Integer(std::numeric_limits<uint32_t>::max());
169
170 TString uuid = TUUID().AsString();
171
172 std::function<std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooArgSet>>(RooAbsPdf *)> genSubPdf;
173
174 genSubPdf = [&](RooAbsPdf *_pdf) {
175 std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooArgSet>> _out;
176 // std::unique_ptr<RooArgSet> _obs(_pdf->getParameters(*pars)); // using this "trick" to get observables can
177 // produce 'error' msg because of RooProdPdf trying to setup partial integrals
178 std::unique_ptr<RooArgSet> _obs(_pdf->getVariables());
179 _obs->remove(fr->constPars(), true, true);
180 _obs->remove(fr->floatParsFinal(), true, true); // use this instead
181
182 if (!_globs->empty()) {
183 RooArgSet *toy_gobs = new RooArgSet(uuid + "_globs");
184 // ensure we use the gobs from the model ...
185 RooArgSet t;
186 t.add(*_globs);
187 std::unique_ptr<RooArgSet> globs(_pdf->getObservables(t));
188 globs->snapshot(*toy_gobs);
189 if (!toy_gobs->empty() &&
190 !dynamic_cast<RooSimultaneous *>(
191 _pdf)) { // if was simPdf will call genSubPdf on each subpdf so no need to generate here
192 if (!expected) {
193 *toy_gobs = *std::unique_ptr<RooDataSet>(_pdf->generate(*globs, 1))->get();
194 } else {
195 // loop over pdfs in top-level prod-pdf,
196 auto pp = dynamic_cast<RooProdPdf *>(_pdf);
197 if (pp) {
198 for (auto thePdf : pp->pdfList()) {
199 auto gob = std::unique_ptr<RooArgSet>(thePdf->getObservables(*globs));
200 if (gob->empty())
201 continue;
202 if (gob->size() > 1) {
203 Warning("generate", "%s contains multiple global obs: %s", thePdf->GetName(),
204 gob->contentsString().c_str());
205 continue;
206 }
207 RooRealVar &rrv = dynamic_cast<RooRealVar &>(*gob->first());
208 std::unique_ptr<RooArgSet> cpars(thePdf->getParameters(*globs));
209
210 bool foundServer = false;
211 // note : this will work only for this type of constraints
212 // expressed as RooPoisson, RooGaussian, RooLognormal, RooGamma
213 // SimpleGaussianConstraint is CMS's own version of a RooGaussian, which also works.
214 TClass *cClass = thePdf->IsA();
215 if (cClass != RooGaussian::Class() && cClass != RooPoisson::Class() &&
216 cClass != RooGamma::Class() && cClass != RooLognormal::Class() &&
217 cClass != RooBifurGauss::Class() &&
218 !(cClass && strcmp(cClass->GetName(), "SimpleGaussianConstraint") == 0)) {
219 TString className = (cClass) ? cClass->GetName() : "undefined";
220 oocoutW((TObject *)nullptr, Generation)
221 << "AsymptoticCalculator::MakeAsimovData:constraint term " << thePdf->GetName()
222 << " of type " << className << " is a non-supported type - result might be not correct "
223 << std::endl;
224 }
225
226 // in case of a Poisson constraint make sure the rounding is not set
227 if (cClass == RooPoisson::Class()) {
228 RooPoisson *pois = dynamic_cast<RooPoisson *>(thePdf);
229 assert(pois);
230 pois->setNoRounding(true);
231 }
232
233 // look at server of the constraint term and check if the global observable is part of the server
234 RooAbsArg *arg = thePdf->findServer(rrv);
235 if (!arg) {
236 // special case is for the Gamma where one might define the global observable n and you have a
237 // Gamma(b, n+1, ...._ in this case n+1 is the server and we don;t have a direct dependency, but
238 // we want to set n to the b value so in case of the Gamma ignore this test
239 if (cClass != RooGamma::Class()) {
240 oocoutE((TObject *)nullptr, Generation)
241 << "AsymptoticCalculator::MakeAsimovData:constraint term " << thePdf->GetName()
242 << " has no direct dependence on global observable- cannot generate it " << std::endl;
243 continue;
244 }
245 }
246
247 // loop on the server of the constraint term
248 // need to treat the Gamma as a special case
249 // the mode of the Gamma is (k-1)*theta where theta is the inverse of the rate parameter.
250 // we assume that the global observable is defined as ngobs = k-1 and the theta parameter has the
251 // name theta otherwise we use other procedure which might be wrong
252 RooAbsReal *thetaGamma = nullptr;
253 if (cClass == RooGamma::Class()) {
254 for (RooAbsArg *a2 : thePdf->servers()) {
255 if (TString(a2->GetName()).Contains("theta")) {
256 thetaGamma = dynamic_cast<RooAbsReal *>(a2);
257 break;
258 }
259 }
260 if (thetaGamma == nullptr) {
261 oocoutI((TObject *)nullptr, Generation)
262 << "AsymptoticCalculator::MakeAsimovData:constraint term " << thePdf->GetName()
263 << " is a Gamma distribution and no server named theta is found. Assume that the Gamma "
264 "scale is 1 "
265 << std::endl;
266 }
267 }
268 for (RooAbsArg *a2 : thePdf->servers()) {
269 RooAbsReal *rrv2 = dynamic_cast<RooAbsReal *>(a2);
270 if (rrv2 && !rrv2->dependsOn(*gob) &&
271 (!rrv2->isConstant() || !rrv2->InheritsFrom("RooConstVar"))) {
272
273 // found server not depending on the gob
274 if (foundServer) {
275 oocoutE((TObject *)nullptr, Generation)
276 << "AsymptoticCalculator::MakeAsimovData:constraint term " << thePdf->GetName()
277 << " constraint term has more server depending on nuisance- cannot generate it "
278 << std::endl;
279 foundServer = false;
280 break;
281 }
282 if (thetaGamma && thetaGamma->getVal() > 0) {
283 rrv.setVal(rrv2->getVal() / thetaGamma->getVal());
284 } else {
285 rrv.setVal(rrv2->getVal());
286 }
287 foundServer = true;
288 }
289 }
290
291 if (!foundServer) {
292 oocoutE((TObject *)nullptr, Generation)
293 << "AsymptoticCalculator::MakeAsimovData - can't find nuisance for constraint term - global "
294 "observables will not be set to Asimov value "
295 << thePdf->GetName() << std::endl;
296 std::cerr << "Parameters: " << std::endl;
297 cpars->Print("V");
298 std::cerr << "Observables: " << std::endl;
299 gob->Print("V");
300 }
301 }
302 } else {
303 Error("generate", "Cannot generate global observables, pdf is: %s::%s", _pdf->ClassName(),
304 _pdf->GetName());
305 }
306 *toy_gobs = *globs;
307 }
308 }
309 _out.second.reset(toy_gobs);
310 } // end of globs generation
311
312 RooRealVar w("weightVar", "weightVar", 1);
313 if (auto s = dynamic_cast<RooSimultaneous *>(_pdf)) {
314 // do subpdf's individually
315 _obs->add(w);
316 _out.first = std::make_unique<RooDataSet>(
317 uuid, TString::Format("%s %s", _pdf->GetTitle(), (expected) ? "Expected" : "Toy"), *_obs,
318 RooFit::WeightVar("weightVar"));
319
320 for (auto &c : s->indexCat()) {
321#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 22, 00)
322 std::string cLabel = c.first;
323#else
324 std::string cLabel = c->GetName();
325#endif
326 auto p = s->getPdf(cLabel.c_str());
327 if (!p)
328 continue;
329 auto toy = genSubPdf(p);
330 if (toy.second && _out.second)
331 *const_cast<RooArgSet *>(_out.second.get()) = *toy.second;
332 _obs->setCatLabel(s->indexCat().GetName(), cLabel.c_str());
333 for (int i = 0; i < toy.first->numEntries(); i++) {
334 *_obs = *toy.first->get(i);
335 _out.first->add(*_obs, toy.first->weight());
336 }
337 }
338 return _out;
339 }
340
341 std::map<RooRealVar *, std::shared_ptr<RooAbsBinning>> binnings;
342
343 for (auto &o : *_obs) {
344 auto r = dynamic_cast<RooRealVar *>(o);
345 if (!r)
346 continue;
347 if (auto res = _pdf->binBoundaries(*r, -std::numeric_limits<double>::infinity(),
348 std::numeric_limits<double>::infinity())) {
349 binnings[r] = std::shared_ptr<RooAbsBinning>(r->getBinning().clone(r->getBinning().GetName()));
350
351 std::vector<double> boundaries;
352 boundaries.reserve(res->size());
353 for (auto &rr : *res) {
354 if (boundaries.empty() || std::abs(boundaries.back() - rr) > 1e-3 ||
355 std::abs(boundaries.back() - rr) > 1e-5 * boundaries.back())
356 boundaries.push_back(rr);
357 } // sometimes get virtual duplicates of boundaries
358 r->setBinning(RooBinning(boundaries.size() - 1, &boundaries[0]));
359 delete res;
360 } else if (r->numBins(r->getBinning().GetName()) == 0 && expected) {
361 // no bins ... in order to generate expected we need to have some bins
362 binnings[r] = std::shared_ptr<RooAbsBinning>(r->getBinning().clone(r->getBinning().GetName()));
363 r->setBinning(RooUniformBinning(r->getMin(), r->getMax(), 100));
364 }
365 }
366
367 // now can generate
368 if (_obs->empty()) {
369 // no observables, create a single dataset with 1 entry ... why 1 entry??
370 _obs->add(w);
371 RooArgSet _tmp;
372 _tmp.add(w);
373 _out.first = std::make_unique<RooDataSet>("", "Toy", _tmp, RooFit::WeightVar("weightVar"));
374 _out.first->add(_tmp);
375 } else {
376 if (_pdf->canBeExtended()) {
377 _out.first =
378 std::unique_ptr<RooDataSet>{_pdf->generate(*_obs, RooFit::Extended(), RooFit::ExpectedData(expected))};
379 } else {
380 if (expected) {
381 // use AsymptoticCalculator because generate expected not working correctly on unextended pdf?
382 // TODO: Can the above code for expected globs be used instead, or what about replace above code with
383 // ObsToExpected?
384 _out.first.reset(RooStats::AsymptoticCalculator::GenerateAsimovData(*_pdf, *_obs));
385 } else {
386 _out.first = std::unique_ptr<RooDataSet>{_pdf->generate(*_obs, RooFit::ExpectedData(expected))};
387 }
388 }
389 }
390 _out.first->SetName(TUUID().AsString());
391
392 for (auto &b : binnings) {
393 auto v = b.first;
394 auto binning = b.second;
395 v->setBinning(*binning);
396 // range of variable in dataset may be less than in the workspace
397 // if e.g. generate for a specific channel. So need to expand ranges to match
398 auto x = dynamic_cast<RooRealVar *>(_out.first->get()->find(v->GetName()));
399 auto r = x->getRange();
400 if (r.first > binning->lowBound())
401 x->setMin(binning->lowBound());
402 if (r.second < binning->highBound())
403 x->setMax(binning->highBound());
404 }
405 return _out;
406 };
407
408 out = genSubPdf(&pdf);
409 out.first->SetName(expected ? (TString(fr->GetName()) + "_asimov") : uuid);
410
411 // from now on we store the globs in the dataset
412 if (out.second) {
413 out.first->setGlobalObservables(*out.second);
414 out.second.reset();
415 }
416
417#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 26, 00)
418 // store fitResult name on the weightVar
419 if (auto w = dynamic_cast<RooDataSet *>(out.first.get())->weightVar()) {
420 w->setStringAttribute("fitResult", fr->GetName());
421 w->setAttribute("expected", expected);
422 }
423#endif
424
425 *_allVars = *_snap;
426
427 // June2023: Added this because found that generation was otherwise getting progressively slower
428 // the RooAbsPdf::generate does a clone, and it seems that the RooCacheManager of the original pdf
429 // is getting polluted on each generate call, causing it to grow larger and therefore the clone of it
430 // to take longer and longer. So sterilize to clear the caches of all components
431#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
432 auto _ws = pdf._myws;
433#else
434 auto _ws = pdf.workspace();
435#endif
436 if (_ws) {
437 // do explicitly rather than via xRooNode sterilize method because don't want to invoke the constructor
438 // workspace tweaking features (which sets poi etc etc)
439 for (auto obj : _ws->components()) {
440 for (int i = 0; i < obj->numCaches(); i++) {
441 if (auto cache = dynamic_cast<RooObjCacheManager *>(obj->getCache(i))) {
442 cache->reset();
443 }
444 }
445 if (RooAbsPdf *p = dynamic_cast<RooAbsPdf *>(obj); p) {
446 p->setNormRange(p->normRange());
447 }
448 obj->setValueDirty();
449 }
450 // xRooNode(pdf.workspace()).sterilize();
451 }
452
453 return out;
454}
455
456std::shared_ptr<RooLinkedList> xRooFit::createNLLOptions()
457{
458 auto out = std::shared_ptr<RooLinkedList>(new RooLinkedList, [](RooLinkedList *l) {
459 l->Delete();
460 delete l;
461 });
462 for (auto opt : *defaultNLLOptions()) {
463 out->Add(opt->Clone(nullptr)); // nullptr needed because accessing Clone via TObject base class puts
464 // "" instead, so doesnt copy names
465 }
466 return out;
467}
468
469std::shared_ptr<RooLinkedList> xRooFit::defaultNLLOptions()
470{
471 if (sDefaultNLLOptions)
472 return sDefaultNLLOptions;
473 sDefaultNLLOptions = std::shared_ptr<RooLinkedList>(new RooLinkedList, [](RooLinkedList *l) {
474 l->Delete();
475 delete l;
476 });
477 sDefaultNLLOptions->Add(RooFit::Offset().Clone());
478 // disable const-optimization at the construction step ... can happen in the minimization though
479 sDefaultNLLOptions->Add(RooFit::Optimize(0).Clone());
480 return sDefaultNLLOptions;
481}
482
483std::shared_ptr<ROOT::Fit::FitConfig> xRooFit::createFitConfig()
484{
485 return std::make_shared<ROOT::Fit::FitConfig>(*defaultFitConfig());
486}
487
488std::shared_ptr<ROOT::Fit::FitConfig> xRooFit::defaultFitConfig()
489{
490 if (sDefaultFitConfig)
491 return sDefaultFitConfig;
492 sDefaultFitConfig = std::make_shared<ROOT::Fit::FitConfig>();
493 auto &fitConfig = *sDefaultFitConfig;
494 fitConfig.SetParabErrors(true); // will use to run hesse after fit
495 fitConfig.MinimizerOptions().SetMinimizerType("Minuit2");
496 fitConfig.MinimizerOptions().SetErrorDef(0.5); // ensures errors are +/- 1 sigma ..IMPORTANT
497 fitConfig.SetParabErrors(true); // runs HESSE
498 fitConfig.SetMinosErrors(true); // computes asymmetric errors on any parameter with the "minos" attribute set
499 fitConfig.MinimizerOptions().SetMaxFunctionCalls(
500 -1); // calls per iteration. if left as 0 will set automatically to 500*nPars below
501 fitConfig.MinimizerOptions().SetMaxIterations(-1); // if left as 0 will set automatically to 500*nPars
502 fitConfig.MinimizerOptions().SetStrategy(-1); // will start at front of StrategySequence (given below)
503 // fitConfig.MinimizerOptions().SetTolerance(
504 // 1); // default is 0.01 (i think) but roominimizer uses 1 as default - use specify with
505 // ROOT::Math::MinimizerOptions::SetDefaultTolerance(..)
506 fitConfig.MinimizerOptions().SetPrintLevel(-2);
507 fitConfig.MinimizerOptions().SetExtraOptions(ROOT::Math::GenAlgoOptions());
508 // have to const cast to set extra options
509 auto extraOpts = const_cast<ROOT::Math::IOptions *>(fitConfig.MinimizerOptions().ExtraOptions());
510 extraOpts->SetValue("OptimizeConst", 2); // if 0 will disable constant term optimization and cache-and-track of the
511 // NLL. 1 = just caching, 2 = cache and track
512#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 29, 00)
513 extraOpts->SetValue("StrategySequence", "0s01s12s2s3m");
514 extraOpts->SetValue("HesseStrategySequence", "23");
515#else
516 extraOpts->SetValue("StrategySequence", "0s01s12s2m");
517 extraOpts->SetValue("HesseStrategySequence", "2");
518#endif
519 extraOpts->SetValue(
520 "HesseStrategy",
521 -1); // when hesse is run after minimization, will use this strategy. -1 means start at begin of strat sequence
522 extraOpts->SetValue("LogSize", 0); // length of log to capture and save
523 extraOpts->SetValue("BoundaryCheck",
524 0.); // if non-zero, warn if any post-fit value is close to boundary (e.g. 0.01 = within 1%)
525 extraOpts->SetValue("TrackProgress", 30); // seconds between output to log of evaluation progress
526 extraOpts->SetValue("xRooFitVersion", GIT_COMMIT_HASH); // not really options but here for logging purposes
527 // extraOpts->SetValue("ROOTVersion",ROOT_VERSION_CODE); - not needed as should by part of the ROOT TFile definition
528
529 // extraOpts->SetValue("HessianStepTolerance",0.);
530 // extraOpts->SetValue("HessianG2Tolerance",0.);
531
532 return sDefaultFitConfig;
533}
534
535ROOT::Math::IOptions *xRooFit::defaultFitConfigOptions()
536{
537 return const_cast<ROOT::Math::IOptions *>(defaultFitConfig()->MinimizerOptions().ExtraOptions());
538}
539
541public:
542 void (*oldHandlerr)(int) = nullptr;
544 static bool fInterrupt;
545 static void interruptHandler(int signum)
546 {
547 if (signum == SIGINT) {
548 std::cout << "Minimization interrupted ... will exit as soon as possible" << std::endl;
549 // TODO: create a global mutex for this
550 fInterrupt = true;
551 } else {
552 if (me)
553 me->oldHandlerr(signum);
554 }
555 };
556 ProgressMonitor(RooAbsReal &f, int interval = 30)
557 : RooAbsReal(Form("progress_%s", f.GetName()), ""),
558 oldHandlerr(signal(SIGINT, interruptHandler)),
559 fFunc("func", "func", this, f),
560 fInterval(interval)
561 {
562 s.Start();
563
564 me = this;
565 vars.reset(std::unique_ptr<RooAbsCollection>(f.getVariables())->selectByAttrib("Constant", false));
566 }
568 {
569 if (oldHandlerr) {
570 signal(SIGINT, oldHandlerr);
571 }
572 if (me == this)
573 me = nullptr;
574 };
575 ProgressMonitor(const ProgressMonitor &other, const char *name = nullptr)
576 : RooAbsReal(other, name), fFunc("func", this, other.fFunc), fInterval(other.fInterval)
577 {
578 }
579 TObject *clone(const char *newname) const override { return new ProgressMonitor(*this, newname); }
580
581 // required forwarding methods for RooEvaluatorWrapper in 6.32 onwards
582 double defaultErrorLevel() const override { return fFunc->defaultErrorLevel(); }
583 bool getParameters(const RooArgSet *observables, RooArgSet &outputSet, bool stripDisconnected) const override
584 {
585 return fFunc->getParameters(observables, outputSet, stripDisconnected);
586 }
587 bool setData(RooAbsData &data, bool cloneData) override { return fFunc->setData(data, cloneData); }
588 double getValV(const RooArgSet *) const override { return evaluate(); }
589 void applyWeightSquared(bool flag) override { fFunc->applyWeightSquared(flag); }
590 void printMultiline(std::ostream &os, Int_t contents, bool verbose = false, TString indent = "") const override
591 {
592 fFunc->printMultiline(os, contents, verbose, indent);
593 }
594 void constOptimizeTestStatistic(ConstOpCode opcode, bool doAlsoTrackingOpt) override
595 {
596 fFunc->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt);
597 }
598
599 double evaluate() const override
600 {
601 if (fInterrupt) {
602 throw std::runtime_error("Keyboard interrupt");
603 return std::numeric_limits<double>::quiet_NaN();
604 }
605 double out = fFunc;
606 if (prevMin == std::numeric_limits<double>::infinity()) {
607 prevMin = out;
609 }
610 if (!std::isnan(out)) {
611 if (out < minVal) {
612 if (minPars.empty())
614 minPars = *vars;
615 }
616 minVal = std::min(minVal, out);
617 }
618 counter++;
619 if (s.RealTime() > fInterval) {
620 double evalRate = (counter - prevCounter) / s.RealTime();
621 s.Reset();
622 std::stringstream sout;
623
624 sout << (counter) << ") (" << evalRate << "Hz) " << TDatime().AsString();
625 if (!fState.empty())
626 sout << " : " << fState;
627 if (counter2) {
628 // doing a hesse step, estimate progress based on evaluations
629 int nRequired = prevPars.size();
630 if (nRequired > 1) {
631 nRequired *= (nRequired - 1) / 2;
632 if (fState == "Hesse3") {
633 nRequired *= 4;
634 }
635 sout << " (~" << int(100.0 * (counter - counter2) / nRequired) << "%)";
636 }
637 }
638 sout << " : " << minVal << " Delta = " << (minVal - prevMin);
639 if (minVal < prevMin) {
640 sout << " : ";
641 // compare minPars and prevPars, print biggest deltas
642 std::vector<std::pair<double, std::string>> parDeltas;
643 parDeltas.reserve(minPars.size());
644 for (auto p : minPars) {
645 parDeltas.emplace_back(std::pair<double, std::string>(
646 dynamic_cast<RooRealVar *>(p)->getVal() - prevPars.getRealValue(p->GetName()), p->GetName()));
647 }
648 std::sort(parDeltas.begin(), parDeltas.end(),
649 [](auto &left, auto &right) { return std::abs(left.first) > std::abs(right.first); });
650 int i;
651 for (i = 0; i < std::min(3, int(parDeltas.size())); i++) {
652 if (parDeltas.at(i).first == 0)
653 break;
654 if (i != 0)
655 sout << ",";
656 sout << parDeltas.at(i).second << (parDeltas.at(i).first >= 0 ? "+" : "-") << "="
657 << std::abs(parDeltas.at(i).first) << "(" << minPars.getRealValue(parDeltas.at(i).second.c_str())
658 << ")";
659 }
660 if (i < int(parDeltas.size()) && parDeltas.at(i).first != 0)
661 sout << " ...";
663 }
664
665 if (gROOT->FromPopUp() && gROOT->GetListOfBrowsers()->At(0)) {
666 auto browser = dynamic_cast<TBrowser *>(gROOT->GetListOfBrowsers()->At(0));
667 std::string status = sout.str();
668 int col = 0;
669 while (col < 4) {
670 std::string status_part;
671 if (status.find(" : ") != std::string::npos) {
672 status_part = status.substr(0, status.find(" : "));
673 status = status.substr(status.find(" : ") + 3);
674 } else {
675 status_part = status;
676 status = "";
677 }
678 browser->SetStatusText(status_part.c_str(), col);
679 col++;
680 }
682 }
683 std::cerr << sout.str() << std::endl;
684
685 prevMin = minVal;
687 } else {
688 s.Continue();
689 }
690 return out;
691 }
692
693 std::string fState;
694 mutable int counter = 0;
695 int counter2 = 0; // used to estimate progress of a Hesse calculation
696
697private:
699 mutable double minVal = std::numeric_limits<double>::infinity();
700 mutable double prevMin = std::numeric_limits<double>::infinity();
703 mutable int prevCounter = 0;
704 mutable int fInterval = 0; // time in seconds before next report
705 mutable TStopwatch s;
706 std::shared_ptr<RooAbsCollection> vars;
707};
708bool ProgressMonitor::fInterrupt = false;
710
711xRooFit::StoredFitResult::StoredFitResult(RooFitResult *_fr) : TNamed(*_fr)
712{
713 fr.reset(_fr);
714}
715
716xRooFit::StoredFitResult::StoredFitResult(const std::shared_ptr<RooFitResult> &_fr) : TNamed(*_fr), fr(_fr) {}
717
718std::shared_ptr<const RooFitResult> xRooFit::minimize(RooAbsReal &nll,
719 const std::shared_ptr<ROOT::Fit::FitConfig> &_fitConfig,
720 const std::shared_ptr<RooLinkedList> &nllOpts)
721{
722
723 auto myFitConfig = _fitConfig ? _fitConfig : createFitConfig();
724 auto &fitConfig = *myFitConfig;
725
726 auto _nll = &nll;
727
728 TString resultTitle = nll.getStringAttribute("fitresultTitle");
729 TString fitName = TUUID().AsString();
730 if (resultTitle == "")
731 resultTitle = TUUID(fitName).GetTime().AsString();
732
733 // extract any user pars from the nll too
734 RooArgList fUserPars;
735 if (nll.getStringAttribute("userPars")) {
736 TStringToken st(nll.getStringAttribute("userPars"), ",");
737 while (st.NextToken()) {
738 TString parName = st;
739 TString parVal = nll.getStringAttribute(parName);
740 if (parVal.IsFloat()) {
741 fUserPars.addClone(RooRealVar(parName, parName, parVal.Atof()));
742 } else {
743 fUserPars.addClone(RooStringVar(parName, parName, parVal));
744 }
745 }
746 }
747
748 auto _nllVars = std::unique_ptr<RooAbsCollection>(_nll->getVariables());
749
750 std::unique_ptr<RooAbsCollection> constPars(_nllVars->selectByAttrib("Constant", true));
751 constPars->add(fUserPars, true); // add here so checked for when loading from cache
752 std::unique_ptr<RooAbsCollection> floatPars(_nllVars->selectByAttrib("Constant", false));
753
754 int _progress = 0;
755 double boundaryCheck = 0;
756 std::string s;
757 std::string hs;
758 int logSize = 0;
759#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 29, 00)
760 int hesseStrategy = 3; // uses most precise hesse settings (step sizes and g2 tolerances)
761#else
762 int hesseStrategy = 2; // uses most precise hesse settings (step sizes and g2 tolerances)
763#endif
764 if (fitConfig.MinimizerOptions().ExtraOptions()) {
765 fitConfig.MinimizerOptions().ExtraOptions()->GetNamedValue("StrategySequence", s);
766 fitConfig.MinimizerOptions().ExtraOptions()->GetIntValue("TrackProgress", _progress);
767 fitConfig.MinimizerOptions().ExtraOptions()->GetRealValue("BoundaryCheck", boundaryCheck);
768 fitConfig.MinimizerOptions().ExtraOptions()->GetIntValue("LogSize", logSize);
769 fitConfig.MinimizerOptions().ExtraOptions()->GetIntValue("HesseStrategy", hesseStrategy);
770 fitConfig.MinimizerOptions().ExtraOptions()->GetNamedValue("HesseStrategySequence", hs);
771 }
772 TString m_strategy = s;
773 TString m_hessestrategy = hs;
774
775 // if fit caching enabled, try to locate a valid fitResult
776 // must have matching constPars
777 TDirectory *cacheDir = gDirectory;
778
779 if (cacheDir) {
780 if (auto nllDir = cacheDir->GetDirectory(nll.GetName()); nllDir) {
781 if (auto keys = nllDir->GetListOfKeys(); keys) {
782 for (auto &&k : *keys) {
783 auto cl = TClass::GetClass((static_cast<TKey *>(k))->GetClassName());
784 if (cl->InheritsFrom("RooFitResult")) {
785 StoredFitResult *storedFr =
786 nllDir->GetList() ? dynamic_cast<StoredFitResult *>(nllDir->GetList()->FindObject(k->GetName()))
787 : nullptr;
788 if (auto cachedFit =
789 (storedFr) ? storedFr->fr.get() : dynamic_cast<TKey *>(k)->ReadObject<RooFitResult>();
790 cachedFit) {
791 if (!storedFr) {
792 storedFr = new StoredFitResult(cachedFit);
793 nllDir->Add(storedFr);
794 // std::cout << "Loaded " << nllDir->GetPath() << "/" << k->GetName() << " : " << k->GetTitle()
795 // << std::endl;
796 }
797 bool match = true;
798 if (!cachedFit->floatParsFinal().equals(*floatPars)) {
799 match = false;
800 } else {
801 for (auto &p : *constPars) {
802 auto v = dynamic_cast<RooAbsReal *>(p);
803 if (!v) {
804 if (auto c = dynamic_cast<RooAbsCategory *>(p)) {
805 if (auto _p =
806 dynamic_cast<RooAbsCategory *>(cachedFit->constPars().find(p->GetName()));
807 _p && !_p->getAttribute("global") &&
808 _p->getCurrentIndex() != c->getCurrentIndex()) {
809 match = false;
810 break;
811 }
812 } else {
813 match = false;
814 break;
815 }
816 };
817 if (auto _p = dynamic_cast<RooAbsReal *>(cachedFit->constPars().find(p->GetName())); _p) {
818 // note: do not need global observable values to match (globals currently added to
819 // constPars list)
820 if (!_p->getAttribute("global") && std::abs(_p->getVal() - v->getVal()) > 1e-12) {
821 match = false;
822 break;
823 }
824 }
825 }
826 }
827 if (match) {
828 return storedFr->fr;
829 // return std::shared_ptr<RooFitResult>(cachedFit,[](RooFitResult*){}); // dir owns the
830 // fitResult - this means dir needs to stay open for fits to be valid return
831 // std::make_shared<RooFitResult>(*cachedFit); // return a copy ... dir doesn't need to stay
832 // open, but fit result isn't shared
833 } else {
834 // delete cachedFit;
835 }
836 }
837 }
838 }
839 }
840 }
841 }
842
843 if (nll.getAttribute("readOnly"))
844 return nullptr;
845
846 int printLevel = fitConfig.MinimizerOptions().PrintLevel();
848 if (printLevel < 0)
850
851 // check how many parameters we have ... if 0 parameters then we wont run a fit, we just evaluate nll and return ...
852 if (floatPars->empty() || fitConfig.MinimizerOptions().MaxFunctionCalls() == 1) {
853 std::shared_ptr<RooFitResult> result;
854 RooArgList parsList;
855 parsList.add(*floatPars);
856 // construct an empty fit result ...
857 result = std::make_shared<RooFitResult>(); // if put name here fitresult gets added to dir, we don't want that
858 result->SetName(TUUID().AsString());
859 result->SetTitle(resultTitle);
860 result->setFinalParList(parsList);
861 result->setInitParList(parsList);
862 result->setConstParList(dynamic_cast<RooArgSet &>(*constPars)); /* RooFitResult takes a snapshot */
864 d.ResizeTo(parsList.size(), parsList.size());
865 result->setCovarianceMatrix(d);
866 result->setCovQual(floatPars->empty() ? 3 : -1);
867 result->setMinNLL(_nll->getVal());
868 result->setEDM(0);
869 result->setStatus(floatPars->empty() ? 0 : 1);
870
871 std::vector<std::pair<std::string, int>> statusHistory;
872 statusHistory.emplace_back(std::make_pair("EVAL", result->status()));
873 result->setStatusHistory(statusHistory);
874
875 if (cacheDir && cacheDir->IsWritable()) {
876 // save a copy of fit result to relevant dir
877 if (!cacheDir->GetDirectory(nll.GetName()))
878 cacheDir->mkdir(nll.GetName());
879 if (auto dir = cacheDir->GetDirectory(nll.GetName()); dir) {
880 // save NLL opts if was given one, unless already present
881 if (nllOpts) {
882 if (strlen(nllOpts->GetName()) == 0) {
883 nllOpts->SetName(TUUID().AsString());
884 }
885 if (!dir->FindKey(nllOpts->GetName())) {
886 dir->WriteObject(nllOpts.get(), nllOpts->GetName());
887 }
888 }
889 dir->WriteObject(result.get(), result->GetName());
890 }
891 }
892
893 if (printLevel < 0)
895 return result;
896 }
897
898 std::shared_ptr<RooFitResult> out;
899
900 // check if any floatPars are categorical .. if so, need to a "discrete minimization" over the permutations
901 RooArgSet floatCats;
902 for (auto p : *floatPars) {
903 if (p->isCategory()) {
904 floatCats.add(*p);
905 }
906 }
907 if (!floatCats.empty()) {
908 RooSuperCategory allCats("floatCats", "Floating categorical parameters", floatCats);
909 std::unique_ptr<RooAbsCollection> _snap(floatCats.snapshot());
910 floatCats.setAttribAll("Constant");
911
912 std::shared_ptr<const RooFitResult> bestFr;
913 for (auto c : allCats) {
914 allCats.setIndex(c.second);
915 Info("minimize", "Minimizing with discrete %s", c.first.c_str());
916 auto fr = minimize(nll, _fitConfig, nllOpts);
917 if (!fr) {
918 Warning("minimize", "Minimization with discrete %s failed", c.first.c_str());
919 continue;
920 }
921 if (!bestFr || fr->minNll() < bestFr->minNll()) {
922 bestFr = fr;
923 }
924 }
925
926 floatCats.setAttribAll("Constant", false);
927
928 if (!bestFr)
929 return out;
930
931 // create a copy of the fit result, give it a new uuid, and move the const categories into the float area
932 out = std::make_shared<RooFitResult>(*bestFr);
933 const_cast<RooArgList &>(out->floatParsFinal())
934 .addClone(*std::unique_ptr<RooAbsCollection>(out->constPars().selectCommon(floatCats)));
935 const_cast<RooArgList &>(out->floatParsInit()).addClone(*_snap);
936 const_cast<RooArgList &>(out->constPars()).remove(floatCats);
937 out->SetName(TUUID().AsString());
938 }
939
940 bool restore = !fitConfig.UpdateAfterFit();
941 std::string logs;
942 if (!out) {
943 int strategy = fitConfig.MinimizerOptions().Strategy();
944 // Note: AsymptoticCalculator enforces not less than 1 on tolerance - should we do so too?
945 if (_progress) {
946 _nll = new ProgressMonitor(*_nll, _progress);
948 }
949 auto logger = (logSize > 0) ? std::make_unique<cout_redirect>(logs, logSize) : nullptr;
950 RooMinimizer _minimizer(*_nll);
951 _minimizer.fitter()->Config() = fitConfig;
952 // if(fitConfig.MinimizerOptions().ExtraOptions()) {
953 // //for loading hesse options
954 // double a;
955 // if(fitConfig.MinimizerOptions().ExtraOptions()->GetValue("HessianStepTolerance",a)) {
956 // ROOT::Math::MinimizerOptions::Default("Minuit2").SetValue("HessianStepTolerance",a);
957 // }
958 // if(fitConfig.MinimizerOptions().ExtraOptions()->GetValue("HessianG2Tolerance",a)) {
959 // ROOT::Math::MinimizerOptions::Default("Minuit2").SetValue("HessianG2Tolerance",a);
960 // }
961 // }
962
963 bool autoMaxCalls = (_minimizer.fitter()->Config().MinimizerOptions().MaxFunctionCalls() == 0);
964 if (autoMaxCalls) {
966 500 * floatPars->size() * floatPars->size()); // hesse requires O(N^2) function calls
967 }
968 if (_minimizer.fitter()->Config().MinimizerOptions().MaxIterations() == 0) {
969 _minimizer.fitter()->Config().MinimizerOptions().SetMaxIterations(500 * floatPars->size());
970 }
971
972 bool hesse = _minimizer.fitter()->Config().ParabErrors();
973 _minimizer.fitter()->Config().SetParabErrors(
974 false); // turn "off" so can run hesse as a separate step, appearing in status
975 bool minos = _minimizer.fitter()->Config().MinosErrors();
976 _minimizer.fitter()->Config().SetMinosErrors(false);
977 _minimizer.fitter()->Config().SetUpdateAfterFit(true); // note: seems to always take effect
978
979 std::vector<std::pair<std::string, int>> statusHistory;
980
981 // gCurrentSampler = this;
982 // gOldHandlerr = signal(SIGINT,toyInterruptHandlerr);
983
984 TString actualFirstMinimizer = _minimizer.fitter()->Config().MinimizerType();
985
986 int status = 0;
987
988 int constOptimize = 2;
989 _minimizer.fitter()->Config().MinimizerOptions().ExtraOptions()->GetValue("OptimizeConst", constOptimize);
990 if (constOptimize) {
991 _minimizer.optimizeConst(constOptimize);
992 // for safety force a refresh of the cache (and tracking) in the nll
993 // DO NOT do a ConfigChange ... this is just a deactivate-reactivate of caching
994 // but it seems like doing this breaks the const optimization and function is badly behaved
995 // so once its turned on never turn it off.
996 // nll.constOptimizeTestStatistic(RooAbsArg::ConfigChange, constOptimize>1 /* do tracking too if >1 */); //
997 // trigger a re-evaluate of which nodes to cache-and-track
998 // the next line seems safe to do but wont bother doing it because not bothering with above
999 // need to understand why turning the cache off and on again breaks it??
1000 // nll.constOptimizeTestStatistic(RooAbsArg::ValueChange, constOptimize>1); // update the cache values -- is
1001 // this needed??
1002 } else {
1003 // disable const optimization
1004 // warning - if the nll was previously activated then it seems like deactivating may break it.
1005 nll.constOptimizeTestStatistic(RooAbsArg::DeActivate);
1006 }
1007
1008 int sIdx = -1;
1009 TString minim = _minimizer.fitter()->Config().MinimizerType();
1010 TString algo = _minimizer.fitter()->Config().MinimizerAlgoType();
1011 if (minim == "Minuit2") {
1012 if (strategy == -1) {
1013 sIdx = 0;
1014 } else {
1015 sIdx = m_strategy.Index('0' + strategy);
1016 }
1017 if (sIdx == -1) {
1018 Warning("minimize", "Strategy %d not specified in StrategySequence %s ... defaulting to start of sequence",
1019 strategy, m_strategy.Data());
1020 sIdx = 0;
1021 }
1022 } else if (minim == "Minuit")
1023 sIdx = m_strategy.Index('m');
1024
1025 int tries = 0;
1026 int maxtries = 4;
1027 bool first = true;
1028 while (tries < maxtries && sIdx != -1) {
1029 if (m_strategy(sIdx) == 'm') {
1030 minim = "Minuit";
1031 algo = "migradImproved";
1032 } else if (m_strategy(sIdx) == 's') {
1033 algo = "Scan";
1034 } else if (m_strategy(sIdx) == 'h') {
1035 break; // jumping straight to a hesse evaluation
1036 } else {
1037 strategy = int(m_strategy(sIdx) - '0');
1038 _minimizer.setStrategy(strategy);
1039 minim = "Minuit2";
1040 algo = "Migrad";
1041 }
1042 if (auto fff = dynamic_cast<ProgressMonitor *>(_nll); fff) {
1043 fff->fState = minim + algo + std::to_string(_minimizer.fitter()->Config().MinimizerOptions().Strategy());
1044 }
1045 try {
1046 status = _minimizer.minimize(minim, algo);
1047 } catch (const std::exception &e) {
1048 std::cerr << "Exception while minimizing: " << e.what() << std::endl;
1049 }
1050 if (first && actualFirstMinimizer != _minimizer.fitter()->Config().MinimizerType())
1051 actualFirstMinimizer = _minimizer.fitter()->Config().MinimizerType();
1052 first = false;
1053 tries++;
1054
1055 if (auto fff = dynamic_cast<ProgressMonitor *>(_nll); fff && fff->fInterrupt) {
1056 delete _nll;
1057 throw std::runtime_error("Keyboard interrupt while minimizing");
1058 }
1059
1060 // RooMinimizer loses the useful status code, so here we will override it
1061 status = _minimizer.fitter()
1062 ->Result()
1063 .Status(); // note: Minuit failure is status code 4, minuit2 that is edm above max
1064 minim = _minimizer.fitter()->Config().MinimizerType(); // may have changed value
1065 statusHistory.emplace_back(_minimizer.fitter()->Config().MinimizerType() +
1066 _minimizer.fitter()->Config().MinimizerAlgoType() +
1067 std::to_string(_minimizer.fitter()->Config().MinimizerOptions().Strategy()),
1068 status);
1069 if (status % 1000 == 0)
1070 break; // fit was good
1071
1072 if (status == 4 && minim != "Minuit") {
1073 if (printLevel >= -1) {
1074 Warning("fitTo", "%s Hit max function calls of %d", fitName.Data(),
1075 _minimizer.fitter()->Config().MinimizerOptions().MaxFunctionCalls());
1076 }
1077 if (autoMaxCalls) {
1078 if (printLevel >= -1)
1079 Warning("fitTo", "will try doubling this");
1081 _minimizer.fitter()->Config().MinimizerOptions().MaxFunctionCalls() * 2);
1083 _minimizer.fitter()->Config().MinimizerOptions().MaxIterations() * 2);
1084 continue;
1085 }
1086 }
1087
1088 // NOTE: minuit2 seems to distort the tolerance in a weird way, so that tol becomes 1000 times smaller than
1089 // specified Also note that if fits are failing because of edm over max, it can be a good idea to activate the
1090 // Offset option when building nll
1091 if (printLevel >= -1) {
1092 Warning("fitTo", "%s %s%s Status=%d (edm=%f, tol=%f, strat=%d), tries=#%d...", fitName.Data(),
1093 _minimizer.fitter()->Config().MinimizerType().c_str(),
1094 _minimizer.fitter()->Config().MinimizerAlgoType().c_str(), status,
1095 _minimizer.fitter()->Result().Edm(), _minimizer.fitter()->Config().MinimizerOptions().Tolerance(),
1096 _minimizer.fitter()->Config().MinimizerOptions().Strategy(), tries);
1097 }
1098
1099 // decide what to do next based on strategy sequence
1100 if (sIdx == m_strategy.Length() - 1) {
1101 break; // done
1102 }
1103
1104 tries--;
1105 sIdx++;
1106 }
1107
1108 /* Minuit2 status codes:
1109 * status = 0 : OK
1110 status = 1 : Covariance was made pos defined
1111 status = 2 : Hesse is invalid
1112 status = 3 : Edm is above max
1113 status = 4 : Reached call limit
1114 status = 5 : Any other failure
1115
1116 For Minuit its basically 0 is OK, 4 is failure, I think?
1117 */
1118
1119 if (printLevel >= -1 && status != 0) {
1120 Warning("fitTo", "%s final status is %d", fitName.Data(), status);
1121 }
1122
1123 // currently dont have a way to access the covariance "dcovar" which is a metric from iterative
1124 // covariance method that is used by minuit2 to say if the covariance is accurate or not
1125 // See MinimumError.h: IsAccurate if Dcovar < 0.1
1126 // Note that if strategy>=2 or (strategy=1 and Dcovar>0.05) then hesse will be forced to be run (see
1127 // VariadicMetricBuilder) So only in Strategy=0 can you skip hesse (even if SetParabErrors false).
1128
1129 double dCovar = std::numeric_limits<double>::quiet_NaN();
1130 // if(auto _minuit2 = dynamic_cast<ROOT::Minuit2::Minuit2Minimizer*>(_minimizer.fitter()->GetMinimizer());
1131 // _minuit2 && _minuit2->fMinimum) {
1132 // dCovar = _minuit2->fMinimum->Error().Dcovar();
1133 // }
1134
1135 // only do hesse if was a valid min and not full accurate cov matrix already (can happen if e.g. ran strat2)
1136 if (hesse &&
1137 (m_strategy(sIdx) == 'h' || ((strategy < 2 || _minimizer.fitter()->GetMinimizer()->CovMatrixStatus() != 3) &&
1138 _minimizer.fitter()->Result().IsValid()))) {
1139
1140 // Note: minima where the covariance was made posdef are deemed 'valid' ...
1141
1142 // remove limits on pars before calculation - CURRENTLY HAS NO EFFECT, minuit still holds the state as
1143 // transformed interesting note: error on pars before hesse can be significantly smaller than after hesse ...
1144 // what is the pre-hesse error corresponding to? - corresponds to approximation of covariance matrix calculated
1145 // with iterative method
1146 /*auto parSettings = _minimizer.fitter()->Config().ParamsSettings();
1147 for (auto &ss : _minimizer.fitter()->Config().ParamsSettings()) {
1148 ss.RemoveLimits();
1149 }
1150
1151 for(auto f : *floatPars) {
1152 auto v = dynamic_cast<RooRealVar*>(f);
1153 if(v->hasRange(nullptr)) v->setRange("backup",v->getMin(),v->getMax());
1154 v->removeRange();
1155 }*/
1156
1157 // std::cout << "nIterations = " << _minimizer.fitter()->GetMinimizer()->NIterations() << std::endl;
1158 // std::cout << "covQual before hesse = " << _minimizer.fitter()->GetMinimizer()->CovMatrixStatus() <<
1159 // std::endl;
1160 sIdx = -1;
1161 if (hesseStrategy == -1) {
1162 sIdx = 0;
1163 } else {
1164 sIdx = m_hessestrategy.Index('0' + hesseStrategy);
1165 }
1166 if (sIdx == -1) {
1167 Warning("minimize",
1168 "HesseStrategy %d not specified in HesseStrategySequence %s ... defaulting to start of sequence",
1169 hesseStrategy, m_hessestrategy.Data());
1170 sIdx = 0;
1171 }
1172 while (sIdx != -1) {
1173 hesseStrategy = int(m_hessestrategy(sIdx) - '0');
1174 _minimizer.fitter()->Config().MinimizerOptions().SetStrategy(hesseStrategy);
1175 // const_cast<ROOT::Math::IOptions*>(_minimizer.fitter()->Config().MinimizerOptions().ExtraOptions())->SetValue("HessianStepTolerance",0.1);
1176 // const_cast<ROOT::Math::IOptions*>(_minimizer.fitter()->Config().MinimizerOptions().ExtraOptions())->SetValue("HessianG2Tolerance",0.02);
1177
1178 if (auto fff = dynamic_cast<ProgressMonitor *>(_nll); fff) {
1179 fff->fState = TString::Format("Hesse%d", _minimizer.fitter()->Config().MinimizerOptions().Strategy());
1180 fff->counter2 = fff->counter;
1181 }
1182
1183 //_nll->getVal(); // for reasons I dont understand, if nll evaluated before hesse call the edm is smaller? -
1184 // and also becomes WRONG :-S
1185
1186 // auto _status = (_minimizer.fitter()->CalculateHessErrors()) ? _minimizer.fitter()->Result().Status() :
1187 // -1;
1188 auto _status = _minimizer.hesse(); // note: I have seen that you can get 'full covariance quality' without
1189 // running hesse ... is that expected?
1190 // note: hesse status will be -1 if hesse failed (no covariance matrix)
1191 // otherwise the status appears to be whatever was the status before
1192 // note that hesse succeeds even if the cov matrix it calculates is forced pos def. Failure is only
1193 // if it cannot calculate a cov matrix at all.
1194 if (_status != -1)
1195 _status = 0; // mark as hesse succeeded, although need to look at covQual to see if was any good
1196
1197 /*for(auto f : *floatPars) {
1198 auto v = dynamic_cast<RooRealVar*>(f);
1199 if(v->hasRange("backup")) {
1200 v->setRange(v->getMin(),v->getMax());
1201 v->removeRange("backup");
1202 }
1203 }
1204 _minimizer.fitter()->Config().SetParamsSettings(parSettings);*/
1205
1206 /*for (auto &ss : _minimizer.fitter()->Config().ParamsSettings()) {
1207 if( ss.HasLowerLimit() || ss.HasUpperLimit() ) std::cout << ss.Name() << " limit restored " <<
1208 ss.LowerLimit() << " - " << ss.UpperLimit() << std::endl;
1209 }*/
1210
1211 statusHistory.push_back(std::pair<std::string, int>(
1212 TString::Format("Hesse%d", _minimizer.fitter()->Config().MinimizerOptions().Strategy()), _status));
1213
1214 if (auto fff = dynamic_cast<ProgressMonitor *>(_nll); fff && fff->fInterrupt) {
1215 delete _nll;
1216 throw std::runtime_error("Keyboard interrupt while hesse calculating");
1217 }
1218 if ((_status != 0 || _minimizer.fitter()->GetMinimizer()->CovMatrixStatus() != 3) && status == 0 &&
1219 printLevel >= -1) {
1220 Warning("fitTo", "%s hesse status is %d, covQual=%d", fitName.Data(), _status,
1221 _minimizer.fitter()->GetMinimizer()->CovMatrixStatus());
1222 }
1223
1224 if (sIdx >= m_hessestrategy.Length() - 1) {
1225 break; // run out of strategies to try, stop
1226 }
1227
1228 if (_status == 0 && _minimizer.fitter()->GetMinimizer()->CovMatrixStatus() == 3) {
1229 // covariance is valid!
1230 break;
1231 } else if (_status == 0) {
1232 // set the statusHistory to the cov status, since that's more informative
1233 statusHistory.back().second = _minimizer.fitter()->GetMinimizer()->CovMatrixStatus();
1234 }
1235 sIdx++;
1236 } // end of hesse attempt loop
1237 }
1238
1239 // call minos if requested on any parameters
1240 if (status == 0 && minos) {
1241 if (std::unique_ptr<RooAbsCollection> mpars(floatPars->selectByAttrib("minos", true)); !mpars->empty()) {
1242 if (auto fff = dynamic_cast<ProgressMonitor *>(_nll); fff) {
1243 fff->fState = "Minos";
1244 }
1245 auto _status = _minimizer.minos(*mpars);
1246 statusHistory.push_back(std::pair("Minos", _status));
1247 }
1248 }
1249
1250 // DO NOT DO THIS - seems to mess with the NLL function in a way that breaks the cache - reactivating wont fix
1251 // if(constOptimize) { _minimizer.optimizeConst(0); } // doing this because saw happens in RooAbsPdf::minimizeNLL
1252 // method
1253
1254 // signal(SIGINT,gOldHandlerr);
1255 out = std::unique_ptr<RooFitResult>{_minimizer.save(fitName, resultTitle)};
1256
1257 // if status is 0 (min succeeded) but the covQual isn't fully accurate but requested hesse, reflect that in the
1258 // status
1259 if (out->status() == 0 && out->covQual() != 3 && hesse) {
1260 if (out->covQual() == 2) { // was made posdef
1261 out->setStatus(1); // indicates covariance made pos-def
1262 } else { // anything else indicates either hessian is approximate or something else wrong (e.g. not pos-def
1263 // return from strat3)
1264 out->setStatus(2); // hesse invalid
1265 }
1266 }
1267
1268 out->setStatusHistory(statusHistory);
1269
1270 // userPars wont have been added to the RooFitResult by RooMinimizer
1271 const_cast<RooArgList &>(out->constPars()).addClone(fUserPars, true);
1272
1273 if (!std::isnan(dCovar)) {
1274 const_cast<RooArgList &>(out->constPars())
1275 .addClone(RooRealVar(".dCovar", "dCovar from minimization", dCovar), true);
1276 }
1277
1278 if (boundaryCheck) {
1279 // check if any of the parameters are at their limits (potentially a problem with fit)
1280 // or their errors go over their limits (just a warning)
1281 int limit_status = 0;
1282 std::string listpars;
1283 for (auto *v : dynamic_range_cast<RooRealVar *>(*floatPars)) {
1284 if (!v)
1285 continue;
1286 double vRange = v->getMax() - v->getMin();
1287 if (v->getMin() > v->getVal() - vRange * boundaryCheck ||
1288 v->getMax() < v->getVal() + vRange * boundaryCheck) {
1289 // within 0.01% of edge
1290
1291 // check if nll actually lower 'at' the boundary, if it is, refine the best fit to the limit value
1292 auto tmp = v->getVal();
1293 v->setVal(v->getMin());
1294 double boundary_nll = _nll->getVal();
1295 if (boundary_nll <= out->minNll()) {
1296 static_cast<RooRealVar *>(out->floatParsFinal().find(v->GetName()))->setVal(v->getMin());
1297 out->setMinNLL(boundary_nll);
1298 // Info("fit","Corrected %s onto minimum @ %g",v->GetName(),v->getMin());
1299 } else {
1300 // not better, so restore value
1301 v->setVal(tmp);
1302 }
1303
1304 // if has a 'physical' range specified, don't warn if near the limit
1305 if (v->hasRange("physical"))
1306 limit_status = 900;
1307 listpars += v->GetName();
1308 listpars += ",";
1309 } else if (hesse &&
1310 (v->getMin() > v->getVal() - v->getError() || v->getMax() < v->getVal() + v->getError())) {
1311 if (printLevel >= 0) {
1312 Info("minimize", "PARLIM: %s (%f +/- %f) range (%f - %f)", v->GetName(), v->getVal(), v->getError(),
1313 v->getMin(), v->getMax());
1314 }
1315 limit_status = 9000;
1316 }
1317 }
1318 if (limit_status == 900) {
1319 if (printLevel >= 0) {
1320 Warning("minimize", "BOUNDCHK: Parameters within %g%% limit in fit result: %s", boundaryCheck * 100,
1321 listpars.c_str());
1322 }
1323 } else if (limit_status > 0) {
1324 if (printLevel >= 0)
1325 Warning("minimize", "BOUNDCHK: Parameters near limit in fit result");
1326 }
1327
1328 // store the limit check result
1329 statusHistory.emplace_back("BOUNDCHK", limit_status);
1330 out->setStatusHistory(statusHistory);
1331 out->setStatus(out->status() + limit_status);
1332 }
1333
1334 // // automatic parameter range adjustment based on errors
1335 // for(auto a : *floatPars) {
1336 // RooRealVar *v = dynamic_cast<RooRealVar *>(a);
1337 // if(v->getMin() > v->getVal() - 3.*v->getError()) {
1338 // v->setMin(v->getVal() - 3.1*v->getError());
1339 // }
1340 // if(v->getMax() < v->getVal() + 3.*v->getError()) {
1341 // v->setMax(v->getVal() + 3.1*v->getError());
1342 // }
1343 // // also make sure the range isn't too big (fits can struggle)
1344 // if(v->getMin() < v->getVal() - 10.*v->getError()) {
1345 // v->setMin(v->getVal() - 9.9*v->getError());
1346 // }
1347 // if(v->getMax() > v->getVal() + 10.*v->getError()) {
1348 // v->setMax(v->getVal() + 9.9*v->getError());
1349 // }
1350 // }
1351
1352 if (printLevel < 0)
1354
1355 // before returning we will override _minLL with the actual NLL value ... offsetting could have messed up the
1356 // value
1357 out->setMinNLL(_nll->getVal());
1358
1359 // ensure no asymm errors on any pars unless had minuitMinos
1360 for (auto o : out->floatParsFinal()) {
1361 if (auto v = dynamic_cast<RooRealVar *>(o);
1362 v && !v->getAttribute("minos") && !v->getAttribute("xminos") && !v->getAttribute("xMinos"))
1363 v->removeAsymError();
1364 }
1365
1366 // minimizer may have slightly altered the fitConfig (e.g. unavailable minimizer etc) so update for that ...
1367 if (fitConfig.MinimizerOptions().MinimizerType() != actualFirstMinimizer) {
1368 fitConfig.MinimizerOptions().SetMinimizerType(actualFirstMinimizer);
1369 }
1370
1371 if (_progress) {
1372 delete _nll;
1373 }
1374
1375 // call minos if requested on any parameters
1376 if (status == 0 && minos) {
1377 for (auto label : {"xminos", "xMinos"}) {
1378 std::unique_ptr<RooAbsCollection> pars(floatPars->selectByAttrib(label, true));
1379 for (auto p : *pars) {
1380 Info("minimize", "Computing xminos error for %s", p->GetName());
1381 xRooFit::minos(nll, *out, p->GetName(), myFitConfig);
1382 }
1383 if (!pars->empty())
1384 *floatPars = out->floatParsFinal(); // put values back to best fit
1385 }
1386 }
1387 }
1388
1389 if (restore) {
1390 *floatPars = out->floatParsInit();
1391 }
1392
1393 if (out && !logs.empty()) {
1394 // save logs to StringVar in constPars list
1395#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 28, 00)
1396 const_cast<RooArgList &>(out->constPars()).addOwned(std::make_unique<RooStringVar>(".log", "log", logs.c_str()));
1397#else
1398 const_cast<RooArgList &>(out->constPars()).addOwned(*new RooStringVar(".log", "log", logs.c_str()));
1399#endif
1400 }
1401
1402 if (out && cacheDir && cacheDir->IsWritable()) {
1403 // std::cout << "Saving " << out->GetName() << " " << out->GetTitle() << " to " << nll.GetName() << std::endl;
1404 // save a copy of fit result to relevant dir
1405 if (!cacheDir->GetDirectory(nll.GetName()))
1406 cacheDir->mkdir(nll.GetName());
1407 if (auto dir = cacheDir->GetDirectory(nll.GetName()); dir) {
1408 // save NLL opts if was given one, unless already present
1409 if (nllOpts) {
1410 if (strlen(nllOpts->GetName()) == 0) {
1411 nllOpts->SetName(TUUID().AsString());
1412 }
1413 if (!dir->FindKey(nllOpts->GetName())) {
1414 dir->WriteObject(nllOpts.get(), nllOpts->GetName());
1415 }
1416 }
1417
1418 // also save the fitConfig ... unless one with same name already present
1419 std::string configName;
1420 if (!fitConfig.MinimizerOptions().ExtraOptions()->GetValue("Name", configName)) {
1421 auto extraOpts = const_cast<ROOT::Math::IOptions *>(fitConfig.MinimizerOptions().ExtraOptions());
1422 configName = TUUID().AsString();
1423 extraOpts->SetValue("Name", configName.data());
1424 }
1425 if (!dir->GetKey(configName.data())) {
1426 dir->WriteObject(&fitConfig, configName.data());
1427 }
1428 // add the fitConfig name into the fit result before writing, so can retrieve in future
1429#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 28, 00)
1430 const_cast<RooArgList &>(out->constPars())
1431 .addOwned(std::make_unique<RooStringVar>(".fitConfigName", "fitConfigName", configName.c_str()));
1432#else
1433 const_cast<RooArgList &>(out->constPars())
1434 .addOwned(*new RooStringVar(".fitConfigName", "fitConfigName", configName.c_str()));
1435#endif
1436 dir->WriteObject(out.get(), out->GetName());
1437 auto sfr = new StoredFitResult(out);
1438 dir->Add(sfr);
1439 return sfr->fr;
1440 // return std::shared_ptr<const RooFitResult>(out, [](const RooFitResult*){}); // disowned shared_ptr
1441 }
1442 }
1443
1444 return out;
1445}
1446
1447// calculate asymmetric errors, if required, on the named parameter that was floating in the fit
1448// returns status code. 0 = all good, 1 = failure, ...
1449int xRooFit::minos(RooAbsReal &nll, const RooFitResult &ufit, const char *parName,
1450 const std::shared_ptr<ROOT::Fit::FitConfig> &_fitConfig)
1451{
1452
1453 auto par = dynamic_cast<RooRealVar *>(std::unique_ptr<RooArgSet>(nll.getVariables())->find(parName));
1454 if (!par)
1455 return 1;
1456
1457 auto par_hat = dynamic_cast<RooRealVar *>(ufit.floatParsFinal().find(parName));
1458 if (!par_hat)
1459 return 1;
1460
1461 auto myFitConfig = _fitConfig ? _fitConfig : createFitConfig();
1462 auto &fitConfig = *myFitConfig;
1463
1464 bool pErrs = fitConfig.ParabErrors();
1465 fitConfig.SetParabErrors(false);
1466 double mErrs = fitConfig.MinosErrors();
1467 fitConfig.SetMinosErrors(false);
1468
1469 double val_best = par_hat->getVal();
1470 double val_err = (par_hat->hasError() ? par_hat->getError() : -1);
1471 double orig_err = val_err;
1472 double nll_min = ufit.minNll();
1473
1474 int status = 0;
1475
1476 bool isConst = par->isConstant();
1477 par->setConstant(true);
1478
1479 auto findValue = [&](double val_guess, double N_sigma = 1, double precision = 0.002, int printLevel = 0) {
1480 double tmu;
1481 int nrItr = 0;
1482 double sigma_guess = std::abs((val_guess - val_best) / N_sigma);
1483 double val_pre =
1484 val_guess -
1485 10 * precision * sigma_guess; // this is just to set value st. guarantees will do at least one iteration
1486 bool lastOverflow = false;
1487 bool lastUnderflow = false;
1488 while (std::abs(val_pre - val_guess) > precision * sigma_guess) {
1489 val_pre = val_guess;
1490 if (val_guess > 0 && par->getMax() < val_guess)
1491 par->setMax(2 * val_guess);
1492 if (val_guess < 0 && par->getMin() > val_guess)
1493 par->setMin(2 * val_guess);
1494 par->setVal(val_guess);
1495 // std::cout << "Guessing " << val_guess << std::endl;
1496 auto result = xRooFit::minimize(nll, myFitConfig);
1497 if (!result) {
1498 status = 1;
1499 return std::numeric_limits<double>::quiet_NaN();
1500 }
1501 double nll_val = result->minNll();
1502 status += result->status() * 10;
1503 tmu = 2 * (nll_val - nll_min);
1504 sigma_guess = std::abs(val_guess - val_best) / sqrt(tmu);
1505
1506 if (tmu <= 0) {
1507 // found an alternative or improved minima
1508 std::cout << "Warning: Alternative best-fit of " << par->GetName() << " @ " << val_guess << " vs "
1509 << val_best << " (delta=" << tmu / 2. << ")" << std::endl;
1510 double new_guess = val_guess + (val_guess - val_best);
1511 val_best = val_guess;
1512 val_guess = new_guess;
1513 sigma_guess = std::abs((val_guess - val_best) / N_sigma);
1514 val_pre = val_guess - 10 * precision * sigma_guess;
1515 status = (status / 10) * 10 + 1;
1516 continue;
1517 }
1518
1519 double corr = /*damping_factor**/ (val_pre - val_best - N_sigma * sigma_guess);
1520
1521 // subtract off the difference in the new and damped correction
1522 val_guess -= corr;
1523
1524 if (printLevel > 1) {
1525 // cout << "nPars: " << nPars << std::endl;
1526 // cout << "NLL: " << nll->GetName() << " = " << nll->getVal() << endl;
1527 // cout << "delta(NLL): " << nll->getVal()-nll_min << endl;
1528 std::cout << "NLL min: " << nll_min << std::endl;
1529 std::cout << "N_sigma*sigma(pre): " << std::abs(val_pre - val_best) << std::endl;
1530 std::cout << "sigma(guess): " << sigma_guess << std::endl;
1531 std::cout << "par(guess): " << val_guess + corr << std::endl;
1532 std::cout << "true val: " << val_best << std::endl;
1533 std::cout << "tmu: " << tmu << std::endl;
1534 std::cout << "Precision: " << sigma_guess * precision << std::endl;
1535 std::cout << "Correction: " << (-corr < 0 ? " " : "") << -corr << std::endl;
1536 std::cout << "N_sigma*sigma(guess): " << std::abs(val_guess - val_best) << std::endl;
1537 std::cout << std::endl;
1538 }
1539 if (val_guess > par->getMax()) {
1540 if (lastOverflow) {
1541 val_guess = par->getMin();
1542 break;
1543 }
1544 lastOverflow = true;
1545 lastUnderflow = false;
1546 val_guess = par->getMax() - 1e-12;
1547 } else if (val_guess < par->getMin()) {
1548 if (lastUnderflow) {
1549 val_guess = par->getMin();
1550 break;
1551 }
1552 lastOverflow = false;
1553 lastUnderflow = true;
1554 val_guess = par->getMin() + 1e-12;
1555 } else {
1556 lastUnderflow = false;
1557 lastOverflow = false;
1558 }
1559
1560 nrItr++;
1561 if (nrItr > 25) {
1562 status = (status / 10) * 10 + 3;
1563 break;
1564 }
1565 }
1566
1567 if (lastOverflow) {
1568 // msg().Error("findSigma","%s at upper limit of %g .. error may be underestimated
1569 // (t=%g)",par->GetName(),par->getMax(),tmu);
1570 status = (status / 10) * 10 + 2;
1571 } else if (lastUnderflow) {
1572 // msg().Error("findSigma","%s at lower limit of %g .. error may be underestimated
1573 // (t=%g)",par->GetName(),par->getMin(),tmu);
1574 status = (status / 10) * 10 + 2;
1575 }
1576
1577 if (printLevel > 1)
1578 std::cout << "Found sigma for nll " << nll.GetName() << ": " << (val_guess - val_best) / N_sigma << std::endl;
1579 if (printLevel > 1)
1580 std::cout << "Finished in " << nrItr << " iterations." << std::endl;
1581 if (printLevel > 1)
1582 std::cout << std::endl;
1583 return (val_guess - val_best) / N_sigma;
1584 };
1585
1586 // determine if asym error defined by temporarily setting error to nan ... will then return non-nan if defined
1587
1588 par_hat->setError(std::numeric_limits<double>::quiet_NaN());
1589 double lo = par_hat->getErrorLo();
1590 double hi = par_hat->getErrorHi();
1591 if (std::isnan(hi)) {
1592 hi = findValue(val_best + val_err, 1) + val_best -
1593 par_hat->getVal(); // put error wrt par_hat value, even if found better min
1594 if (hi > val_err)
1595 val_err = hi; // in case val_err was severe underestimate, don't want to waste time being too 'near' min
1596 }
1597 if (std::isnan(lo)) {
1598 lo = -findValue(val_best - val_err, -1) + val_best -
1599 par_hat->getVal(); // put error wrt par_hat value, even if found better min
1600 }
1601 dynamic_cast<RooRealVar *>(ufit.floatParsFinal().find(parName))->setAsymError(lo, hi);
1602 par_hat->setError(orig_err);
1603
1604 fitConfig.SetParabErrors(pErrs);
1605 fitConfig.SetMinosErrors(mErrs);
1606 par->setConstant(isConst);
1607
1608 std::vector<std::pair<std::string, int>> statusHistory;
1609 for (unsigned int i = 0; i < ufit.numStatusHistory(); i++) {
1610 statusHistory.emplace_back(ufit.statusLabelHistory(i), ufit.statusCodeHistory(i));
1611 }
1612 statusHistory.emplace_back(TString::Format("xMinos:%s", parName), status);
1613 const_cast<RooFitResult &>(ufit).setStatusHistory(statusHistory);
1614 const_cast<RooFitResult &>(ufit).setStatus(ufit.status() + status);
1615
1616 return status;
1617}
1618
1619TCanvas *
1620xRooFit::hypoTest(RooWorkspace &w, int nToysNull, int /*nToysAlt*/, const xRooFit::Asymptotics::PLLType &pllType)
1621{
1622 TCanvas *out = nullptr;
1623
1624 // 1. Determine pdf: use top-level, if more than 1 then exit and tell user they need to flag
1625 RooAbsPdf *model = nullptr;
1626 std::deque<RooAbsArg *> topPdfs;
1627 int flagCount = 0;
1628 for (auto p : w.allPdfs()) {
1629 if (p->hasClients())
1630 continue;
1631 flagCount += p->getAttribute("hypoTest");
1632 if (p->getAttribute("hypoTest")) {
1633 topPdfs.push_front(p);
1634 } else {
1635 topPdfs.push_back(p);
1636 }
1637 }
1638 if (topPdfs.empty()) {
1639 Error("hypoTest", "Cannot find top-level pdf in workspace");
1640 return nullptr;
1641 } else if (topPdfs.size() > 1) {
1642 // should be one flagged
1643 if (flagCount == 0) {
1644 Error("hypoTest", "Multiple top-level pdfs. Flag which one to test with "
1645 "w->pdf(\"pdfName\")->setAttribute(\"hypoTest\",true)");
1646 return out;
1647 } else if (flagCount != 1) {
1648 Error("hypoTest", "Multiple top-level pdfs flagged for hypoTest -- pick one.");
1649 return out;
1650 }
1651 }
1652 model = dynamic_cast<RooAbsPdf *>(topPdfs.front());
1653
1654 Info("hypoTest", "Using PDF: %s", model->GetName());
1655
1656 double CL = 0.95; // TODO: make configurable
1657
1658 // 2. Determine the data (including globs). if more than 1 then exit and tell user they need to flag
1659 RooAbsData *obsData = nullptr;
1660 std::shared_ptr<RooArgSet> obsGlobs = nullptr;
1661
1662 for (auto p : w.allData()) {
1663 if (obsData) {
1664 Error("hypoTest", "Multiple datasets in workspace. Flag which one to test with "
1665 "w->data(\"dataName\")->setAttribute(\"hypoTest\",true)");
1666 return out;
1667 }
1668 obsData = p;
1669 }
1670
1671 if (!obsData) {
1672 Error("hypoTest", "No data -- cannot determine observables");
1673 return nullptr;
1674 }
1675
1676 Info("hypoTest", "Using Dataset: %s", obsData->GetName());
1677
1678 {
1679 auto _globs = xRooNode(w).datasets()[obsData->GetName()]->globs(); // keep alive because may own the globs
1680 obsGlobs = std::make_shared<RooArgSet>();
1681 obsGlobs->addClone(_globs.argList());
1682 Info("hypoTest", "Using Globs: %s", (obsGlobs->empty()) ? " <NONE>" : obsGlobs->contentsString().c_str());
1683 }
1684
1685 // 3. Determine the POI and args - look for model pars with "hypoPoints" binning, if none then cannot scan
1686 // args are const, poi are floating - exception is if only one then assume it is the POI
1687 auto _vars = std::unique_ptr<RooArgSet>(model->getVariables());
1688 RooArgSet poi;
1689 RooArgSet args;
1690 for (auto _v : *_vars) {
1691 if (auto v = dynamic_cast<RooRealVar *>(_v); v && v->hasBinning("hypoPoints")) {
1692 poi.add(*v);
1693 }
1694 }
1695 if (poi.size() > 1) {
1696 auto _const = std::unique_ptr<RooAbsCollection>(poi.selectByAttrib("Constant", true));
1697 args.add(*_const);
1698 poi.remove(*_const);
1699 }
1700 if (!args.empty()) {
1701 Info("hypoTest", "Using Arguments: %s", args.contentsString().c_str());
1702 }
1703 if (poi.empty()) {
1704 Error("hypoTest", "No POI detected: add the hypoPoints binning to at least one non-const model parameter e.g.:\n "
1705 "w->var(\"mu\")->setBinning(RooUniformBinning(0.5,10.5,10),\"hypoPoints\"))");
1706 return nullptr;
1707 }
1708
1709 Info("hypoTest", "Using Parameters of Interest: %s", poi.contentsString().c_str());
1710
1712
1713 // should check if exist in workspace
1714 auto nllOpts = createNLLOptions();
1715 auto fitConfig = createFitConfig();
1716
1717 xRooNLLVar nll(*model, std::make_pair(obsData, obsGlobs.get()), *nllOpts);
1718 nll.SetFitConfig(fitConfig);
1719
1720 if (poi.size() == 1) {
1721 auto mu = dynamic_cast<RooRealVar *>(poi.first());
1722
1723 double altVal = (mu->getStringAttribute("altVal")) ? TString(mu->getStringAttribute("altVal")).Atof()
1724 : std::numeric_limits<double>::quiet_NaN();
1725
1726 if (std::isnan(altVal) && mu->hasRange("physical")) {
1727 // use the smallest absolute value for the altValue
1728 altVal = mu->getMin("physical");
1729 Info("hypoTest", "No altVal specified - using min of given physical range = %g", altVal);
1730 } else {
1731 if (!std::isnan(altVal)) {
1732 Info("hypoTest", "alt hypo: %g - CLs activated", altVal);
1733 } else {
1734 Info("hypoTest", "No altVal found - to specify setStringAttribute(\"altVal\",\"<value>\") on POI or set "
1735 "the physical range");
1736 }
1737 }
1738 bool doCLs = !std::isnan(altVal) && std::abs(mu->getMin("hypoPoints")) > altVal &&
1739 std::abs(mu->getMax("hypoPoints")) > altVal;
1740
1741 const char *sCL = (doCLs) ? "CLs" : "null";
1742 Info("hypoTest", "%s testing active", sCL);
1743
1744 auto obs_ts = new TGraphErrors;
1745 obs_ts->SetNameTitle("obs_ts", TString::Format("Observed TestStat;%s", mu->GetTitle()));
1746 auto obs_pcls = new TGraphErrors;
1747 obs_pcls->SetNameTitle(TString::Format("obs_p%s", sCL),
1748 TString::Format("Observed p_{%s};%s", sCL, mu->GetTitle()));
1749 auto obs_cls = new TGraphErrors;
1750 obs_cls->SetNameTitle(TString::Format("obs_%s", sCL), TString::Format("Observed %s;%s", sCL, mu->GetTitle()));
1751
1752 std::vector<int> expSig = {-2, -1, 0, 1, 2};
1753 if (std::isnan(altVal))
1754 expSig.clear();
1755 std::map<int, TGraphErrors> exp_pcls;
1756 std::map<int, TGraphErrors> exp_cls;
1757 for (auto &s : expSig) {
1758 exp_pcls[s].SetNameTitle(TString::Format("exp%d_p%s", s, sCL),
1759 TString::Format("Expected (%d#sigma) p_{%s};%s", s, sCL, mu->GetTitle()));
1760 exp_cls[s].SetNameTitle(TString::Format("exp%d_%s", s, sCL),
1761 TString::Format("Expected (%d#sigma) %s;%s", s, sCL, mu->GetTitle()));
1762 }
1763
1764 auto getLimit = [CL](TGraphErrors &pValues) {
1765 double _out = std::numeric_limits<double>::quiet_NaN();
1766 bool lastAbove = false;
1767 for (int i = 0; i < pValues.GetN(); i++) {
1768 bool thisAbove = pValues.GetPointY(i) >= (1. - CL);
1769 if (i != 0 && thisAbove != lastAbove) {
1770 // crossed over ... find limit by interpolation
1771 // using linear interpolation so far
1772 _out = pValues.GetPointX(i - 1) + (pValues.GetPointX(i) - pValues.GetPointX(i - 1)) *
1773 ((1. - CL) - pValues.GetPointY(i - 1)) /
1774 (pValues.GetPointY(i) - pValues.GetPointY(i - 1));
1775 }
1776 lastAbove = thisAbove;
1777 }
1778 return _out;
1779 };
1780
1781 auto testPoint = [&](double testVal) {
1782 auto hp = nll.hypoPoint(mu->GetName(), testVal, altVal, pllType);
1783 obs_ts->SetPoint(obs_ts->GetN(), testVal, hp.pll().first);
1784 obs_ts->SetPointError(obs_ts->GetN() - 1, 0, hp.pll().second);
1785
1786 if (nToysNull > 0) {
1787 }
1788
1789 obs_pcls->SetPoint(obs_pcls->GetN(), testVal, (doCLs) ? hp.pCLs_asymp().first : hp.pNull_asymp().first);
1790 obs_pcls->SetPointError(obs_pcls->GetN() - 1, 0, (doCLs) ? hp.pCLs_asymp().second : hp.pNull_asymp().second);
1791 for (auto &s : expSig) {
1792 exp_pcls[s].SetPoint(exp_pcls[s].GetN(), testVal,
1793 (doCLs) ? hp.pCLs_asymp(s).first : hp.pNull_asymp(s).first);
1794 }
1795 if (doCLs) {
1796 Info("hypoTest", "%s=%g: %s=%g sigma_mu=%g %s=%g", mu->GetName(), testVal, obs_ts->GetName(),
1797 obs_ts->GetPointY(obs_ts->GetN() - 1), hp.sigma_mu().first, obs_pcls->GetName(),
1798 obs_pcls->GetPointY(obs_pcls->GetN() - 1));
1799 } else {
1800 Info("hypoTest", "%s=%g: %s=%g %s=%g", mu->GetName(), testVal, obs_ts->GetName(),
1801 obs_ts->GetPointY(obs_ts->GetN() - 1), obs_pcls->GetName(), obs_pcls->GetPointY(obs_pcls->GetN() - 1));
1802 }
1803 };
1804
1805 if (mu->getBins("hypoPoints") <= 0) {
1806 // autoTesting
1807 // evaluate min and max points
1808 testPoint(mu->getMin("hypoPoints"));
1809 testPoint(mu->getMax("hypoPoints"));
1810 testPoint((mu->getMax("hypoPoints") + mu->getMin("hypoPoints")) / 2.);
1811
1812 while (std::abs(obs_pcls->GetPointY(obs_pcls->GetN() - 1) - (1. - CL)) > 0.01) {
1813 obs_pcls->Sort();
1814 double nextTest = getLimit(*obs_pcls);
1815 if (std::isnan(nextTest))
1816 break;
1817 testPoint(nextTest);
1818 }
1819 for (auto s : expSig) {
1820 while (std::abs(exp_pcls[s].GetPointY(exp_pcls[s].GetN() - 1) - (1. - CL)) > 0.01) {
1821 exp_pcls[s].Sort();
1822 double nextTest = getLimit(exp_pcls[s]);
1823 if (std::isnan(nextTest))
1824 break;
1825 testPoint(nextTest);
1826 }
1827 }
1828 obs_ts->Sort();
1829 obs_pcls->Sort();
1830 for (auto &s : expSig)
1831 exp_pcls[s].Sort();
1832
1833 } else {
1834 for (int i = 0; i <= mu->getBins("hypoPoints"); i++) {
1835 testPoint((i == mu->getBins("hypoPoints")) ? mu->getBinning("hypoPoints").binHigh(i - 1)
1836 : mu->getBinning("hypoPoints").binLow(i));
1837 }
1838 }
1839
1840 obs_cls->SetPoint(obs_cls->GetN(), getLimit(*obs_pcls), 0.05);
1841 for (auto &s : expSig) {
1842 exp_cls[s].SetPoint(exp_cls[s].GetN(), getLimit(exp_pcls[s]), 0.05);
1843 }
1844
1845 // if more than two hypoPoints, visualize as bands
1846 if (exp_pcls[2].GetN() > 1) {
1847 TGraph *band2 = new TGraph;
1848 band2->SetNameTitle(".pCLs_2sigma", "2 sigma band");
1849 TGraph *band2up = new TGraph;
1850 band2up->SetNameTitle(".pCLs_2sigma_upUncert", "");
1851 TGraph *band2down = new TGraph;
1852 band2down->SetNameTitle(".pCLs_2sigma_downUncert", "");
1853 band2->SetFillColor(kYellow);
1854 band2up->SetFillColor(kYellow);
1855 band2down->SetFillColor(kYellow);
1856 band2up->SetFillStyle(3005);
1857 band2down->SetFillStyle(3005);
1858 for (int i = 0; i < exp_pcls[2].GetN(); i++) {
1859 band2->SetPoint(band2->GetN(), exp_pcls[2].GetPointX(i),
1860 exp_pcls[2].GetPointY(i) - exp_pcls[2].GetErrorYlow(i));
1861 band2up->SetPoint(band2up->GetN(), exp_pcls[2].GetPointX(i),
1862 exp_pcls[2].GetPointY(i) + exp_pcls[2].GetErrorYhigh(i));
1863 }
1864 for (int i = exp_pcls[2].GetN() - 1; i >= 0; i--) {
1865 band2up->SetPoint(band2up->GetN(), exp_pcls[2].GetPointX(i),
1866 exp_pcls[2].GetPointY(i) - exp_pcls[2].GetErrorYlow(i));
1867 }
1868 for (int i = 0; i < exp_pcls[-2].GetN(); i++) {
1869 band2down->SetPoint(band2down->GetN(), exp_pcls[-2].GetPointX(i),
1870 exp_pcls[-2].GetPointY(i) + exp_pcls[-2].GetErrorYhigh(i));
1871 }
1872 for (int i = exp_pcls[-2].GetN() - 1; i >= 0; i--) {
1873 band2->SetPoint(band2->GetN(), exp_pcls[-2].GetPointX(i),
1874 exp_pcls[-2].GetPointY(i) + exp_pcls[-2].GetErrorYhigh(i));
1875 band2down->SetPoint(band2down->GetN(), exp_pcls[-2].GetPointX(i),
1876 exp_pcls[-2].GetPointY(i) - exp_pcls[-2].GetErrorYlow(i));
1877 }
1878 band2->SetBit(kCanDelete);
1879 band2up->SetBit(kCanDelete);
1880 band2down->SetBit(kCanDelete);
1881 auto ax = static_cast<TNamed *>(band2->Clone(".axis"));
1882 ax->SetTitle(TString::Format("Hypothesis Test;%s", mu->GetTitle()));
1883 ax->Draw("AF");
1884 band2->Draw("F");
1885 band2up->Draw("F");
1886 band2down->Draw("F");
1887 }
1888
1889 if (exp_pcls[1].GetN() > 1) {
1890 TGraph *band2 = new TGraph;
1891 band2->SetNameTitle(".pCLs_1sigma", "1 sigma band");
1892 TGraph *band2up = new TGraph;
1893 band2up->SetNameTitle(".pCLs_1sigma_upUncert", "");
1894 TGraph *band2down = new TGraph;
1895 band2down->SetNameTitle(".pCLs_1sigma_downUncert", "");
1896 band2->SetFillColor(kGreen);
1897 band2up->SetFillColor(kGreen);
1898 band2down->SetFillColor(kGreen);
1899 band2up->SetFillStyle(3005);
1900 band2down->SetFillStyle(3005);
1901 for (int i = 0; i < exp_pcls[1].GetN(); i++) {
1902 band2->SetPoint(band2->GetN(), exp_pcls[1].GetPointX(i),
1903 exp_pcls[1].GetPointY(i) - exp_pcls[1].GetErrorYlow(i));
1904 band2up->SetPoint(band2up->GetN(), exp_pcls[1].GetPointX(i),
1905 exp_pcls[1].GetPointY(i) + exp_pcls[1].GetErrorYhigh(i));
1906 }
1907 for (int i = exp_pcls[1].GetN() - 1; i >= 0; i--) {
1908 band2up->SetPoint(band2up->GetN(), exp_pcls[1].GetPointX(i),
1909 exp_pcls[1].GetPointY(i) - exp_pcls[1].GetErrorYlow(i));
1910 }
1911 for (int i = 0; i < exp_pcls[-1].GetN(); i++) {
1912 band2down->SetPoint(band2down->GetN(), exp_pcls[-1].GetPointX(i),
1913 exp_pcls[-1].GetPointY(i) + exp_pcls[-1].GetErrorYhigh(i));
1914 }
1915 for (int i = exp_pcls[-1].GetN() - 1; i >= 0; i--) {
1916 band2->SetPoint(band2->GetN(), exp_pcls[-1].GetPointX(i),
1917 exp_pcls[-1].GetPointY(i) + exp_pcls[-1].GetErrorYhigh(i));
1918 band2down->SetPoint(band2down->GetN(), exp_pcls[-1].GetPointX(i),
1919 exp_pcls[-1].GetPointY(i) - exp_pcls[-1].GetErrorYlow(i));
1920 }
1921 band2->SetBit(kCanDelete);
1922 band2up->SetBit(kCanDelete);
1923 band2down->SetBit(kCanDelete);
1924 band2->Draw("F");
1925 band2up->Draw("F");
1926 band2down->Draw("F");
1927 }
1928
1929 TObject *expPlot = nullptr;
1930 if (exp_cls[0].GetN() > 0) {
1931 exp_pcls[0].SetLineStyle(2);
1932 exp_pcls[0].SetFillColor(kGreen);
1933 exp_pcls[0].SetMarkerStyle(0);
1934 expPlot = exp_pcls[0].DrawClone("L");
1935 }
1936 obs_pcls->SetBit(kCanDelete);
1937 obs_pcls->Draw(gPad->GetListOfPrimitives()->IsEmpty() ? "ALP" : "LP");
1938
1939 obs_ts->SetLineColor(kRed);
1940 obs_ts->SetMarkerColor(kRed);
1941 obs_ts->SetBit(kCanDelete);
1942 obs_ts->Draw("LP");
1943
1944 auto l = new TLegend(0.5, 0.6, 1. - gPad->GetRightMargin(), 1. - gPad->GetTopMargin());
1945 l->SetName("legend");
1946 l->AddEntry(obs_ts, obs_ts->GetTitle(), "LPE");
1947 l->AddEntry(obs_pcls, obs_pcls->GetTitle(), "LPE");
1948 if (expPlot)
1949 l->AddEntry(expPlot, "Expected", "LFE");
1951 l->Draw();
1952
1953 obs_cls->SetMarkerStyle(29);
1954 obs_cls->SetEditable(false);
1955 obs_cls->Draw("LP");
1956 for (auto s : expSig) {
1957 exp_cls[s].SetMarkerStyle(29);
1958 exp_cls[s].SetEditable(false);
1959 exp_cls[s].DrawClone("LP");
1960 }
1961 }
1962
1963 if (out)
1964 out->RedrawAxis();
1965
1966 return out;
1967}
1968
1969double round_to_digits(double value, int digits)
1970{
1971 if (value == 0.0)
1972 return 0.0;
1973 double factor = pow(10.0, digits - ceil(log10(std::abs(value))));
1974 return std::round(value * factor) / factor;
1975}
1976double round_to_decimal(double value, int decimal_places)
1977{
1978 const double multiplier = std::pow(10.0, decimal_places);
1979 return std::round(value * multiplier) / multiplier;
1980}
1981
1982// rounds error to 1 or 2 sig fig and round value to match that precision
1983std::pair<double, double> xRooFit::matchPrecision(const std::pair<double, double> &in)
1984{
1985 auto out = in;
1986 if (!std::isinf(out.second)) {
1987 auto tmp = out.second;
1988 out.second = round_to_digits(out.second, 2);
1989 int expo = (out.second == 0) ? 0 : (int)std::floor(std::log10(std::abs(out.second)));
1990 if (TString::Format("%e", out.second)(0) != '1') {
1991 out.second = round_to_digits(tmp, 1);
1992 out.first = (expo >= 0) ? round(out.first) : round_to_decimal(out.first, -expo);
1993 } else if (out.second != 0) {
1994 out.first = (expo >= 0) ? round(out.first) : round_to_decimal(out.first, -expo + 1);
1995 }
1996 }
1997 return out;
1998}
1999
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define e(i)
Definition RSha256.hxx:103
#define oocoutW(o, a)
#define oocoutE(o, a)
#define oocoutI(o, a)
@ kRed
Definition Rtypes.h:66
@ kGreen
Definition Rtypes.h:66
@ kYellow
Definition Rtypes.h:66
static void indent(ostringstream &buf, int indent_level)
#define gDirectory
Definition TDirectory.h:384
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition TError.cxx:218
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:229
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
char name[80]
Definition TGX11.cxx:110
#define hi
@ kCanDelete
Definition TObject.h:367
#define gROOT
Definition TROOT.h:407
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
R__EXTERN TSystem * gSystem
Definition TSystem.h:555
#define gPad
double getValV(const RooArgSet *) const override
Return value of object.
Definition xRooFit.cxx:588
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition xRooFit.cxx:599
std::string fState
Definition xRooFit.cxx:693
static bool fInterrupt
Definition xRooFit.cxx:544
bool setData(RooAbsData &data, bool cloneData) override
Definition xRooFit.cxx:587
ProgressMonitor(const ProgressMonitor &other, const char *name=nullptr)
Definition xRooFit.cxx:575
~ProgressMonitor() override
Definition xRooFit.cxx:567
TStopwatch s
Definition xRooFit.cxx:705
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Interface for detailed printing of object.
Definition xRooFit.cxx:590
TObject * clone(const char *newname) const override
Definition xRooFit.cxx:579
std::shared_ptr< RooAbsCollection > vars
Definition xRooFit.cxx:706
bool getParameters(const RooArgSet *observables, RooArgSet &outputSet, bool stripDisconnected) const override
Fills a list with leaf nodes in the arg tree starting with ourself as top node that don't match any o...
Definition xRooFit.cxx:583
RooArgList minPars
Definition xRooFit.cxx:701
RooRealProxy fFunc
Definition xRooFit.cxx:698
void applyWeightSquared(bool flag) override
Disables or enables the usage of squared weights.
Definition xRooFit.cxx:589
void constOptimizeTestStatistic(ConstOpCode opcode, bool doAlsoTrackingOpt) override
Interface function signaling a request to perform constant term optimization.
Definition xRooFit.cxx:594
ProgressMonitor(RooAbsReal &f, int interval=30)
Definition xRooFit.cxx:556
static ProgressMonitor * me
Definition xRooFit.cxx:543
double defaultErrorLevel() const override
Definition xRooFit.cxx:582
RooArgList prevPars
Definition xRooFit.cxx:702
static void interruptHandler(int signum)
Definition xRooFit.cxx:545
void(* oldHandlerr)(int)
Definition xRooFit.cxx:542
static int minos(RooAbsReal &nll, const RooFitResult &ufit, const char *parName="", const std::shared_ptr< ROOT::Fit::FitConfig > &_fitConfig=nullptr)
Definition xRooFit.cxx:1449
static std::shared_ptr< const RooFitResult > minimize(RooAbsReal &nll, const std::shared_ptr< ROOT::Fit::FitConfig > &fitConfig=nullptr, const std::shared_ptr< RooLinkedList > &nllOpts=nullptr)
Definition xRooFit.cxx:718
static std::shared_ptr< RooLinkedList > createNLLOptions()
Definition xRooFit.cxx:456
static TCanvas * hypoTest(RooWorkspace &w, const xRooFit::Asymptotics::PLLType &pllType=xRooFit::Asymptotics::Unknown)
Definition xRooFit.h:226
static std::shared_ptr< ROOT::Fit::FitConfig > createFitConfig()
Definition xRooFit.cxx:483
static std::pair< double, double > matchPrecision(const std::pair< double, double > &in)
Definition xRooFit.cxx:1983
This xRooNLLVar object has several special methods, e.g.
Definition xRooNLLVar.h:59
xRooFitResult minimize(const std::shared_ptr< ROOT::Fit::FitConfig > &=nullptr)
Class describing the configuration of the fit, options and parameter settings using the ROOT::Fit::Pa...
Definition FitConfig.h:47
void SetMinosErrors(bool on=true)
set Minos errors computation to be performed after fitting
Definition FitConfig.h:229
const std::string & MinimizerAlgoType() const
return type of minimizer algorithms
Definition FitConfig.h:192
bool ParabErrors() const
do analysis for parabolic errors
Definition FitConfig.h:205
void SetUpdateAfterFit(bool on=true)
Update configuration after a fit using the FitResult.
Definition FitConfig.h:243
void SetParabErrors(bool on=true)
set parabolic errors
Definition FitConfig.h:226
const std::string & MinimizerType() const
return type of minimizer package
Definition FitConfig.h:187
ROOT::Math::MinimizerOptions & MinimizerOptions()
access to the minimizer control parameter (non const method)
Definition FitConfig.h:167
bool MinosErrors() const
do minos errors analysis on the parameters
Definition FitConfig.h:208
bool IsValid() const
True if fit successful, otherwise false.
Definition FitResult.h:105
double Edm() const
Expected distance from minimum.
Definition FitResult.h:117
int Status() const
minimizer status code
Definition FitResult.h:128
ROOT::Math::Minimizer * GetMinimizer() const
return pointer to last used minimizer (is NULL in case fit is not yet done) This pointer is guarantee...
Definition Fitter.h:444
const FitResult & Result() const
get fit result
Definition Fitter.h:394
const FitConfig & Config() const
access to the fit configuration (const method)
Definition Fitter.h:422
class implementing generic options for a numerical algorithm Just store the options in a map of strin...
Generic interface for defining configuration options of a numerical algorithm.
Definition IOptions.h:28
void SetValue(const char *name, double val)
generic methods for retrieving options
Definition IOptions.h:42
bool GetValue(const char *name, T &t) const
Definition IOptions.h:54
void SetMaxFunctionCalls(unsigned int maxfcn)
set maximum of function calls
void SetStrategy(int stra)
set the strategy
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
const IOptions * ExtraOptions() const
return extra options (NULL pointer if they are not present)
double Tolerance() const
absolute tolerance
unsigned int MaxIterations() const
max iterations
unsigned int MaxFunctionCalls() const
max number of function calls
virtual int CovMatrixStatus() const
return status of covariance matrix using Minuit convention {0 not calculated 1 approximated 2 made po...
Definition Minimizer.h:271
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
bool dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=nullptr, bool valueOnly=false) const
Test whether we depend on (ie, are served by) any object in the specified collection.
RooWorkspace * _myws
Prevent 'AlwaysDirty' mode for this node.
Definition RooAbsArg.h:732
bool isConstant() const
Check if the "Constant" attribute is set.
Definition RooAbsArg.h:361
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
RooWorkspace * workspace() const
Definition RooAbsArg.h:550
virtual void constOptimizeTestStatistic(ConstOpCode opcode, bool doAlsoTrackingOpt=true)
Interface function signaling a request to perform constant term optimization.
bool getAttribute(const Text_t *name) const
Check if a named attribute is set. By default, all attributes are unset.
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expression tree)
virtual void applyWeightSquared(bool flag)
Disables or enables the usage of squared weights.
RooAbsArg * findServer(const char *name) const
Return server of this with name name. Returns nullptr if not found.
Definition RooAbsArg.h:210
A space to attach TBranches.
RooAbsCollection * selectByAttrib(const char *name, bool value) const
Create a subset of the current collection, consisting only of those elements with the specified attri...
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
Storage_t const & get() const
Const access to the underlying stl container.
void assignFast(const RooAbsCollection &other, bool setValDirty=true) const
Functional equivalent of assign() but assumes this and other collection have same layout.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void setAttribAll(const Text_t *name, bool value=true)
Set given attribute in each element of the collection by calling each elements setAttribute() functio...
Storage_t::size_type size() const
RooAbsArg * first() const
virtual RooAbsArg * addClone(const RooAbsArg &var, bool silent=false)
Add a clone of the specified argument to list.
std::string contentsString() const
Return comma separated list of contained object names as STL string.
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
virtual double defaultErrorLevel() const
Definition RooAbsReal.h:248
virtual bool setData(RooAbsData &, bool=true)
Definition RooAbsReal.h:374
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Structure printing.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:178
static TClass * Class()
Implements a RooAbsBinning in terms of an array of boundary values, posing no constraints on the choi...
Definition RooBinning.h:27
Named container for two doubles, two integers two object points and three string pointers that can be...
Definition RooCmdArg.h:26
Container class to hold unbinned data.
Definition RooDataSet.h:57
RooRealVar * weightVar() const
Returns a pointer to the weight variable (if set).
Definition RooDataSet.h:129
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Int_t statusCodeHistory(UInt_t icycle) const
const char * statusLabelHistory(UInt_t icycle) const
const RooArgList & floatParsFinal() const
Return list of floating parameters after fit.
Int_t status() const
Return MINUIT status code.
UInt_t numStatusHistory() const
double minNll() const
Return minimized -log(L) value.
static TClass * Class()
static TClass * Class()
Collection class for internal use, storing a collection of RooAbsArg pointers in a doubly linked list...
static TClass * Class()
Wrapper class around ROOT::Fit:Fitter that provides a seamless interface between the minimizer functi...
void optimizeConst(int flag)
If flag is true, perform constant term optimization on function being minimized.
RooFit::OwningPtr< RooFitResult > save(const char *name=nullptr, const char *title=nullptr)
Save and return a RooFitResult snapshot of current minimizer status.
int minos()
Execute MINOS.
int hesse()
Execute HESSE.
int minimize(const char *type, const char *alg=nullptr)
Minimise the function passed in the constructor.
ROOT::Fit::Fitter * fitter()
Return underlying ROOT fitter object.
void setStrategy(int istrat)
Change MINUIT strategy to istrat.
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
Implementation of a RooCacheManager<RooAbsCacheElement> that specializes in the storage of cache elem...
Poisson pdf.
Definition RooPoisson.h:19
static TClass * Class()
void setNoRounding(bool flag=true)
Switch off/on rounding of x to the nearest integer.
Definition RooPoisson.h:36
Efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:33
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:48
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
void setError(double value)
Definition RooRealVar.h:60
bool hasBinning(const char *name) const override
Returns true if variable has a binning named 'name'.
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
static RooAbsData * GenerateAsimovData(const RooAbsPdf &pdf, const RooArgSet &observables)
generate the asimov data for the observables (not the global ones) need to deal with the case of a si...
A RooAbsArg implementing string values.
Joins several RooAbsCategoryLValue objects into a single category.
bool setIndex(value_type index, bool printError=true) override
Set the value of the super category to the specified index.
Implementation of RooAbsBinning that provides a uniform binning in 'n' bins between the range end poi...
Persistable container for RooFit projects.
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:37
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition TAttFill.h:39
Using a TBrowser one can browse all ROOT objects.
Definition TBrowser.h:37
void SetStatusText(const char *txt, Int_t col)
Definition TBrowser.h:112
The Canvas class.
Definition TCanvas.h:23
static TCanvas * MakeDefCanvas()
Static function to build a default canvas.
Definition TCanvas.cxx:1516
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition TClass.h:81
TClass * IsA() const override
Definition TClass.h:618
static TClass * GetClass(const char *name, Bool_t load=kTRUE, Bool_t silent=kFALSE)
Static method returning pointer to TClass of the specified class name.
Definition TClass.cxx:2968
This class stores the date and time with a precision of one second in an unsigned 32 bit word (950130...
Definition TDatime.h:37
const char * AsString() const
Return the date & time as a string (ctime() format).
Definition TDatime.cxx:102
Describe directory structure in memory.
Definition TDirectory.h:45
virtual TDirectory * GetDirectory(const char *namecycle, Bool_t printError=false, const char *funcname="GetDirectory")
Find a directory using apath.
virtual Bool_t IsWritable() const
Definition TDirectory.h:237
virtual TDirectory * mkdir(const char *name, const char *title="", Bool_t returnExistingDirectory=kFALSE)
Create a sub-directory "a" or a hierarchy of sub-directories "a/b/c/...".
void SetName(const char *newname) override
Set the name for directory If the directory name is changed after the directory was written once,...
A TGraphErrors is a TGraph with error bars.
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2325
Int_t GetN() const
Definition TGraph.h:131
void Draw(Option_t *chopt="") override
Draw this graph with its current attributes.
Definition TGraph.cxx:814
void SetNameTitle(const char *name="", const char *title="") override
Set graph name and title.
Definition TGraph.cxx:2400
Book space in a file, create I/O buffers, to fill them, (un)compress them.
Definition TKey.h:28
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
TObject * Clone(const char *newname="") const override
Make a clone of an object using the Streamer facility.
Definition TNamed.cxx:74
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Mother of all ROOT objects.
Definition TObject.h:41
virtual const char * GetName() const
Returns name of object.
Definition TObject.cxx:438
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad with: gROOT->SetSelectedPad(c1).
Definition TObject.cxx:299
virtual void Delete(Option_t *option="")
Delete this object.
Definition TObject.cxx:248
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:780
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition TObject.cxx:524
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition TRandom.cxx:615
virtual UInt_t Integer(UInt_t imax)
Returns a random integer uniformly distributed on the interval [ 0, imax-1 ].
Definition TRandom.cxx:361
Stopwatch class.
Definition TStopwatch.h:28
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void Continue()
Resume a stopped stopwatch.
void Reset()
Definition TStopwatch.h:52
Provides iteration through tokens of a given string.
Definition TPRegexp.h:143
Bool_t NextToken()
Get the next token, it is stored in this TString.
Basic string class.
Definition TString.h:139
Ssiz_t Length() const
Definition TString.h:417
Double_t Atof() const
Return floating-point value contained in string.
Definition TString.cxx:2054
Bool_t IsFloat() const
Returns kTRUE if string contains a floating point or integer number.
Definition TString.cxx:1858
const char * Data() const
Definition TString.h:376
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2378
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition TString.h:651
virtual Bool_t ProcessEvents()
Process pending events (GUI, timers, sockets).
Definition TSystem.cxx:416
This class defines a UUID (Universally Unique IDentifier), also known as GUIDs (Globally Unique IDent...
Definition TUUID.h:42
TDatime GetTime() const
Get time from UUID.
Definition TUUID.cxx:670
const char * AsString() const
Return UUID as string. Copy string immediately since it will be reused.
Definition TUUID.cxx:571
RooCmdArg WeightVar(const char *name="weight", bool reinterpretAsWeight=false)
RooCmdArg Offset(std::string const &mode)
RooCmdArg Optimize(Int_t flag=2)
RooCmdArg Extended(bool flag=true)
RooCmdArg ExpectedData(bool flag=true)
Double_t x[n]
Definition legend1.C:17
double nll(double pdf, double weight, int binnedL, int doBinOffset)
Definition MathFuncs.h:321
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Sort the n elements of the array a of generic templated type Element.
Definition TMathBase.h:431
#define BEGIN_XROOFIT_NAMESPACE
Definition Config.h:24
#define END_XROOFIT_NAMESPACE
Definition Config.h:25
TLine l
Definition textangle.C:4
#define GIT_COMMIT_DATE
#define GIT_COMMIT_HASH
double round_to_decimal(double value, int decimal_places)
Definition xRooFit.cxx:1976
double round_to_digits(double value, int digits)
Definition xRooFit.cxx:1969