Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
xRooHypoSpace.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
14
15#include "RooArgSet.h"
16#include "RooArgList.h"
17#include "RooRealVar.h"
18#include "RooFitResult.h"
19#include "RooConstVar.h"
20#include "RooAbsPdf.h"
21
22#include "TCanvas.h"
23#include "TStopwatch.h"
24#include "TSystem.h"
25#include "TPRegexp.h"
26#include "TMemFile.h"
27#include "TROOT.h"
28#include "RooDataSet.h"
29#include "TKey.h"
30#include "TFile.h"
31#include "TGraphErrors.h"
32#include "TMultiGraph.h"
33#include "TH1F.h"
34#include "TStyle.h"
35#include "TLegend.h"
36#include "TLine.h"
38#include "TEnv.h"
39
41
42xRooNLLVar::xRooHypoSpace::xRooHypoSpace(const char *name, const char *title)
43 : TNamed(name, title), fPars(std::make_shared<RooArgSet>())
44{
45 if (name == nullptr || strlen(name) == 0) {
46 SetName(TUUID().AsString());
47 }
48}
49
51 : fPars(std::make_shared<RooArgSet>())
52{
53 if (!result)
54 return;
55
57
58 fPars->addClone(*std::unique_ptr<RooAbsCollection>(result->GetParameters()));
59 double spaceSize = 1;
60 for (auto p : *fPars) {
61 auto v = dynamic_cast<RooRealVar *>(p);
62 if (!v)
63 continue;
64 spaceSize *= (v->getMax() - v->getMin());
65 }
66 for (int i = 0; i < result->ArraySize(); i++) {
67 auto point = result->GetResult(i);
68 double xVal = result->GetXValue(i);
69 double ssize = spaceSize;
70 for (auto p : *fPars) {
71 auto v = dynamic_cast<RooRealVar *>(p);
72 if (!v)
73 continue;
74 ssize /= (v->getMax() - v->getMin());
75 double remain = std::fmod(xVal, ssize);
76 v->setVal((xVal - remain) / ssize);
77 xVal = remain;
78 }
79 emplace_back(xRooHypoPoint(std::make_shared<RooStats::HypoTestResult>(*point), fPars.get()));
80 }
81 // add any pars we might have missed
82 for (auto &p : *this) {
83 for (auto a : *p.coords) {
84 if (!fPars->find(a->GetName()))
85 fPars->addClone(*a);
86 }
87 }
88}
89
90std::shared_ptr<xRooNode> xRooNLLVar::xRooHypoSpace::pdf(const char *parValues) const
91{
92 return pdf(toArgs(parValues));
93}
94
95std::shared_ptr<xRooNode> xRooNLLVar::xRooHypoSpace::pdf(const RooAbsCollection &parValues) const
96{
97 RooArgList rhs;
98 rhs.add(parValues);
99 rhs.sort();
100
101 std::shared_ptr<xRooNode> out = nullptr;
102
103 for (auto &[_range, _pdf] : fPdfs) {
104 // any pars not in rhs are assumed to have infinite range in rhs
105 // and vice versa
106 bool collision = true;
107 for (auto &_lhs : *_range) {
108 auto _rhs = rhs.find(*_lhs);
109 if (!_rhs)
110 continue;
111 if (auto v = dynamic_cast<RooRealVar *>(_rhs); v) {
112 if (auto v2 = dynamic_cast<RooRealVar *>(_lhs)) {
113 if (!(v->getMin() <= v2->getMax() && v2->getMin() <= v->getMax())) {
114 collision = false;
115 break;
116 }
117 } else if (auto c2 = dynamic_cast<RooConstVar *>(_lhs)) {
118 if (!(v->getMin() <= c2->getVal() && c2->getVal() <= v->getMax())) {
119 collision = false;
120 break;
121 }
122 }
123 } else if (auto c = dynamic_cast<RooConstVar *>(_rhs); c) {
124 if (auto v2 = dynamic_cast<RooRealVar *>(_lhs)) {
125 if (!(c->getVal() <= v2->getMax() && v2->getMin() <= c->getVal())) {
126 collision = false;
127 break;
128 }
129 } else if (auto c2 = dynamic_cast<RooConstVar *>(_lhs)) {
130 if (!(c->getVal() == c2->getVal())) {
131 collision = false;
132 break;
133 }
134 }
135 }
136 }
137 if (collision) {
138 if (out) {
139 throw std::runtime_error("Multiple pdf possibilities");
140 }
141 out = _pdf;
142 }
143 }
144
145 return out;
146}
147
149{
150
151 RooArgList out;
152
153 TStringToken pattern(str, ";");
154 while (pattern.NextToken()) {
155 TString s = pattern;
156 // split by "=" sign
157 auto _idx = s.Index('=');
158 if (_idx == -1)
159 continue;
160 TString _name = s(0, _idx);
161 TString _val = s(_idx + 1, s.Length());
162
163 if (_val.IsFloat()) {
164 out.addClone(RooConstVar(_name, _name, _val.Atof()));
165 } else if (_val.BeginsWith('[')) {
166 _idx = _val.Index(',');
167 if (_idx == -1)
168 continue;
169 TString _min = _val(0, _idx);
170 TString _max = _val(_idx + 1, _val.Length() - _idx - 2);
171 out.addClone(RooRealVar(_name, _name, _min.Atof(), _max.Atof()));
172 }
173 }
174
175 return out;
176}
177
178int xRooNLLVar::xRooHypoSpace::AddPoints(const char *parName, size_t nPoints, double low, double high)
179{
180 if (nPoints == 0)
181 return nPoints;
182
183 auto _par = dynamic_cast<RooAbsRealLValue *>(fPars->find(parName));
184 if (!_par)
185 throw std::runtime_error("Unknown parameter");
186 _par->setAttribute("axis");
187
188 if (nPoints == 1) {
189 _par->setVal((high + low) * 0.5);
190 AddPoint();
191 return nPoints;
192 }
193
194 double step = (high - low) / (nPoints - 1);
195 if (step <= 0)
196 throw std::runtime_error("Invalid steps");
197
198 for (size_t i = 0; i < nPoints; i++) {
199 _par->setVal((i == nPoints - 1) ? high : (low + step * i));
200 AddPoint();
201 }
202 return nPoints;
203}
204
206{
207 if (axes().empty()) {
208 // set the first poi as the axis variable to scan
209 if (poi().empty()) {
210 throw std::runtime_error("No POI to scan");
211 } else {
212 poi().first()->setAttribute("axis");
213 }
214 }
215
216 if (empty()) {
217 // promote all axes to being poi and demote all non-axes to non-poi
218 poi().setAttribAll("poi", false);
219 axes().setAttribAll("poi");
220 }
221
222 return AddPoint(TString::Format("%s=%f", axes().first()->GetName(), value));
223}
224
225int xRooNLLVar::xRooHypoSpace::scan(const char *type, size_t nPoints, double low, double high,
226 const std::vector<double> &nSigmas, double relUncert)
227{
228
229 TString sType(type);
230 sType.ToLower();
231 if (sType.Contains("cls") && !sType.Contains("pcls"))
232 sType.ReplaceAll("cls", "pcls");
233 if (!sType.Contains("pcls") && !sType.Contains("ts") && !sType.Contains("pnull") && !sType.Contains("plr")) {
234 throw std::runtime_error("scan type must be equal to one of: plr, cls, ts, pnull");
235 }
236
237 // will scan the first axes variable ... if there is none, specify the first poi as the axis var
238 if (axes().empty()) {
239 // set the first poi as the axis variable to scan
240 if (poi().empty()) {
241 throw std::runtime_error("No POI to scan");
242 } else {
243 poi().first()->setAttribute("axis");
244 }
245 }
246
247 // promote all axes to being poi and demote all non-axes to non-poi
248 poi().setAttribAll("poi", false);
249 axes().setAttribAll("poi");
250
251 auto p = dynamic_cast<RooRealVar *>(axes().first());
252 if (!p) {
253 throw std::runtime_error(TString::Format("%s not scannable", axes().first()->GetName()));
254 }
255
256 if (sType.Contains("cls")) {
257 if (empty() && relUncert == std::numeric_limits<double>::infinity()) {
258 // use default uncertainty precision of 10%
259 ::Info("xRooHypoSpace::scan", "Using default precision of 10%% for auto-scan");
260 relUncert = 0.1;
261 }
262 for (auto a : axes()) {
263 if (!a->hasRange("physical")) {
264 ::Info("xRooHypoSpace::scan", "No physical range set for %s, setting to [0,inf]", p->GetName());
265 dynamic_cast<RooRealVar *>(a)->setRange("physical", 0, std::numeric_limits<double>::infinity());
266 }
267 if (!a->getStringAttribute("altVal") || !strlen(p->getStringAttribute("altVal"))) {
268 ::Info("xRooHypoSpace::scan", "No altVal set for %s, setting to 0", a->GetName());
269 a->setStringAttribute("altVal", "0");
270 }
271 // ensure range straddles altVal
272 double altVal = TString(a->getStringAttribute("altVal")).Atof();
273 auto v = dynamic_cast<RooRealVar *>(a);
274 if (v->getMin() >= altVal) {
275 ::Info("xRooHypoSpace::scan", "range of POI does not straddle alt value, adjusting minimum to %g",
276 altVal - 1e-5);
277 v->setMin(altVal - 1e-5);
278 }
279 if (v->getMax() <= altVal) {
280 ::Info("xRooHypoSpace::scan", "range of POI does not straddle alt value, adjusting maximum to %g",
281 altVal + 1e-5);
282 v->setMax(altVal + 1e-5);
283 }
284 for (auto &[pdf, nll] : fNlls) {
285 if (auto _v = dynamic_cast<RooRealVar *>(nll->pars()->find(*a))) {
286 _v->setRange(v->getMin(), v->getMax());
287 }
288 }
289 }
290 } else if (sType.Contains("plr")) {
291 // force use of two-sided test statistic for any new points
292 fTestStatType = xRooFit::Asymptotics::TwoSided;
293 sType.ReplaceAll("plr", "ts");
294 }
295
296 if (p && high <= low) {
297 // take from parameter
298 low = p->getMin("scan");
299 high = p->getMax("scan");
300 ::Info("xRooHypoSpace::scan", "Using %s range: %g - %g", p->GetName(), low, high);
301 }
302 if (!std::isnan(low) && !std::isnan(high) && !(std::isinf(low) && std::isinf(high))) {
303 p->setRange("scan", low, high);
304 }
305
306 bool doObs = false;
307 for (auto nSigma : nSigmas) {
308 if (std::isnan(nSigma)) {
309 doObs = true;
310 break;
311 }
312 }
313
314 if (fNlls.empty()) {
315 // this happens when loaded hypoSpace from a hypoSpaceInverterResult
316 // set relUncert to infinity so that we don't test any new points
317 relUncert = std::numeric_limits<double>::infinity(); // no NLL available so just get whatever limit we can
318
319 // if any of the defined points are 'expected' data don't do obs
320 // for(auto& hp : *this) {
321 // if(hp.isExpected) {
322 // doObs = false; break;
323 // }
324 // }
325 } else {
326#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 26, 00)
327// bool allGen = true;
328// for (auto &[pdf, nll]: fNlls) {
329// auto _d = dynamic_cast<RooDataSet *>(nll->data());
330// if (!_d || !_d->weightVar() || !_d->weightVar()->getStringAttribute("fitResult") ||
331// !_d->weightVar()->getAttribute("expected")) {
332// allGen = false;
333// break;
334// }
335// }
336// if (allGen)
337// doObs = false;
338#endif
339 }
340
341 // create a fitDatabase if required
342 TDirectory *origDir = gDirectory;
343 if (!gDirectory || !gDirectory->IsWritable()) {
344 // locate a TMemFile in the open list of files and move to that
345 // or create one if cannot find
346 for (auto file : *gROOT->GetListOfFiles()) {
347 if (auto f = dynamic_cast<TMemFile *>(file)) {
348 f->cd();
349 break;
350 }
351 }
352 if (!gDirectory || !gDirectory->IsWritable()) {
353 new TMemFile("fitDatabase", "RECREATE");
354 }
355 }
356
357 int out = 0;
358
359 if (nPoints == 0) {
360 // automatic scan
361 if (sType.Contains("cls")) {
362 for (double nSigma : nSigmas) {
363 xValueWithError res(std::make_pair(0., 0.));
364 if (std::isnan(nSigma)) {
365 if (!doObs)
366 continue;
367 res = findlimit(TString::Format("%s obs", sType.Data()), relUncert);
368 } else {
369 res =
370 findlimit(TString::Format("%s exp%s%d", sType.Data(), nSigma > 0 ? "+" : "", int(nSigma)), relUncert);
371 }
372 if (std::isnan(res.first) || std::isnan(res.second)) {
373 out = 1;
374 } else if (std::isinf(res.second)) {
375 out = 2;
376 }
377 }
378 } else {
379 throw std::runtime_error(TString::Format("Automatic scanning not yet supported for %s", type));
380 }
381 } else {
382 // add the required points and then compute the required value
383 if (nPoints == 1) {
384 AddPoint(TString::Format("%s=%g", poi().first()->GetName(), (high + low) / 2.));
385 graphs(sType); // triggers computation
386 } else {
387 double step = (high - low) / (nPoints - 1);
388 for (size_t i = 0; i < nPoints; i++) {
389 AddPoint(TString::Format("%s=%g", poi().first()->GetName(), low + step * i));
390 graphs(sType); // triggers computation
391 }
392 }
393 }
394
395 if (origDir)
396 origDir->cd();
397
398 return out;
399}
400
401std::map<std::string, std::pair<double, double>>
402xRooNLLVar::xRooHypoSpace::limits(const char *opt, const std::vector<double> &nSigmas, double relUncert)
403{
404
405 if (fNlls.empty()) {
406 // this happens when loaded hypoSpace from a hypoSpaceInverterResult
407 // set relUncert to infinity so that we don't test any new points
408 relUncert = std::numeric_limits<double>::infinity(); // no NLL available so just get whatever limit we can
409 }
410
411 scan(opt, nSigmas, relUncert);
412
413 std::map<std::string, std::pair<double, double>> out;
414 for (auto nSigma : nSigmas) {
415 auto lim = limit(opt, nSigma);
416 if (lim.second < 0)
417 lim.second = -lim.second; // make errors positive for this method
418 out[std::isnan(nSigma) ? "obs" : TString::Format("%d", int(nSigma)).Data()] = xRooFit::matchPrecision(lim);
419 }
420 return out;
421}
422
424{
425 // move to given coords, if any ... will mark them const too
426 std::unique_ptr<RooAbsCollection, std::function<void(RooAbsCollection *)>> _snap(fPars->snapshot(),
427 [&](RooAbsCollection *c) {
428 *fPars = *c;
429 delete c;
430 });
431 TStringToken pattern(coords, ";");
432 while (pattern.NextToken()) {
433 TString s = pattern;
434 // split by "=" sign
435 auto _idx = s.Index('=');
436 if (_idx == -1)
437 continue;
438 TString _name = s(0, _idx);
439 TString _val = s(_idx + 1, s.Length());
440 auto _v = dynamic_cast<RooRealVar *>(fPars->find(_name));
441 if (!_v)
442 continue;
443
444 if (_val.IsFloat()) {
445 _v->setConstant();
446 _v->setVal(_val.Atof());
447 }
448 }
449
450 auto _pdf = pdf();
451
452 if (!_pdf)
453 throw std::runtime_error("no model at coordinates");
454
455 // if (std::unique_ptr<RooAbsCollection>(fPars->selectByAttrib("poi", true))->size() == 0) {
456 // throw std::runtime_error(
457 // "No pars designated as POI - set with pars()->find(<parName>)->setAttribute(\"poi\",true)");
458 // }
459
460 if (fNlls.find(_pdf) == fNlls.end()) {
461 fNlls[_pdf] = std::make_shared<xRooNLLVar>(_pdf->nll("" /*TODO:allow change dataset name and nll opts*/, {}));
462 }
463
464 xRooHypoPoint out;
465
466 out.nllVar = fNlls[_pdf];
467 out.fData = fNlls[_pdf]->getData();
468 out.isExpected = dynamic_cast<RooDataSet *>(out.fData.first.get()) &&
469 dynamic_cast<RooDataSet *>(out.fData.first.get())->weightVar()->getAttribute("expected");
470 // TODO: need to access the genfit of the data and add that to the point, somehow ...
471
472 out.coords.reset(fPars->snapshot()); // should already have altVal prop on poi, and poi labelled
473 // ensure all poi are marked const ... required by xRooHypoPoint behaviour
474 out.poi().setAttribAll("Constant");
475 // and now remove anything that's marked floating
476
477 // do to bug in remove have to ensure not using the hash map otherwise will be doing an invalid read after the
478 // deletion of the owned pars
479 const_cast<RooAbsCollection *>(out.coords.get())->useHashMapForFind(false);
480 const_cast<RooAbsCollection *>(out.coords.get())
481 ->remove(*std::unique_ptr<RooAbsCollection>(out.coords->selectByAttrib("Constant", false)), true, true);
482
483 double value = out.fNullVal();
484 double alt_value = out.fAltVal();
485
486 auto _type = fTestStatType;
487 if (_type == xRooFit::Asymptotics::Unknown) {
488 // decide based on values
489 if (std::isnan(alt_value)) {
491 } else if (value >= alt_value) {
493 } else {
495 }
496 }
497
498 out.fPllType = _type;
499
500 // look for a matching point
501 for (auto &p : *this) {
502 if (p.nllVar != out.nllVar)
503 continue;
504 if (p.fData != out.fData)
505 continue;
506 if (!p.alt_poi().equals(out.alt_poi()))
507 continue;
508 bool match = true;
509 for (auto c : p.alt_poi()) {
510 if (auto v = dynamic_cast<RooAbsReal *>(c);
511 v && std::abs(v->getVal() - out.alt_poi().getRealValue(v->GetName())) > 1e-12) {
512 match = false;
513 break;
514 }
515 }
516 if (!match)
517 continue;
518 if (!p.coords->equals(*out.coords))
519 continue;
520 for (auto c : *p.coords) {
521 if (c->getAttribute("poi")) {
522 continue; // check poi below
523 }
524 if (auto v = dynamic_cast<RooAbsReal *>(c);
525 v && std::abs(v->getVal() - out.coords->getRealValue(v->GetName())) > 1e-12) {
526 match = false;
527 break;
528 }
529 }
530 if (!match)
531 continue;
532 // if reached here we can copy over the asimov dataset to save re-generating it
533 // first copy over cfit_alt (if its there) because that is also the same since coords and alt_poi values same
534 if (auto cfit = p.cfit_alt(true)) {
535 out.fAlt_cfit = cfit;
536 }
537 if (p.asimov(true) && p.asimov(true)->fData.first && (!out.asimov(true) || !out.asimov(true)->fData.first)) {
538 out.asimov()->fData = p.asimov(true)->fData;
539 }
540 if (!p.poi().equals(out.poi()))
541 continue;
542 for (auto c : p.poi()) {
543 if (auto v = dynamic_cast<RooAbsReal *>(c);
544 v && std::abs(v->getVal() - out.poi().getRealValue(v->GetName())) > 1e-12) {
545 match = false;
546 break;
547 }
548 }
549 if (match) {
550 // found a duplicate point, return that!
551 return p;
552 }
553 }
554
555 std::string coordString;
556 for (auto a : axes()) {
557 coordString += TString::Format("%s=%g", a->GetName(), out.coords->getRealValue(a->GetName()));
558 coordString += ",";
559 }
560 coordString.erase(coordString.end() - 1);
561
562 ::Info("xRooHypoSpace::AddPoint", "Added new point @ %s", coordString.c_str());
563 return emplace_back(out);
564}
565
566bool xRooNLLVar::xRooHypoSpace::AddModel(const xRooNode &_pdf, const char *validity)
567{
568
569 if (!_pdf.get<RooAbsPdf>()) {
570 throw std::runtime_error("Not a pdf");
571 }
572
573 auto pars = _pdf.pars().argList();
574
575 // replace any pars with validity pars and add new pars
576 auto vpars = toArgs(validity);
577 pars.replace(vpars);
578 vpars.remove(pars, true, true);
579 pars.add(vpars);
580
581 if (auto existing = pdf(pars)) {
582 throw std::runtime_error(std::string("Clashing model: ") + existing->GetName());
583 }
584
585 auto myPars = std::shared_ptr<RooArgList>(dynamic_cast<RooArgList *>(pars.snapshot()));
586 myPars->sort();
587
588 pars.remove(*fPars, true, true);
589
590 fPars->addClone(pars);
591
592 fPdfs.insert(std::make_pair(myPars, std::make_shared<xRooNode>(_pdf)));
593
594 return true;
595}
596
598{
599 // determine which pars are the minimal set to distinguish all points in the space
600 RooArgList out;
601 out.setName("axes");
602
603 out.add(*std::unique_ptr<RooAbsCollection>(
604 fPars->selectByAttrib("axis", true))); // start with any pars explicitly designated as axes
605
606 bool clash;
607 do {
608 clash = false;
609
610 std::set<std::vector<double>> coords;
611 for (auto &p : *this) {
612 std::vector<double> p_coords;
613 for (auto o : out) {
614 auto _v = dynamic_cast<RooRealVar *>(p.coords->find(o->GetName()));
615 p_coords.push_back(
616 (_v && _v->isConstant())
617 ? _v->getVal()
618 : std::numeric_limits<double>::infinity()); // non-const coords are treating as non-existent
619 // p_coords.push_back(p.coords->getRealValue(o->GetName(), std::numeric_limits<double>::quiet_NaN()));
620 }
621 if (coords.find(p_coords) != coords.end()) {
622 clash = true;
623 break;
624 }
625 coords.insert(p_coords);
626 }
627
628 if (clash) {
629 // add next best coordinate
630 std::map<std::string, std::unordered_set<double>> values;
631 for (auto &par : *pars()) {
632 if (out.find(*par))
633 continue;
634 for (auto p : *this) {
635 auto _v = dynamic_cast<RooRealVar *>(p.coords->find(par->GetName()));
636 values[par->GetName()].insert(
637 (_v && _v->isConstant())
638 ? _v->getVal()
639 : std::numeric_limits<double>::infinity()); // non-const coords are treating as non-existent
640 // values[par->GetName()].insert(
641 // p.coords->getRealValue(par->GetName(), std::numeric_limits<double>::quiet_NaN()));
642 }
643 }
644
645 std::string bestVar;
646 size_t maxDiff = 0;
647 bool isPOI = false;
648 for (auto &[k, v] : values) {
649 if (v.size() > maxDiff || (v.size() == maxDiff && !isPOI && pars()->find(k.c_str())->getAttribute("poi"))) {
650 bestVar = k;
651 isPOI = pars()->find(k.c_str())->getAttribute("poi");
652 maxDiff = std::max(maxDiff, v.size());
653 }
654 }
655 if (bestVar.empty()) {
656 break;
657 }
658
659 out.add(*pars()->find(bestVar.c_str()));
660 }
661 } while (clash);
662
663 // ensure poi are at the end
664 std::unique_ptr<RooAbsCollection> poi(out.selectByAttrib("poi", true));
665 out.remove(*poi);
666 out.add(*poi);
667
668 return out;
669}
670
672{
673 RooArgList out;
674 out.setName("poi");
675 out.add(*std::unique_ptr<RooAbsCollection>(pars()->selectByAttrib("poi", true)));
676 return out;
677}
678
680{
681
682 if (!gDirectory)
683 return;
684 auto dir = gDirectory->GetDirectory(apath);
685 if (!dir) {
686 // try open file first
687 TString s(apath);
688 auto f = TFile::Open(s.Contains(":") ? TString(s(0, s.Index(":"))) : s);
689 if (f) {
690 if (!s.Contains(":"))
691 s += ":";
692 dir = gDirectory->GetDirectory(s);
693 if (dir) {
694 LoadFits(s);
695 return;
696 }
697 }
698 if (!dir) {
699 Error("LoadFits", "Path not found %s", apath);
700 return;
701 }
702 }
703
704 // assume for now all fits in given path will have the same pars
705 // so can just look at the float and const pars of first fit result to get all of them
706 // tuple is: parName, parValue, parAltValue (blank if nan)
707 // key represents the ufit values, value represents the sets of poi for the available cfits (subfits of the ufit)
708
709 std::map<std::set<std::tuple<std::string, double, std::string>>, std::set<std::set<std::string>>> cfits;
710 std::set<std::string> allpois;
711
712 int nFits = 0;
713 std::function<void(TDirectory *)> processDir;
714 processDir = [&](TDirectory *_dir) {
715 std::cout << "Processing " << _dir->GetName() << std::endl;
716 if (auto keys = _dir->GetListOfKeys(); keys) {
717 // first check if dir doesn't contain any RooLinkedList ... this identifies it as not an nll dir
718 // so treat any sub-dirs as new nll
719 bool isNllDir = false;
720 for (auto &&k : *keys) {
721 TKey *key = dynamic_cast<TKey *>(k);
722 if (strcmp(key->GetClassName(), "RooLinkedList") == 0) {
723 isNllDir = true;
724 break;
725 }
726 }
727
728 for (auto &&k : *keys) {
729 if (auto subdir = _dir->GetDirectory(k->GetName()); subdir) {
730 if (!isNllDir) {
731 LoadFits(subdir->GetPath());
732 } else {
733 processDir(subdir);
734 }
735 continue;
736 }
737 auto cl = TClass::GetClass((static_cast<TKey *>(k))->GetClassName());
738 if (cl->InheritsFrom("RooFitResult")) {
739 if (auto cachedFit = _dir->Get<RooFitResult>(k->GetName()); cachedFit) {
740 nFits++;
741 if (nFits == 1) {
742 // for first fit add any missing float pars
743 std::unique_ptr<RooAbsCollection> snap(cachedFit->floatParsFinal().snapshot());
744 snap->remove(*fPars, true, true);
745 fPars->addClone(*snap);
746 // add also the non-string const pars
747 for (auto &p : cachedFit->constPars()) {
748 if (p->getAttribute("global"))
749 continue; // don't consider globals
750 auto v = dynamic_cast<RooAbsReal *>(p);
751 if (!v) {
752 continue;
753 };
754 if (!fPars->contains(*v))
755 fPars->addClone(*v);
756 }
757 }
758 // get names of all the floats
759 std::set<std::string> floatPars;
760 for (auto &p : cachedFit->floatParsFinal())
761 floatPars.insert(p->GetName());
762 // see if
763
764 // build a set of the const par values
765 std::set<std::tuple<std::string, double, std::string>> constPars;
766 for (auto &p : cachedFit->constPars()) {
767 if (p->getAttribute("global"))
768 continue; // don't consider globals when looking for cfits
769 auto v = dynamic_cast<RooAbsReal *>(p);
770 if (!v) {
771 continue;
772 };
773 constPars.insert(
774 std::make_tuple(v->GetName(), v->getVal(),
775 v->getStringAttribute("altVal") ? v->getStringAttribute("altVal") : ""));
776 }
777 // now see if this is a subset of any existing cfit ...
778 for (auto &&[key, value] : cfits) {
779 if (constPars == key)
780 continue; // ignore cases where we already recorded this list of constPars
781 if (std::includes(constPars.begin(), constPars.end(), key.begin(), key.end())) {
782 // usual case ... cachedFit has more constPars than one of the fits we have already encountered
783 // (the ufit)
784 // => cachedFit is a cfit of key fr ...
785 std::set<std::string> pois;
786 for (auto &&par : constPars) {
787 if (key.find(par) == key.end()) {
788 pois.insert(std::get<0>(par));
789 allpois.insert(std::get<0>(par));
790 }
791 }
792 if (!pois.empty()) {
793 cfits[constPars].insert(pois);
794 // std::cout << cachedFit->GetName() << " ";
795 // for(auto ff: constPars) std::cout << ff.first << "=" <<
796 // ff.second << " "; std::cout << std::endl;
797 }
798 }
799 /* FOR NOW we will skip cases where we encounter the cfit before the ufit - usually should eval the
800 ufit first
801 * else if (std::includes(key.begin(), key.end(), constPars.begin(), constPars.end())) {
802 // constPars are subset of key
803 // => key is a ufit of the cachedFit
804 // add all par names of key that aren't in constPars ... these are the poi
805 std::set<std::string> pois;
806 for (auto &&par: key) {
807 if (constPars.find(par) == constPars.end()) {
808 pois.insert(std::get<0>(par));
809 allpois.insert(std::get<0>(par));
810 }
811 }
812 if (!pois.empty()) {
813 std::cout << "found cfit BEFORE ufit??" << std::endl;
814 value.insert(pois);
815 }
816 } */
817 }
818 // ensure that this combination of constPars has entry in map,
819 // even if it doesn't end up with any poi identified from cfits to it
820 cfits[constPars];
821 delete cachedFit;
822 }
823 }
824 }
825 }
826 };
827 processDir(dir);
828 ::Info("xRooHypoSpace::xRooHypoSpace", "%s - Loaded %d fits", apath, nFits);
829
830 if (allpois.size() == 1) {
831 ::Info("xRooHypoSpace::xRooHypoSpace", "Detected POI: %s", allpois.begin()->c_str());
832
833 auto nll = std::make_shared<xRooNLLVar>(nullptr, nullptr);
834 auto dummyNll = std::make_shared<RooRealVar>(apath, "Dummy NLL", std::numeric_limits<double>::quiet_NaN());
835 nll->std::shared_ptr<RooAbsReal>::operator=(dummyNll);
836 dummyNll->setAttribute("readOnly");
837 // add pars as 'servers' on the dummy NLL
838 if (fPars) {
839 for (auto &&p : *fPars) {
840 dummyNll->addServer(
841 *p); // this is ok provided fPars (i.e. hypoSpace) stays alive as long as the hypoPoint ...
842 }
843 // flag poi
844 for (auto &p : allpois) {
845 fPars->find(p.c_str())->setAttribute("poi", true);
846 }
847 }
848 nll->reinitialize(); // triggers filling of par lists etc
849
850 for (auto &&[key, value] : cfits) {
851 if (value.find(allpois) != value.end()) {
852 // get the value of the poi in the key set
853 auto _coords = std::make_shared<RooArgSet>();
854 for (auto &k : key) {
855 auto v = _coords->addClone(RooRealVar(std::get<0>(k).c_str(), std::get<0>(k).c_str(), std::get<1>(k)));
856 v->setAttribute("poi", allpois.find(std::get<0>(k)) != allpois.end());
857 if (!std::get<2>(k).empty()) {
858 v->setStringAttribute("altVal", std::get<2>(k).c_str());
859 }
860 }
862 // hp.fPOIName = allpois.begin()->c_str();
863 // hp.fNullVal = _coords->getRealValue(hp.fPOIName.c_str());
864 hp.coords = _coords;
865 hp.nllVar = nll;
866
867 // auto altVal =
868 // hp.null_cfit()->constPars().find(hp.fPOIName.c_str())->getStringAttribute("altVal");
869 // if(altVal) hp.fAltVal = TString(altVal).Atof();
870 // else hp.fAltVal = std::numeric_limits<double>::quiet_NaN();
871
872 // decide based on values
873 if (std::isnan(hp.fAltVal())) {
875 } else if (hp.fNullVal() >= hp.fAltVal()) {
877 } else {
879 }
880
881 emplace_back(hp);
882 }
883 }
884 } else if (nFits > 0) {
885 std::cout << "possible POI: ";
886 for (auto p : allpois)
887 std::cout << p << ",";
888 std::cout << std::endl;
889 }
890}
891
893{
894
895 auto _axes = axes();
896
897 size_t badFits = 0;
898
899 for (size_t i = 0; i < size(); i++) {
900 std::cout << i << ") ";
901 for (auto a : _axes) {
902 if (a != _axes.first())
903 std::cout << ",";
904 std::cout << a->GetName() << "="
905 << at(i).coords->getRealValue(a->GetName(), std::numeric_limits<double>::quiet_NaN());
906 }
907 std::cout << " status=[ufit:";
908 auto ufit = const_cast<xRooHypoPoint &>(at(i)).ufit(true);
909 if (!ufit) {
910 std::cout << "-";
911 } else {
912 std::cout << ufit->status();
913 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(ufit->status()) == 0);
914 }
915 std::cout << ",cfit_null:";
916 auto cfit = const_cast<xRooHypoPoint &>(at(i)).cfit_null(true);
917 if (!cfit) {
918 std::cout << "-";
919 } else {
920 std::cout << cfit->status();
921 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(cfit->status()) == 0);
922 }
923 std::cout << ",cfit_alt:";
924 auto afit = const_cast<xRooHypoPoint &>(at(i)).cfit_alt(true);
925 if (!afit) {
926 std::cout << "-";
927 } else {
928 std::cout << afit->status();
929 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(afit->status()) == 0);
930 }
931 if (auto asiPoint = const_cast<xRooHypoPoint &>(at(i)).asimov(true)) {
932 std::cout << ",asimov.ufit:";
933 auto asi_ufit = asiPoint->ufit(true);
934 if (!asi_ufit) {
935 std::cout << "-";
936 } else {
937 std::cout << asi_ufit->status();
938 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(asi_ufit->status()) == 0);
939 }
940 std::cout << ",asimov.cfit_null:";
941 auto asi_cfit = asiPoint->cfit_null(true);
942 if (!asi_cfit) {
943 std::cout << "-";
944 } else {
945 std::cout << asi_cfit->status();
946 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(asi_cfit->status()) == 0);
947 }
948 }
949 std::cout << "]";
950 auto sigma_mu = const_cast<xRooHypoPoint &>(at(i)).sigma_mu(true);
951 if (!std::isnan(sigma_mu.first)) {
952 std::cout << " sigma_mu=" << sigma_mu.first;
953 if (sigma_mu.second)
954 std::cout << " +/- " << sigma_mu.second;
955 }
956 std::cout << std::endl;
957 }
958 std::cout << "--------------------------" << std::endl;
959 std::cout << "Number of bad fits: " << badFits << std::endl;
960}
961
962std::shared_ptr<TGraphErrors> xRooNLLVar::xRooHypoSpace::graph(
963 const char *opt /*, const std::function<void(xRooNLLVar::xRooHypoSpace*)>& progress*/) const
964{
965
966 TString sOpt(opt);
967 sOpt.ToLower();
968
969 bool doCLs = sOpt.Contains("cls");
970 bool readOnly = sOpt.Contains("readonly");
971 bool visualize = sOpt.Contains("visualize") && !readOnly;
972
973 double nSigma =
974 (sOpt.Contains("exp"))
975 ? (TString(sOpt(sOpt.Index("exp") + 3,
976 sOpt.Index(" ", sOpt.Index("exp")) == -1 ? sOpt.Length() : sOpt.Index(" ", sOpt.Index("exp"))))
977 .Atof())
978 : std::numeric_limits<double>::quiet_NaN();
979 bool expBand =
980 !std::isnan(nSigma) && nSigma && !(sOpt(sOpt.Index("exp") + 3) == '+' || sOpt(sOpt.Index("exp") + 3) == '-');
981
982 auto _axes = axes();
983 if (_axes.size() != 1)
984 return nullptr;
985
986 auto out = std::make_shared<TGraphErrors>();
987 out->SetName(GetName());
988 out->SetEditable(false);
989 const char *sCL = (doCLs) ? "CLs" : "null";
990
991 TString title =
992 TString::Format("%s;%s;p_{%s}",
993 (std::isnan(nSigma))
994 ? "Observed"
995 : (!nSigma ? "Expected"
996 : TString::Format("%s%d#sigma Expected",
997 expBand || !nSigma ? "" : ((nSigma < 0) ? "-" : "+"), int(nSigma))
998 .Data()),
999 _axes.at(0)->GetTitle(), sCL);
1000
1001 if (std::isnan(nSigma)) {
1002 out->SetNameTitle(TString::Format("obs_p%s", sCL), title);
1003 out->SetMarkerStyle(20);
1004 out->SetMarkerSize(0.5);
1005 if (sOpt.Contains("ts")) {
1006 out->SetNameTitle("obs_ts", TString::Format("Observed;%s;%s", _axes.at(0)->GetTitle(),
1007 (empty() ? "" : front().tsTitle(true).Data())));
1008 }
1009 } else {
1010 out->SetNameTitle(TString::Format("exp%d_p%s", int(nSigma), sCL), title);
1011 out->SetMarkerStyle(0);
1012 out->SetMarkerSize(0);
1013 out->SetLineStyle(2 + int(nSigma));
1014 if (expBand && nSigma) {
1015 out->SetFillColor((nSigma == 2) ? kYellow : kGreen);
1016 // out->SetFillStyle(3005);
1017 out->SetLineStyle(0);
1018 out->SetLineWidth(0);
1019 auto x = out->Clone("up");
1020 x->SetBit(kCanDelete);
1021 dynamic_cast<TAttFill *>(x)->SetFillStyle(nSigma == 2 ? 3005 : 3004);
1022 dynamic_cast<TAttFill *>(x)->SetFillColor(kBlack);
1023 out->GetListOfFunctions()->Add(x, "F");
1024 x = out->Clone("down");
1025 x->SetBit(kCanDelete);
1026 // dynamic_cast<TAttFill*>(x)->SetFillColor((nSigma==2) ? kYellow : kGreen);
1027 // dynamic_cast<TAttFill*>(x)->SetFillStyle(1001);
1028 out->GetListOfFunctions()->Add(x, "F");
1029 }
1030 if (sOpt.Contains("ts")) {
1031 out->SetNameTitle(TString::Format("exp_ts%d", int(nSigma)),
1032 TString::Format("Expected;%s;%s", _axes.at(0)->GetTitle(), front().tsTitle(true).Data()));
1033 }
1034 }
1035
1036 auto badPoints = [&]() {
1037 auto badPoints2 = dynamic_cast<TGraph *>(out->GetListOfFunctions()->FindObject("badPoints"));
1038 if (!badPoints2) {
1039 badPoints2 = new TGraph;
1040 badPoints2->SetBit(kCanDelete);
1041 badPoints2->SetName("badPoints");
1042 badPoints2->SetMarkerStyle(5);
1043 badPoints2->SetMarkerColor(std::isnan(nSigma) ? kRed : kBlue);
1044 badPoints2->SetMarkerSize(1);
1045 out->GetListOfFunctions()->Add(badPoints2, "P");
1046 }
1047 return badPoints2;
1048 };
1049 int nPointsDown = 0;
1050 bool above = true;
1051 TStopwatch s;
1052 s.Start();
1053 size_t nDone = 0;
1054 for (auto &p : *this) {
1055 if (s.RealTime() > 5) {
1056 if (visualize) {
1057 // draw readonly version of the graph
1058 auto gra = graph(sOpt + " readOnly");
1059 if (gra && gra->GetN()) {
1060 if (!gPad && gROOT->GetSelectedPad())
1061 gROOT->GetSelectedPad()->cd();
1062 if (gPad)
1063 gPad->Clear();
1064 gra->DrawClone(expBand ? "AF" : "ALP")->SetBit(kCanDelete);
1065#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1066 if (auto pad = gROOT->GetSelectedPad()) {
1067 pad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1068 }
1069#endif
1071 }
1072 } else {
1073 ::Info("xRooHypoSpace::graph", "Completed %d/%d points for %s", int(nDone), int(size()), sOpt.Data());
1074 }
1075 s.Start();
1076 } else {
1077 s.Continue();
1078 }
1079 double _x = p.coords->getRealValue(_axes.at(0)->GetName(), std::numeric_limits<double>::quiet_NaN());
1080 auto pval = const_cast<xRooHypoPoint &>(p).getVal(sOpt);
1081 auto idx = out->GetN() - nPointsDown;
1082
1083 if (std::isnan(pval.first)) {
1084 if (p.status() != 0) { // if status is 0 then bad pval is really just absence of fits, not bad fits
1085 badPoints()->SetPoint(badPoints()->GetN(), _x, 0);
1086 }
1087 } else {
1088 out->InsertPointBefore(idx, _x, pval.first);
1089 out->SetPointError(idx, 0, pval.second);
1090 }
1091
1092 if (expBand && nSigma) {
1093 TString sOpt2 = sOpt;
1094 sOpt2.ReplaceAll("exp", "exp-");
1095 pval = const_cast<xRooHypoPoint &>(p).getVal(sOpt2);
1096 if (std::isnan(pval.first)) {
1097 if (p.status() != 0) { // if status is 0 then bad pval is really just absence of fits, not bad fits
1098 badPoints()->SetPoint(badPoints()->GetN(), _x, 0);
1099 }
1100 } else {
1101 out->InsertPointBefore(idx + 1, _x, pval.first);
1102 out->SetPointError(idx + 1, 0, pval.second);
1103 nPointsDown++;
1104 if (out->GetPointY(idx) < pval.first)
1105 above = false; // the -sigma points are actually above +sigma
1106 }
1107 }
1108 nDone++;
1109 }
1110
1111 if (out->GetN() == 0)
1112 return out;
1113
1114 if (!expBand) {
1115 out->Sort();
1116 if (out->GetListOfFunctions()->FindObject("badPoints")) {
1117 // try to interpolate the points
1118 for (int i = 0; i < badPoints()->GetN(); i++) {
1119 badPoints()->SetPointY(i, out->Eval(badPoints()->GetPointX(i)));
1120 }
1121 }
1122 } else {
1123 out->Sort(&TGraph::CompareX, true, 0, out->GetN() - nPointsDown - 1); // sort first half
1124 out->Sort(&TGraph::CompareX, false, out->GetN() - nPointsDown, out->GetN() - 1); // reverse sort second half
1125
1126 // now populate the up and down values
1127 auto up = dynamic_cast<TGraph *>(out->GetListOfFunctions()->FindObject("up"));
1128 auto down = dynamic_cast<TGraph *>(out->GetListOfFunctions()->FindObject("down"));
1129
1130 for (int i = 0; i < out->GetN(); i++) {
1131 if (i < out->GetN() - nPointsDown) {
1132 up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) + out->GetErrorY(i) * (above ? 1. : -1.));
1133 down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.));
1134 } else {
1135 up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.));
1136 down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) + out->GetErrorY(i) * (above ? 1. : -1.));
1137 }
1138 }
1139 }
1140
1141 if (visualize) {
1142 // draw result
1143 if (!gPad && gROOT->GetSelectedPad())
1144 gROOT->GetSelectedPad()->cd();
1145 if (gPad)
1146 gPad->Clear();
1147 out->DrawClone(expBand ? "AF" : "ALP")->SetBit(kCanDelete);
1148#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1149 if (auto pad = gROOT->GetSelectedPad()) {
1150 pad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1151 }
1152#endif
1154 }
1155
1156 return out;
1157}
1158
1159std::shared_ptr<TMultiGraph> xRooNLLVar::xRooHypoSpace::graphs(const char *opt)
1160{
1161 TString sOpt(opt);
1162 sOpt.ToLower();
1163 std::shared_ptr<TMultiGraph> out;
1164 if (sOpt.Contains("pcls") || sOpt.Contains("pnull") || sOpt.Contains("ts")) {
1165
1166 bool visualize = sOpt.Contains("visualize");
1167 sOpt.ReplaceAll("visualize", "");
1168
1169 auto exp2 = graph(sOpt + " exp2");
1170 auto exp1 = graph(sOpt + " exp1");
1171 auto exp = graph(sOpt + " exp");
1172 bool doObs = true;
1173 // for(auto& hp : *this) { if(hp.isExpected) {doObs=false; break;} }
1174 auto obs = (doObs) ? graph(sOpt) : nullptr;
1175
1176 out = std::make_shared<TMultiGraph>(GetName(), GetTitle());
1177 if (exp2 && exp2->GetN() > 1)
1178 out->Add(static_cast<TGraph *>(exp2->Clone()), "FP");
1179 if (exp1 && exp1->GetN() > 1)
1180 out->Add(static_cast<TGraph *>(exp1->Clone()), "FP");
1181 if (exp && exp->GetN() > 1)
1182 out->Add(static_cast<TGraph *>(exp->Clone()), "LP");
1183 if (obs && obs->GetN() > 1)
1184 out->Add(static_cast<TGraph *>(obs->Clone()), "LP");
1185
1186 if (!out->GetListOfGraphs()) {
1187 return nullptr;
1188 }
1189
1190 TGraph *testedPoints = nullptr;
1191 if (sOpt.Contains("pcls")) {
1192 TGraph *line = new TGraph;
1193 line->SetName("alpha");
1194 line->SetLineStyle(2);
1195 line->SetEditable(false);
1196 line->SetPoint(line->GetN(), out->GetHistogram()->GetXaxis()->GetXmin() - 10, 0.05);
1197 testedPoints = new TGraph;
1198 testedPoints->SetName("hypoPoints");
1199 testedPoints->SetEditable(false);
1200 testedPoints->SetMarkerStyle(24);
1201 testedPoints->SetMarkerSize(0.4); // use line to indicate tested points
1202 if (exp) {
1203 for (int i = 0; i < exp->GetN(); i++) {
1204 testedPoints->SetPoint(testedPoints->GetN(), exp->GetPointX(i), 0.05);
1205 }
1206 }
1207 line->SetPoint(line->GetN(), out->GetHistogram()->GetXaxis()->GetXmax() + 10, 0.05);
1209 out->GetListOfFunctions()->Add(line, "L");
1210 }
1211 if (exp) {
1212 out->GetHistogram()->GetXaxis()->SetTitle(exp->GetHistogram()->GetXaxis()->GetTitle());
1213 out->GetHistogram()->GetYaxis()->SetTitle(exp->GetHistogram()->GetYaxis()->GetTitle());
1214 }
1215 auto leg = new TLegend(1. - gStyle->GetPadRightMargin() - 0.3, 1. - gStyle->GetPadTopMargin() - 0.35,
1216 1. - gStyle->GetPadRightMargin() - 0.05, 1. - gStyle->GetPadTopMargin() - 0.05);
1217 leg->SetName("legend");
1218 leg->SetBit(kCanDelete);
1219
1220 out->GetListOfFunctions()->Add(leg);
1221 // out->GetListOfFunctions()->Add(out->GetHistogram()->Clone(".axis"),"sameaxis"); // redraw axis
1222
1223 for (auto g : *out->GetListOfGraphs()) {
1224 if (auto o = dynamic_cast<TGraph *>(g)->GetListOfFunctions()->FindObject("down")) {
1225 leg->AddEntry(o, "", "F");
1226 } else {
1227 leg->AddEntry(g, "", "LPE");
1228 }
1229 }
1230
1231 if (sOpt.Contains("pcls")) {
1232 // add current limit estimates to legend
1233 if (exp2 && exp2->GetN() > 1) {
1234 auto l = xRooFit::matchPrecision(GetLimit(*graph(sOpt + "exp-2")));
1235 leg->AddEntry((TObject *)nullptr, TString::Format("-2#sigma: %g +/- %g", l.first, l.second), "");
1236 }
1237 if (exp1 && exp1->GetN() > 1) {
1238 auto l = xRooFit::matchPrecision(GetLimit(*graph(sOpt + "exp-1")));
1239 leg->AddEntry((TObject *)nullptr, TString::Format("-1#sigma: %g +/- %g", l.first, l.second), "");
1240 }
1241 if (exp && exp->GetN() > 1) {
1242 auto l = xRooFit::matchPrecision(GetLimit(*exp));
1243 leg->AddEntry((TObject *)nullptr, TString::Format("0#sigma: %g +/- %g", l.first, l.second), "");
1244 }
1245 if (exp1 && exp1->GetN() > 1) {
1246 auto l = xRooFit::matchPrecision(GetLimit(*graph(sOpt + "exp+1")));
1247 leg->AddEntry((TObject *)nullptr, TString::Format("+1#sigma: %g +/- %g", l.first, l.second), "");
1248 }
1249 if (exp2 && exp2->GetN() > 1) {
1250 auto l = xRooFit::matchPrecision(GetLimit(*graph(sOpt + "exp+2")));
1251 leg->AddEntry((TObject *)nullptr, TString::Format("+2#sigma: %g +/- %g", l.first, l.second), "");
1252 }
1253 if (obs && obs->GetN() > 1) {
1254 auto l = xRooFit::matchPrecision(GetLimit(*obs));
1255 leg->AddEntry((TObject *)nullptr, TString::Format("Observed: %g +/- %g", l.first, l.second), "");
1256 }
1257 }
1258 if (testedPoints)
1259 out->Add(testedPoints, "P");
1260
1261 if (visualize) {
1262 if (!gPad && gROOT->GetSelectedPad())
1263 gROOT->GetSelectedPad()->cd();
1264 if (gPad)
1265 gPad->Clear();
1266 auto gra2 = static_cast<TMultiGraph *>(out->DrawClone("A"));
1267 gra2->SetBit(kCanDelete);
1268 if (sOpt.Contains("pcls") || sOpt.Contains("pnull")) {
1269 gra2->GetHistogram()->SetMinimum(1e-6);
1270 }
1271 if (gPad) {
1272 gPad->RedrawAxis();
1273 gPad->GetCanvas()->Paint();
1274 gPad->GetCanvas()->Update();
1275#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1276 gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1277#endif
1278 }
1280 }
1281 }
1282
1283 return out;
1284}
1285
1286std::pair<double, double> xRooNLLVar::xRooHypoSpace::GetLimit(const TGraph &pValues, double target)
1287{
1288
1289 if (std::isnan(target)) {
1290 target = 1. - gEnv->GetValue("xRooHypoSpace.CL", 95.) / 100.;
1291 }
1292
1293 auto gr = std::make_shared<TGraph>(pValues);
1294 // remove any nan points and duplicates
1295 int i = 0;
1296 std::set<double> existingX;
1297 while (i < gr->GetN()) {
1298 if (std::isnan(gr->GetPointY(i))) {
1299 gr->RemovePoint(i);
1300 } else if (existingX.find(gr->GetPointX(i)) != existingX.end()) {
1301 gr->RemovePoint(i);
1302 } else {
1303 existingX.insert(gr->GetPointX(i));
1304 // convert to log ....
1305 gr->SetPointY(i, log(std::max(gr->GetPointY(i), 1e-10)));
1306 i++;
1307 }
1308 }
1309
1310 gr->Sort();
1311
1312 // simple linear extrapolation to critical value ... return nan if problem
1313 if (gr->GetN() < 2) {
1314 return std::pair(std::numeric_limits<double>::quiet_NaN(), 0);
1315 }
1316
1317 double alpha = log(target);
1318
1319 bool above = gr->GetPointY(0) > alpha;
1320 for (int ii = 1; ii < gr->GetN(); ii++) {
1321 if ((above && (gr->GetPointY(ii) <= alpha)) || (!above && (gr->GetPointY(ii) >= alpha))) {
1322 // found the limit ... return linearly extrapolated point
1323 double lim = gr->GetPointX(ii - 1) + (gr->GetPointX(ii) - gr->GetPointX(ii - 1)) *
1324 (alpha - gr->GetPointY(ii - 1)) /
1325 (gr->GetPointY(ii) - gr->GetPointY(ii - 1));
1326 // use points either side as error
1327 double err = std::max(lim - gr->GetPointX(ii - 1), gr->GetPointX(ii) - lim);
1328 // give err negative sign to indicate if error due to negative side
1329 if ((lim - gr->GetPointX(ii - 1)) > (gr->GetPointX(ii) - lim))
1330 err *= -1;
1331 return std::pair(lim, err);
1332 }
1333 }
1334 // if reach here need to extrapolate ...
1335 if ((above && gr->GetPointY(gr->GetN() - 1) <= gr->GetPointY(0)) ||
1336 (!above && gr->GetPointY(gr->GetN() - 1) >= gr->GetPointY(0))) {
1337 // extrapolating above based on last two points
1338 // in fact, if 2nd last point is a p=1 (log(p)=0) then go back
1339 int offset = 2;
1340 while (offset < gr->GetN() && gr->GetPointY(gr->GetN() - offset) == 0)
1341 offset++;
1342 double x1 = gr->GetPointX(gr->GetN() - offset);
1343 double y1 = gr->GetPointY(gr->GetN() - offset);
1344 double m = (gr->GetPointY(gr->GetN() - 1) - y1) / (gr->GetPointX(gr->GetN() - 1) - x1);
1345 if (m == 0.)
1346 return std::pair(std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity());
1347 return std::pair((alpha - y1) / m + x1, std::numeric_limits<double>::infinity());
1348 } else {
1349 // extrapolating below based on first two points
1350 double x1 = gr->GetPointX(0);
1351 double y1 = gr->GetPointY(0);
1352 double m = (gr->GetPointY(1) - y1) / (gr->GetPointX(1) - x1);
1353 if (m == 0.)
1354 return std::pair(-std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity());
1355 return std::pair((alpha - y1) / m + x1, -std::numeric_limits<double>::infinity());
1356 }
1357}
1358
1360{
1361 TString sOpt = TString::Format("p%s", type);
1362 if (std::isnan(nSigma)) {
1363 sOpt += "obs";
1364 } else {
1365 sOpt += TString::Format("exp%s%d", nSigma > 0 ? "+" : "", int(nSigma));
1366 }
1367 return GetLimit(*graph(sOpt + " readonly"));
1368}
1369
1371xRooNLLVar::xRooHypoSpace::findlimit(const char *opt, double relUncert, unsigned int maxTries)
1372{
1373 TString sOpt(opt);
1374 bool visualize = sOpt.Contains("visualize");
1375 sOpt.ReplaceAll("visualize", "");
1376 std::shared_ptr<TGraphErrors> gr = graph(sOpt + " readonly");
1377 if (visualize) {
1378 auto gra = graphs(sOpt.Contains("toys") ? "pcls readonly toys" : "pcls readonly");
1379 if (gra) {
1380 if (!gPad)
1381 gra->Draw(); // in 6.28 DrawClone wont make the gPad defined :( ... so Draw then clear and Draw Clone
1382 gPad->Clear();
1383 gra->DrawClone("A")->SetBit(kCanDelete);
1384 gPad->RedrawAxis();
1385 gra->GetHistogram()->SetMinimum(1e-9);
1386 gra->GetHistogram()->GetYaxis()->SetRangeUser(1e-9, 1);
1387 gPad->Modified();
1388#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1389 gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1390#endif
1392 }
1393 }
1394
1395 // resync parameter boundaries from nlls (may have been modified by fits)
1396 for (auto p : axes()) {
1397 for (auto &[pdf, nll] : fNlls) {
1398 if (auto _v = dynamic_cast<RooRealVar *>(nll->pars()->find(*p))) {
1399 dynamic_cast<RooRealVar *>(p)->setRange(_v->getMin(), _v->getMax());
1400 }
1401 }
1402 }
1403
1404 if (!gr || gr->GetN() < 2) {
1405 auto v = (axes().empty()) ? nullptr : dynamic_cast<RooRealVar *>(*axes().rbegin());
1406 if (!v)
1407 return std::pair(std::numeric_limits<double>::quiet_NaN(), 0.);
1408 double muMax = std::min(std::min(v->getMax("physical"), v->getMax()), v->getMax("scan"));
1409 double muMin = std::max(std::max(v->getMin("physical"), v->getMin()), v->getMin("scan"));
1410 if (!gr || gr->GetN() < 1) {
1411 if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), muMin)).getVal(sOpt).first)) {
1412 // first point failed ... give up
1413 Error("findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), muMin);
1414 return std::pair(std::numeric_limits<double>::quiet_NaN(), 0.);
1415 }
1416 gr.reset();
1417 return findlimit(opt, relUncert, maxTries - 1); // do this to resync parameter limits
1418 }
1419
1420 // can approximate expected limit using
1421 // mu_hat + sigma_mu*ROOT::Math::gaussian_quantile(1.-alpha/2.,1) for cls
1422 // or mu_hat + sigma_mu*ROOT::Math::gaussian_quantile((1.-alpha),1) for cls+b
1423 // get a very first estimate of sigma_mu from ufit to expected data, take error on mu as sigma_mu
1424 double nextPoint = muMin + (muMax - muMin) / 50;
1425
1426 // if done an expected limit, assume data is like expected and choose expected limit point as first test point
1427 if (sOpt.Contains("obs")) {
1428 TString sOpt2 = sOpt;
1429 sOpt2.ReplaceAll("obs", "exp");
1430 auto expLim = findlimit(sOpt2, std::numeric_limits<double>::infinity(), 0);
1431 if (!std::isnan(expLim.first) && expLim.first < nextPoint)
1432 nextPoint = expLim.first;
1433 }
1434
1435 auto point =
1436 (sOpt.Contains("exp")) ? back().asimov() : std::shared_ptr<xRooHypoPoint>(&back(), [](xRooHypoPoint *) {});
1437 point = nullptr;
1438 if (point && point->ufit()) {
1439 double rough_sigma_mu = point->mu_hat().getError();
1440 double another_estimate = point->mu_hat().getVal() + rough_sigma_mu * ROOT::Math::gaussian_quantile(0.95, 1);
1441 // if (another_estimate < nextPoint) {
1442 nextPoint = another_estimate;
1443 ::Info("xRooHypoSpace::findlimit", "Guessing %g based on rough sigma_mu = %g", nextPoint, rough_sigma_mu);
1444 //}
1445 }
1446
1447 if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), nextPoint)).getVal(sOpt).first)) {
1448 // second point failed ... give up
1449 Error("findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), nextPoint);
1450 return std::pair(std::numeric_limits<double>::quiet_NaN(), 0.);
1451 }
1452 gr.reset();
1453 return findlimit(opt, relUncert, maxTries - 1);
1454 }
1455
1456 auto lim = GetLimit(*gr);
1457
1458 if (std::isnan(lim.first)) {
1459 return lim;
1460 }
1461
1462 auto v = dynamic_cast<RooRealVar *>(*axes().rbegin());
1463 double maxMu = std::min(std::min(v->getMax("physical"), v->getMax()), v->getMax("scan"));
1464 double minMu = std::max(std::max(v->getMin("physical"), v->getMin()), v->getMin("scan"));
1465
1466 // static double MIN_LIMIT_UNCERT = 1e-4; // stop iterating once uncert gets this small
1467 if (lim.first > -std::numeric_limits<double>::infinity() && lim.first < std::numeric_limits<double>::infinity() &&
1468 (std::abs(lim.second) <= relUncert * std::abs(lim.first) /* || std::abs(lim.second)<MIN_LIMIT_UNCERT*/))
1469 return lim;
1470
1471 double nextPoint;
1472
1473 if (lim.second == std::numeric_limits<double>::infinity()) {
1474 // limit was found by extrapolating to right
1475 nextPoint = lim.first;
1476 if (nextPoint == std::numeric_limits<double>::infinity() || nextPoint > maxMu) {
1477 nextPoint = gr->GetPointX(gr->GetN() - 1) + (maxMu - minMu) / 50;
1478 }
1479
1480 // prefer extrapolation with sigma_mu, if available, if it takes us further
1481 // as shape of p-value curve is usually
1482 auto point =
1483 (sOpt.Contains("exp")) ? back().asimov() : std::shared_ptr<xRooHypoPoint>(&back(), [](xRooHypoPoint *) {});
1484 point = nullptr;
1485 if (point && point->ufit()) {
1486 double rough_sigma_mu = point->mu_hat().getError();
1487 double another_estimate = point->mu_hat().getVal() + rough_sigma_mu * ROOT::Math::gaussian_quantile(0.95, 1);
1488 // if (another_estimate < nextPoint) {
1489 nextPoint = std::max(nextPoint, another_estimate);
1490 ::Info("xRooHypoSpace::findlimit", "Guessing %g based on rough sigma_mu = %g", nextPoint, rough_sigma_mu);
1491 //}
1492 }
1493 nextPoint = std::min(nextPoint + nextPoint * relUncert * 0.99, maxMu); // ensure we step over location if possible
1494
1495 if (nextPoint > maxMu)
1496 return lim;
1497 } else if (lim.second == -std::numeric_limits<double>::infinity()) {
1498 // limit from extrapolating to left
1499 nextPoint = lim.first;
1500 if (nextPoint < minMu)
1501 nextPoint = gr->GetPointX(0) - (maxMu - minMu) / 50;
1502 if (nextPoint < minMu)
1503 return lim;
1504 } else {
1505 nextPoint = lim.first + lim.second * relUncert * 0.99;
1506 }
1507
1508 // got here need a new point .... evaluate the estimated lim location +/- the relUncert (signed error takes care of
1509 // direction)
1510
1511 ::Info("xRooHypoSpace::findlimit", "%s -- Testing new point @ %s=%g (delta=%g)", sOpt.Data(), v->GetName(),
1512 nextPoint, lim.second);
1513 if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), nextPoint)).getVal(sOpt).first)) {
1514 if (maxTries == 0) {
1515 Warning("findlimit", "Reached max number of point evaluations");
1516 } else {
1517 Error("findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), nextPoint);
1518 }
1519 return lim;
1520 }
1521 gr.reset();
1522 return findlimit(opt, relUncert, maxTries - 1);
1523}
1524
1526{
1527
1528 TString sOpt(opt);
1529 sOpt.ToLower();
1530
1531 if ((sOpt == "" || sOpt == "same") && !empty()) {
1532 if (front().fPllType == xRooFit::Asymptotics::OneSidedPositive) {
1533 sOpt += "pcls"; // default to showing cls p-value scan if drawing a limit
1534 for (auto &hp : *this) {
1535 if (!hp.nullToys.empty() || !hp.altToys.empty()) {
1536 sOpt += " toys";
1537 break; // default to toys if done toys
1538 }
1539 }
1540 } else if (front().fPllType == xRooFit::Asymptotics::TwoSided) {
1541 sOpt += "ts";
1542 }
1543 }
1544
1545 // split up by ; and call Draw for each (with 'same' appended)
1546 auto _axes = axes();
1547 if (_axes.empty())
1548 return;
1549
1550 if (sOpt == "status") {
1551 // draw the points in the space
1552 if (_axes.size() <= 2) {
1553 TGraphErrors *out = new TGraphErrors;
1554 out->SetBit(kCanDelete);
1555 out->SetName("points");
1556 out->SetMarkerSize(0.5);
1557 TGraph *tsAvail = new TGraph;
1558 tsAvail->SetName("ts");
1559 tsAvail->SetBit(kCanDelete);
1560 tsAvail->SetMarkerStyle(20);
1561 TGraph *expAvail = new TGraph;
1562 expAvail->SetName("exp");
1563 expAvail->SetBit(kCanDelete);
1564 expAvail->SetMarkerStyle(25);
1565 expAvail->SetMarkerSize(out->GetMarkerSize() * 1.5);
1566 TGraph *badPoints = new TGraph;
1567 badPoints->SetName("bad_ufit");
1568 badPoints->SetBit(kCanDelete);
1569 badPoints->SetMarkerStyle(5);
1570 badPoints->SetMarkerColor(kRed);
1571 badPoints->SetMarkerSize(out->GetMarkerSize());
1572 TGraph *badPoints2 = new TGraph;
1573 badPoints2->SetName("bad_cfit_null");
1574 badPoints2->SetBit(kCanDelete);
1575 badPoints2->SetMarkerStyle(2);
1576 badPoints2->SetMarkerColor(kRed);
1577 badPoints2->SetMarkerSize(out->GetMarkerSize());
1578
1579 out->SetTitle(TString::Format("%s;%s;%s", GetTitle(), _axes.at(0)->GetTitle(),
1580 (_axes.size() == 1) ? "" : _axes.at(1)->GetTitle()));
1581 for (auto &p : *this) {
1582 bool _readOnly = p.nllVar ? p.nllVar->get()->getAttribute("readOnly") : false;
1583 if (p.nllVar)
1584 p.nllVar->get()->setAttribute("readOnly", true);
1585 double x = p.coords->getRealValue(_axes.at(0)->GetName());
1586 double y = _axes.size() == 1 ? p.ts_asymp().first : p.coords->getRealValue(_axes.at(1)->GetName());
1587 out->SetPoint(out->GetN(), x, y);
1588 if (!std::isnan(p.ts_asymp().first)) {
1589 if (_axes.size() == 1)
1590 out->SetPointError(out->GetN() - 1, 0, p.ts_asymp().second);
1591 tsAvail->SetPoint(tsAvail->GetN(), x, y);
1592 } else if (p.fUfit && (std::isnan(p.fUfit->minNll()) ||
1593 xRooNLLVar::xRooHypoPoint::allowedStatusCodes.find(p.fUfit->status()) ==
1595 badPoints->SetPoint(badPoints->GetN(), x, y);
1596 } else if (p.fNull_cfit && (std::isnan(p.fNull_cfit->minNll()) ||
1597 xRooNLLVar::xRooHypoPoint::allowedStatusCodes.find(p.fNull_cfit->status()) ==
1599 badPoints2->SetPoint(badPoints2->GetN(), x, y);
1600 }
1601 if (!std::isnan(p.ts_asymp(0).first)) {
1602 expAvail->SetPoint(expAvail->GetN(), x, y);
1603 } else if (p.asimov() && p.asimov()->fUfit &&
1604 (std::isnan(p.asimov()->fUfit->minNll()) ||
1605 xRooNLLVar::xRooHypoPoint::allowedStatusCodes.find(p.asimov()->fUfit->status()) ==
1607
1608 } else if (p.asimov() && p.asimov()->fNull_cfit &&
1609 (std::isnan(p.asimov()->fNull_cfit->minNll()) ||
1610 xRooNLLVar::xRooHypoPoint::allowedStatusCodes.find(p.asimov()->fNull_cfit->status()) ==
1612 }
1613 if (p.nllVar)
1614 p.nllVar->get()->setAttribute("readOnly", _readOnly);
1615 }
1616
1617 if (_axes.size() == 1) {
1618 TGraph tmp;
1619 for (int i = 0; i < out->GetN(); i++) {
1620 if (!std::isnan(out->GetPointY(i)))
1621 tmp.SetPoint(tmp.GetN(), out->GetPointX(i), out->GetPointY(i));
1622 }
1623 auto fixPoints = [&](TGraph *g) {
1624 for (int i = 0; i < g->GetN(); i++) {
1625 if (std::isnan(g->GetPointY(i)))
1626 g->SetPointY(i, std::isnan(tmp.Eval(g->GetPointX(i))) ? 0. : tmp.Eval(g->GetPointX(i)));
1627 }
1628 };
1629 fixPoints(out);
1630 fixPoints(tsAvail);
1631 fixPoints(expAvail);
1632 fixPoints(badPoints);
1633 fixPoints(badPoints2);
1634 }
1635
1636 out->SetMarkerStyle(4);
1637 out->Draw("AP");
1638 auto leg = new TLegend(1. - gPad->GetRightMargin() - 0.3, 1. - gPad->GetTopMargin() - 0.35,
1639 1. - gPad->GetRightMargin() - 0.05, 1. - gPad->GetTopMargin() - 0.05);
1640 leg->SetName("legend");
1641 leg->AddEntry(out, "Uncomputed", "P");
1642
1643 if (tsAvail->GetN()) {
1644 out->GetListOfFunctions()->Add(tsAvail, "P");
1645 leg->AddEntry(tsAvail, "Computed", "P");
1646 } else {
1647 delete tsAvail;
1648 }
1649 if (expAvail->GetN()) {
1650 out->GetListOfFunctions()->Add(expAvail, "P");
1651 leg->AddEntry(expAvail, "Expected computed", "P");
1652 } else {
1653 delete expAvail;
1654 }
1655 if (badPoints->GetN()) {
1656 out->GetListOfFunctions()->Add(badPoints, "P");
1657 leg->AddEntry(badPoints, "Bad ufit", "P");
1658 } else {
1659 delete badPoints;
1660 }
1661 if (badPoints2->GetN()) {
1662 out->GetListOfFunctions()->Add(badPoints2, "P");
1663 leg->AddEntry(badPoints2, "Bad null cfit", "P");
1664 } else {
1665 delete badPoints2;
1666 }
1668 leg->Draw();
1669 // if(_axes.size()==1) out->GetHistogram()->GetYaxis()->SetRangeUser(0,1);
1670 gPad->SetGrid(true, _axes.size() > 1);
1671 if (_axes.size() == 1)
1672 gPad->SetLogy(false);
1673 }
1674
1676
1677 return;
1678 }
1679
1680 if (sOpt.Contains("pcls") || sOpt.Contains("pnull") || sOpt.Contains("ts")) {
1681 auto gra = graphs(sOpt + " readonly");
1682 if (!gPad && gROOT->GetSelectedPad())
1683 gROOT->GetSelectedPad()->cd();
1684 if (!sOpt.Contains("same") && gPad) {
1685 gPad->Clear();
1686 }
1687 if (gra) {
1688 auto gra2 = static_cast<TMultiGraph *>(gra->DrawClone(sOpt.Contains("same") ? "" : "A"));
1689 gra2->SetBit(kCanDelete);
1690 if (sOpt.Contains("pcls") || sOpt.Contains("pnull")) {
1691 gra2->GetHistogram()->SetMinimum(1e-6);
1692 }
1693 if (gPad) {
1694 gPad->RedrawAxis();
1695 gPad->GetCanvas()->Paint();
1696 gPad->GetCanvas()->Update();
1697#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1698 gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1699#endif
1701 }
1702 }
1703 if (!sOpt.Contains("same") && gPad) {
1704 // auto mg = static_cast<TMultiGraph*>(gPad->GetPrimitive(gra->GetName()));
1705 // mg->GetHistogram()->SetMinimum(1e-9);
1706 // mg->GetHistogram()->GetYaxis()->SetRangeUser(1e-9,1);
1707 // gPad->SetGrid(0, 0);
1708 // gPad->SetLogy(1);
1709 }
1710
1712
1713 return;
1714 }
1715
1716 // graphs("ts visualize");return;
1717
1718 TGraphErrors *out = new TGraphErrors;
1719 out->SetName(GetName());
1720
1721 TString title = (!axes().empty()) ? TString::Format(";%s", axes().first()->GetTitle()) : "";
1722
1723 auto pllType = xRooFit::Asymptotics::TwoSided;
1724 if (!empty() && axes().size() == 1) {
1725 for (auto &p : *this) {
1726 if (p.fPllType != xRooFit::Asymptotics::TwoSided) {
1727 pllType = p.fPllType;
1728 }
1729 }
1730 title += ";";
1731 title += front().tsTitle(true);
1732 }
1733
1734 out->SetTitle(title);
1735 *dynamic_cast<TAttFill *>(out) = *this;
1736 *dynamic_cast<TAttLine *>(out) = *this;
1737 *dynamic_cast<TAttMarker *>(out) = *this;
1738 out->SetBit(kCanDelete);
1739
1740 if (!gPad)
1742 TVirtualPad *basePad = gPad;
1743 if (!sOpt.Contains("same"))
1744 basePad->Clear();
1745
1746 bool doFits = false;
1747 if (sOpt.Contains("fits")) {
1748 doFits = true;
1749 sOpt.ReplaceAll("fits", "");
1750 }
1751
1752 auto mainPad = gPad;
1753
1754 out->SetEditable(false);
1755
1756 if (doFits) {
1757 gPad->Divide(1, 2);
1758 gPad->cd(1);
1759 gPad->SetBottomMargin(gPad->GetBottomMargin() * 2.); // increase margin to be same as before
1760 gPad->SetGrid();
1761 out->Draw(sOpt);
1762 basePad->cd(2);
1763 mainPad = basePad->GetPad(1);
1764 } else {
1765 gPad->SetGrid();
1766 out->Draw("ALP");
1767 }
1768
1769 std::pair<double, double> minMax(std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity());
1770 for (auto &p : *this) {
1771 if (p.fPllType != pllType)
1772 continue; // must all have same pll type
1773 auto val = p.pll(true).first;
1774 if (std::isnan(val))
1775 continue;
1776 minMax.first = std::min(minMax.first, val);
1777 minMax.second = std::max(minMax.second, val);
1778 }
1779 if (minMax.first < std::numeric_limits<double>::infinity())
1780 out->GetHistogram()->SetMinimum(minMax.first);
1781 if (minMax.second > -std::numeric_limits<double>::infinity())
1782 out->GetHistogram()->SetMaximum(minMax.second);
1783
1784 TGraph *badPoints = nullptr;
1785
1786 TStopwatch s;
1787 s.Start();
1788 std::shared_ptr<const RooFitResult> ufr;
1789 for (auto &p : *this) {
1790 if (p.fPllType != pllType)
1791 continue; // must all have same pll type
1792 auto val = p.pll().first;
1793 if (!ufr)
1794 ufr = p.ufit();
1795 if (out->GetN() == 0 && ufr && ufr->status() == 0) {
1796 out->SetPoint(out->GetN(),
1797 ufr->floatParsFinal().getRealValue(axes().first()->GetName(),
1798 ufr->constPars().getRealValue(axes().first()->GetName())),
1799 0.);
1800 out->SetPointError(out->GetN() - 1, 0, ufr->edm());
1801 }
1802 if (auto fr = p.fNull_cfit;
1803 fr && doFits) { // access member to avoid unnecessarily creating fit result if wasnt needed
1804 // create a new subpad and draw fitResult on it
1805 auto _pad = gPad;
1806 auto pad =
1807 new TPad(fr->GetName(), TString::Format("%s = %g", poi().first()->GetTitle(), p.fNullVal()), 0, 0, 1., 1);
1808 pad->SetNumber(out->GetN() + 1); // can't use "0" for a subpad
1809 pad->cd();
1810 xRooNode(fr).Draw("goff");
1811 _pad->cd();
1812 //_pad->GetListOfPrimitives()->AddFirst(pad);
1813 pad->AppendPad();
1814 }
1815 if (std::isnan(val) && p.status() != 0) {
1816 if (!badPoints) {
1817 badPoints = new TGraph;
1818 badPoints->SetBit(kCanDelete);
1819 badPoints->SetName("badPoints");
1820 badPoints->SetMarkerStyle(5);
1821 badPoints->SetMarkerColor(kRed);
1822 badPoints->SetMarkerSize(1);
1823 out->GetListOfFunctions()->Add(badPoints, "P");
1824 }
1825 badPoints->SetPoint(badPoints->GetN(), p.fNullVal(), out->Eval(p.fNullVal()));
1826 mainPad->Modified();
1827 } else if (!std::isnan(val)) {
1828 out->SetPoint(out->GetN(), p.coords->getRealValue(axes().first()->GetName()), p.pll().first);
1829 out->SetPointError(out->GetN() - 1, 0, p.pll().second);
1830 out->Sort();
1831
1832 // reposition bad points
1833 if (badPoints) {
1834 for (int i = 0; i < badPoints->GetN(); i++)
1835 badPoints->SetPointY(i, out->Eval(badPoints->GetPointX(i)));
1836 }
1837
1838 mainPad->Modified();
1839 }
1840 if (s.RealTime() > 3) { // stops the clock
1841 basePad->GetCanvas()->Paint();
1842 basePad->GetCanvas()->Update();
1844 s.Reset();
1845 s.Start();
1846 }
1847 s.Continue();
1848 }
1849 basePad->GetCanvas()->Paint();
1850 basePad->GetCanvas()->Update();
1851#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1852 basePad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1853#endif
1855
1856 // finish by overlaying ufit
1857 if (ufr && doFits) {
1858 auto _pad = gPad;
1859 auto pad = new TPad(ufr->GetName(), "unconditional fit", 0, 0, 1., 1.);
1860 pad->SetNumber(-1);
1861 pad->cd();
1862 xRooNode(ufr).Draw("goff");
1863 _pad->cd();
1864 pad->AppendPad();
1865
1866 // draw one more pad to represent the selected, and draw the ufit pad onto that pad
1867 pad = new TPad("selected", "selected", 0, 0, 1, 1);
1868 pad->Draw();
1869 pad->cd();
1870 basePad->GetPad(2)->GetPad(-1)->AppendPad();
1871 pad->Modified();
1872 pad->Update();
1874
1875 basePad->cd();
1876 }
1877
1878 if (doFits) {
1879 if (!xRooNode::gIntObj) {
1881 }
1882 gPad->GetCanvas()->Connect("Highlighted(TVirtualPad*,TObject*,Int_t,Int_t)", "xRooNode::InteractiveObject",
1883 xRooNode::gIntObj, "Interactive_PLLPlot(TVirtualPad*,TObject*,Int_t,Int_t)");
1884 }
1885
1886 return;
1887}
1888
1890{
1891
1892 RooStats::HypoTestInverterResult *out = nullptr;
1893
1894 auto _axes = axes();
1895 if (_axes.empty())
1896 return out;
1897
1898 out = new RooStats::HypoTestInverterResult(GetName(), *dynamic_cast<RooRealVar *>(_axes.at(0)), 0.95);
1899 out->SetTitle(GetTitle());
1900
1901 for (auto &p : *this) {
1902 double _x = p.coords->getRealValue(_axes.at(0)->GetName(), std::numeric_limits<double>::quiet_NaN());
1903 out->Add(_x, p.result());
1904 }
1905
1906 return out;
1907}
1908
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
const char Option_t
Definition RtypesCore.h:66
@ kRed
Definition Rtypes.h:66
@ kBlack
Definition Rtypes.h:65
@ kGreen
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
@ kYellow
Definition Rtypes.h:66
#define gDirectory
Definition TDirectory.h:384
R__EXTERN TEnv * gEnv
Definition TEnv.h:170
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 SetFillStyle
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 Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h offset
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 Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t target
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
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t SetFillColor
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 Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
Option_t Option_t TPoint TPoint const char y1
char name[80]
Definition TGX11.cxx:110
@ kCanDelete
Definition TObject.h:369
#define gROOT
Definition TROOT.h:406
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
R__EXTERN TSystem * gSystem
Definition TSystem.h:555
#define gPad
static std::pair< double, double > matchPrecision(const std::pair< double, double > &in)
Definition xRooFit.cxx:1863
std::shared_ptr< const RooAbsCollection > coords
Definition xRooNLLVar.h:247
xValueWithError limit(const char *type="cls", double nSigma=std::numeric_limits< double >::quiet_NaN()) const
bool AddModel(const xRooNode &pdf, const char *validity="")
std::shared_ptr< TMultiGraph > graphs(const char *opt)
xValueWithError findlimit(const char *opt, double relUncert=std::numeric_limits< double >::infinity(), unsigned int maxTries=20)
int AddPoints(const char *parName, size_t nPoints, double low, double high)
static std::pair< double, double > GetLimit(const TGraph &pValues, double target=std::numeric_limits< double >::quiet_NaN())
std::shared_ptr< TGraphErrors > graph(const char *opt) const
void Print(Option_t *opt="") const override
Print TNamed name and title.
xRooHypoSpace(const char *name="", const char *title="")
int scan(const char *type, size_t nPoints, double low=std::numeric_limits< double >::quiet_NaN(), double high=std::numeric_limits< double >::quiet_NaN(), const std::vector< double > &nSigmas={0, 1, 2, -1, -2, std::numeric_limits< double >::quiet_NaN()}, double relUncert=0.1)
std::shared_ptr< xRooNode > pdf(const RooAbsCollection &parValues) const
void Draw(Option_t *opt="") override
Default Draw method for all objects.
std::map< std::string, std::pair< double, double > > limits(const char *opt="cls", const std::vector< double > &nSigmas={0, 1, 2, -1, -2, std::numeric_limits< double >::quiet_NaN()}, double relUncert=std::numeric_limits< double >::infinity())
std::shared_ptr< RooAbsPdf > pdf() const
Definition xRooNLLVar.h:429
std::shared_ptr< RooArgSet > pars(bool stripGlobalObs=true) const
The xRooNode class is designed to wrap over a TObject and provide functionality to aid with interacti...
Definition xRooNode.h:51
static InteractiveObject * gIntObj
Definition xRooNode.h:465
void Draw(Option_t *opt="") override
Default Draw method for all objects.
xRooNode pars() const
List of parameters (non-observables) of this node.
void setAttribute(const Text_t *name, bool value=true)
Set (default) or clear a named boolean attribute of this object.
Abstract container object that can hold multiple RooAbsArg objects.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void sort(bool reverse=false)
Sort collection using std::sort and name comparison.
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
void setConstant(bool value=true)
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Represents a constant real-valued object.
Definition RooConstVar.h:23
Container class to hold unbinned data.
Definition RooDataSet.h:57
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
int ArraySize() const
number of entries in the results array
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results
RooArgSet * GetParameters() const override
return a cloned list with the parameter of interest
Fill Area Attributes class.
Definition TAttFill.h:19
Line Attributes class.
Definition TAttLine.h:18
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:42
Marker Attributes class.
Definition TAttMarker.h:19
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition TAttMarker.h:38
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition TAttMarker.h:45
void ResetUpdated()
Definition TCanvas.h:193
static TCanvas * MakeDefCanvas()
Static function to build a default canvas.
Definition TCanvas.cxx:1514
void Paint(Option_t *option="") override
Paint canvas.
Definition TCanvas.cxx:1541
void Update() override
Update canvas pad buffers.
Definition TCanvas.cxx:2483
static TClass * GetClass(const char *name, Bool_t load=kTRUE, Bool_t silent=kFALSE)
Static method returning pointer to TClass of the specified class name.
Definition TClass.cxx:2968
Describe directory structure in memory.
Definition TDirectory.h:45
virtual Bool_t cd()
Change current directory to "this" directory.
virtual Int_t GetValue(const char *name, Int_t dflt) const
Returns the integer value for a resource.
Definition TEnv.cxx:491
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition TFile.cxx:4067
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 Double_t GetPointX(Int_t i) const
Get x value for point i.
Definition TGraph.cxx:1522
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2319
static Bool_t CompareX(const TGraph *gr, Int_t left, Int_t right)
Return kTRUE if fX[left] > fX[right]. Can be used by Sort.
Definition TGraph.cxx:683
Int_t GetN() const
Definition TGraph.h:131
virtual void Sort(Bool_t(*greater)(const TGraph *, Int_t, Int_t)=&TGraph::CompareX, Bool_t ascending=kTRUE, Int_t low=0, Int_t high=-1111)
Sorts the points of this TGraph using in-place quicksort (see e.g.
Definition TGraph.cxx:2463
TList * GetListOfFunctions() const
Definition TGraph.h:125
virtual Int_t RemovePoint()
Delete point close to the mouse position Returns index of removed point (or -1 if nothing was changed...
Definition TGraph.cxx:2015
void SetName(const char *name="") override
Set graph name.
Definition TGraph.cxx:2358
virtual Double_t GetPointY(Int_t i) const
Get y value for point i.
Definition TGraph.cxx:1533
virtual void SetPointY(Int_t i, Double_t y)
Set y value for point i.
Definition TGraph.cxx:2351
virtual void SetEditable(Bool_t editable=kTRUE)
if editable=kFALSE, the graph cannot be modified with the mouse by default a TGraph is editable
Definition TGraph.cxx:2278
Book space in a file, create I/O buffers, to fill them, (un)compress them.
Definition TKey.h:28
virtual const char * GetClassName() const
Definition TKey.h:75
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
TObject * FindObject(const char *name) const override
Find an object in this list using its name.
Definition TList.cxx:576
A TMemFile is like a normal TFile except that it reads and writes only from memory.
Definition TMemFile.h:19
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition TMultiGraph.h:34
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:48
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition TNamed.cxx:140
virtual void SetNameTitle(const char *name, const char *title)
Set all the TNamed parameters (name and title).
Definition TNamed.cxx:154
Mother of all ROOT objects.
Definition TObject.h:41
virtual const char * GetName() const
Returns name of object.
Definition TObject.cxx:439
virtual void AppendPad(Option_t *option="")
Append graphics object to current pad.
Definition TObject.cxx:184
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:780
The most important graphics class in the ROOT system.
Definition TPad.h:28
Bool_t Connect(const char *signal, const char *receiver_class, void *receiver, const char *slot)
Non-static method is used to connect from the signal of this object to the receiver slot.
Definition TQObject.cxx:869
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.
Definition TPRegexp.cxx:975
Basic string class.
Definition TString.h:139
Ssiz_t Length() const
Definition TString.h:417
void ToLower()
Change string to lower-case.
Definition TString.cxx:1182
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
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition TString.h:704
Bool_t BeginsWith(const char *s, ECaseCompare cmp=kExact) const
Definition TString.h:623
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
Float_t GetPadRightMargin() const
Definition TStyle.h:212
Float_t GetPadTopMargin() const
Definition TStyle.h:210
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
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition TVirtualPad.h:51
virtual void SetGrid(Int_t valuex=1, Int_t valuey=1)=0
virtual TVirtualPad * cd(Int_t subpadnumber=0)=0
virtual TVirtualPad * GetPad(Int_t subpadnumber) const =0
virtual TCanvas * GetCanvas() const =0
void Clear(Option_t *option="") override=0
TLine * line
double gaussian_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
TGraphErrors * gr
Definition legend1.C:25
leg
Definition legend1.C:34
return c2
Definition legend2.C:14
Definition graph.py:1
#define END_XROOFIT_NAMESPACE
Definition Config.h:25
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4
BEGIN_XROOFIT_NAMESPACE