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