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