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 (high < low || (high == low && nPoints != 1)) {
297 // take from parameter
298 low = p->getMin("scan");
299 high = p->getMax("scan");
300 }
301 if (!std::isnan(low) && !std::isnan(high) && !(std::isinf(low) && std::isinf(high))) {
302 p->setRange("scan", low, high);
303 }
304 if (p->hasRange("scan")) {
305 ::Info("xRooHypoSpace::scan", "Using %s scan range: %g - %g", p->GetName(), p->getMin("scan"), p->getMax("scan"));
306 }
307
308 bool doObs = false;
309 for (auto nSigma : nSigmas) {
310 if (std::isnan(nSigma)) {
311 doObs = true;
312 break;
313 }
314 }
315
316 if (fNlls.empty()) {
317 // this happens when loaded hypoSpace from a hypoSpaceInverterResult
318 // set relUncert to infinity so that we don't test any new points
319 relUncert = std::numeric_limits<double>::infinity(); // no NLL available so just get whatever limit we can
320
321 // if any of the defined points are 'expected' data don't do obs
322 // for(auto& hp : *this) {
323 // if(hp.isExpected) {
324 // doObs = false; break;
325 // }
326 // }
327 } else {
328#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 26, 00)
329// bool allGen = true;
330// for (auto &[pdf, nll]: fNlls) {
331// auto _d = dynamic_cast<RooDataSet *>(nll->data());
332// if (!_d || !_d->weightVar() || !_d->weightVar()->getStringAttribute("fitResult") ||
333// !_d->weightVar()->getAttribute("expected")) {
334// allGen = false;
335// break;
336// }
337// }
338// if (allGen)
339// doObs = false;
340#endif
341 }
342
343 // create a fitDatabase if required
344 TDirectory *origDir = gDirectory;
345 if (!gDirectory || !gDirectory->IsWritable()) {
346 // locate a TMemFile in the open list of files and move to that
347 // or create one if cannot find
348 for (auto file : *gROOT->GetListOfFiles()) {
349 if (auto f = dynamic_cast<TMemFile *>(file)) {
350 f->cd();
351 break;
352 }
353 }
354 if (!gDirectory || !gDirectory->IsWritable()) {
355 new TMemFile("fitDatabase", "RECREATE");
356 }
357 }
358
359 int out = 0;
360
361 if (nPoints == 0) {
362 // automatic scan
363 if (sType.Contains("cls")) {
364 for (double nSigma : nSigmas) {
365 xValueWithError res(std::make_pair(0., 0.));
366 if (std::isnan(nSigma)) {
367 if (!doObs)
368 continue;
369 res = findlimit(TString::Format("%s obs", sType.Data()), relUncert);
370 } else {
371 res =
372 findlimit(TString::Format("%s exp%s%d", sType.Data(), nSigma > 0 ? "+" : "", int(nSigma)), relUncert);
373 }
374 if (std::isnan(res.first) || std::isnan(res.second)) {
375 out = 1;
376 } else if (std::isinf(res.second)) {
377 out = 2;
378 }
379 }
380 } else {
381 throw std::runtime_error(TString::Format("Automatic scanning not yet supported for %s", type));
382 }
383 } else {
384 // add the required points and then compute the required value
385 if (nPoints == 1) {
386 AddPoint(TString::Format("%s=%g", poi().first()->GetName(), (high + low) / 2.));
387 graphs(sType); // triggers computation
388 } else {
389 double step = (high - low) / (nPoints - 1);
390 for (size_t i = 0; i < nPoints; i++) {
391 AddPoint(TString::Format("%s=%g", poi().first()->GetName(), low + step * i));
392 graphs(sType); // triggers computation
393 }
394 }
395 }
396
397 if (origDir)
398 origDir->cd();
399
400 return out;
401}
402
403std::map<std::string, xRooNLLVar::xValueWithError>
404xRooNLLVar::xRooHypoSpace::limits(const char *opt, const std::vector<double> &nSigmas, double relUncert)
405{
406
407 if (fNlls.empty()) {
408 // this happens when loaded hypoSpace from a hypoSpaceInverterResult
409 // set relUncert to infinity so that we don't test any new points
410 relUncert = std::numeric_limits<double>::infinity(); // no NLL available so just get whatever limit we can
411 }
412
413 scan(opt, nSigmas, relUncert);
414
415 std::map<std::string, xRooNLLVar::xValueWithError> out;
416 for (auto nSigma : nSigmas) {
417 auto lim = limit(opt, nSigma);
418 if (lim.second < 0)
419 lim.second = -lim.second; // make errors positive for this method
420 out[std::isnan(nSigma) ? "obs" : TString::Format("%d", int(nSigma)).Data()] = xRooFit::matchPrecision(lim);
421 }
422 return out;
423}
424
426{
427 // move to given coords, if any ... will mark them const too
428 std::unique_ptr<RooAbsCollection, std::function<void(RooAbsCollection *)>> _snap(fPars->snapshot(),
429 [&](RooAbsCollection *c) {
430 *fPars = *c;
431 delete c;
432 });
433 TStringToken pattern(coords, ";");
434 while (pattern.NextToken()) {
435 TString s = pattern;
436 // split by "=" sign
437 auto _idx = s.Index('=');
438 if (_idx == -1)
439 continue;
440 TString _name = s(0, _idx);
441 TString _val = s(_idx + 1, s.Length());
442 auto _v = dynamic_cast<RooRealVar *>(fPars->find(_name));
443 if (!_v)
444 continue;
445
446 if (_val.IsFloat()) {
447 _v->setConstant();
448 _v->setVal(_val.Atof());
449 }
450 }
451
452 auto _pdf = pdf();
453
454 if (!_pdf)
455 throw std::runtime_error("no model at coordinates");
456
457 // if (std::unique_ptr<RooAbsCollection>(fPars->selectByAttrib("poi", true))->size() == 0) {
458 // throw std::runtime_error(
459 // "No pars designated as POI - set with pars()->find(<parName>)->setAttribute(\"poi\",true)");
460 // }
461
462 if (fNlls.find(_pdf) == fNlls.end()) {
463 fNlls[_pdf] = std::make_shared<xRooNLLVar>(_pdf->nll("" /*TODO:allow change dataset name and nll opts*/, {}));
464 }
465
466 xRooHypoPoint out;
467
468 out.nllVar = fNlls[_pdf];
469 out.fData = fNlls[_pdf]->getData();
470 out.isExpected = dynamic_cast<RooDataSet *>(out.fData.first.get()) &&
471 dynamic_cast<RooDataSet *>(out.fData.first.get())->weightVar()->getAttribute("expected");
472 // TODO: need to access the genfit of the data and add that to the point, somehow ...
473
474 out.coords.reset(fPars->snapshot()); // should already have altVal prop on poi, and poi labelled
475 // ensure all poi are marked const ... required by xRooHypoPoint behaviour
476 out.poi().setAttribAll("Constant");
477 // and now remove anything that's marked floating
478
479 // do to bug in remove have to ensure not using the hash map otherwise will be doing an invalid read after the
480 // deletion of the owned pars
481 const_cast<RooAbsCollection *>(out.coords.get())->useHashMapForFind(false);
482 const_cast<RooAbsCollection *>(out.coords.get())
483 ->remove(*std::unique_ptr<RooAbsCollection>(out.coords->selectByAttrib("Constant", false)), true, true);
484
485 double value = out.fNullVal();
486 double alt_value = out.fAltVal();
487
488 auto _type = fTestStatType;
489 if (_type == xRooFit::Asymptotics::Unknown) {
490 // decide based on values
491 if (std::isnan(alt_value)) {
493 } else if (value >= alt_value) {
495 } else {
497 }
498 }
499
500 out.fPllType = _type;
501
502 // look for a matching point
503 for (auto &p : *this) {
504 if (p.nllVar != out.nllVar)
505 continue;
506 if (p.fData != out.fData)
507 continue;
508 if (!p.alt_poi().equals(out.alt_poi()))
509 continue;
510 bool match = true;
511 for (auto c : p.alt_poi()) {
512 if (auto v = dynamic_cast<RooAbsReal *>(c);
513 v && std::abs(v->getVal() - out.alt_poi().getRealValue(v->GetName())) > 1e-12) {
514 match = false;
515 break;
516 } else if (auto cat = dynamic_cast<RooAbsCategory *>(c);
517 cat && cat->getCurrentIndex() ==
518 out.alt_poi().getCatIndex(cat->GetName(), std::numeric_limits<int>().max())) {
519 match = false;
520 break;
521 }
522 }
523 if (!match)
524 continue;
525 if (!p.coords->equals(*out.coords))
526 continue;
527 for (auto c : *p.coords) {
528 if (c->getAttribute("poi")) {
529 continue; // check poi below
530 }
531 if (auto v = dynamic_cast<RooAbsReal *>(c);
532 v && std::abs(v->getVal() - out.coords->getRealValue(v->GetName())) > 1e-12) {
533 match = false;
534 break;
535 } else if (auto cat = dynamic_cast<RooAbsCategory *>(c);
536 cat && cat->getCurrentIndex() ==
537 out.alt_poi().getCatIndex(cat->GetName(), std::numeric_limits<int>().max())) {
538 match = false;
539 break;
540 }
541 }
542 if (!match)
543 continue;
544 // if reached here we can copy over the asimov dataset to save re-generating it
545 // first copy over cfit_alt (if its there) because that is also the same since coords and alt_poi values same
546 if (auto cfit = p.cfit_alt(true)) {
547 out.fAlt_cfit = cfit;
548 }
549 if (p.asimov(true) && p.asimov(true)->fData.first && (!out.asimov(true) || !out.asimov(true)->fData.first)) {
550 out.asimov()->fData = p.asimov(true)->fData;
551 }
552 if (!p.poi().equals(out.poi()))
553 continue;
554 for (auto c : p.poi()) {
555 if (auto v = dynamic_cast<RooAbsReal *>(c);
556 v && std::abs(v->getVal() - out.poi().getRealValue(v->GetName())) > 1e-12) {
557 match = false;
558 break;
559 }
560 }
561 if (match) {
562 // found a duplicate point, return that!
563 return p;
564 }
565 }
566
567 std::string coordString;
568 for (auto a : axes()) {
569 coordString += TString::Format("%s=%g", a->GetName(), out.coords->getRealValue(a->GetName()));
570 coordString += ",";
571 }
572 coordString.erase(coordString.end() - 1);
573
574 ::Info("xRooHypoSpace::AddPoint", "Added new point @ %s", coordString.c_str());
575 return emplace_back(out);
576}
577
578bool xRooNLLVar::xRooHypoSpace::AddModel(const xRooNode &_pdf, const char *validity)
579{
580
581 if (!_pdf.get<RooAbsPdf>()) {
582 throw std::runtime_error("Not a pdf");
583 }
584
585 auto pars = _pdf.pars().argList();
586
587 // replace any pars with validity pars and add new pars
588 auto vpars = toArgs(validity);
589 pars.replace(vpars);
590 vpars.remove(pars, true, true);
591 pars.add(vpars);
592
593 if (auto existing = pdf(pars)) {
594 throw std::runtime_error(std::string("Clashing model: ") + existing->GetName());
595 }
596
597 auto myPars = std::shared_ptr<RooArgList>(dynamic_cast<RooArgList *>(pars.snapshot()));
598 myPars->sort();
599
600 pars.remove(*fPars, true, true);
601
602 fPars->addClone(pars);
603
604 fPdfs.insert(std::make_pair(myPars, std::make_shared<xRooNode>(_pdf)));
605
606 return true;
607}
608
610{
611 // determine which pars are the minimal set to distinguish all points in the space
612 RooArgList out;
613 out.setName("axes");
614
615 out.add(*std::unique_ptr<RooAbsCollection>(
616 fPars->selectByAttrib("axis", true))); // start with any pars explicitly designated as axes
617
618 bool clash;
619 do {
620 clash = false;
621
622 std::set<std::vector<double>> coords;
623 for (auto &p : *this) {
624 std::vector<double> p_coords;
625 for (auto o : out) {
626 auto _v = dynamic_cast<RooRealVar *>(p.coords->find(o->GetName()));
627 p_coords.push_back(
628 (_v && _v->isConstant())
629 ? _v->getVal()
630 : std::numeric_limits<double>::infinity()); // non-const coords are treating as non-existent
631 // p_coords.push_back(p.coords->getRealValue(o->GetName(), std::numeric_limits<double>::quiet_NaN()));
632 }
633 if (coords.find(p_coords) != coords.end()) {
634 clash = true;
635 break;
636 }
637 coords.insert(p_coords);
638 }
639
640 if (clash) {
641 // add next best coordinate
642 std::map<std::string, std::unordered_set<double>> values;
643 for (auto &par : *pars()) {
644 if (out.find(*par))
645 continue;
646 for (auto p : *this) {
647 auto _v = dynamic_cast<RooRealVar *>(p.coords->find(par->GetName()));
648 values[par->GetName()].insert(
649 (_v && _v->isConstant())
650 ? _v->getVal()
651 : std::numeric_limits<double>::infinity()); // non-const coords are treating as non-existent
652 // values[par->GetName()].insert(
653 // p.coords->getRealValue(par->GetName(), std::numeric_limits<double>::quiet_NaN()));
654 }
655 }
656
657 std::string bestVar;
658 size_t maxDiff = 0;
659 bool isPOI = false;
660 for (auto &[k, v] : values) {
661 if (v.size() > maxDiff || (v.size() == maxDiff && !isPOI && pars()->find(k.c_str())->getAttribute("poi"))) {
662 bestVar = k;
663 isPOI = pars()->find(k.c_str())->getAttribute("poi");
664 maxDiff = std::max(maxDiff, v.size());
665 }
666 }
667 if (bestVar.empty()) {
668 break;
669 }
670
671 out.add(*pars()->find(bestVar.c_str()));
672 }
673 } while (clash);
674
675 // ensure poi are at the end
676 std::unique_ptr<RooAbsCollection> poi(out.selectByAttrib("poi", true));
677 out.remove(*poi);
678 out.add(*poi);
679
680 return out;
681}
682
684{
685 RooArgList out;
686 out.setName("poi");
687 out.add(*std::unique_ptr<RooAbsCollection>(pars()->selectByAttrib("poi", true)));
688 return out;
689}
690
692{
693
694 if (!gDirectory)
695 return;
696 auto dir = gDirectory->GetDirectory(apath);
697 if (!dir) {
698 // try open file first
699 TString s(apath);
700 auto f = TFile::Open(s.Contains(":") ? TString(s(0, s.Index(":"))) : s);
701 if (f) {
702 if (!s.Contains(":"))
703 s += ":";
704 dir = gDirectory->GetDirectory(s);
705 if (dir) {
706 LoadFits(s);
707 return;
708 }
709 }
710 if (!dir) {
711 Error("LoadFits", "Path not found %s", apath);
712 return;
713 }
714 }
715
716 // assume for now all fits in given path will have the same pars
717 // so can just look at the float and const pars of first fit result to get all of them
718 // tuple is: parName, parValue, parAltValue (blank if nan)
719 // key represents the ufit values, value represents the sets of poi for the available cfits (subfits of the ufit)
720
721 std::map<std::set<std::tuple<std::string, double, std::string>>, std::set<std::set<std::string>>> cfits;
722 std::set<std::string> allpois;
723
724 int nFits = 0;
725 std::function<void(TDirectory *)> processDir;
726 processDir = [&](TDirectory *_dir) {
727 std::cout << "Processing " << _dir->GetName() << std::endl;
728 if (auto keys = _dir->GetListOfKeys(); keys) {
729 // first check if dir doesn't contain any RooLinkedList ... this identifies it as not an nll dir
730 // so treat any sub-dirs as new nll
731 bool isNllDir = false;
732 for (auto &&k : *keys) {
733 TKey *key = dynamic_cast<TKey *>(k);
734 if (strcmp(key->GetClassName(), "RooLinkedList") == 0) {
735 isNllDir = true;
736 break;
737 }
738 }
739
740 for (auto &&k : *keys) {
741 if (auto subdir = _dir->GetDirectory(k->GetName()); subdir) {
742 if (!isNllDir) {
743 LoadFits(subdir->GetPath());
744 } else {
745 processDir(subdir);
746 }
747 continue;
748 }
749 auto cl = TClass::GetClass((static_cast<TKey *>(k))->GetClassName());
750 if (cl->InheritsFrom("RooFitResult")) {
751 if (auto cachedFit = _dir->Get<RooFitResult>(k->GetName()); cachedFit) {
752 nFits++;
753 if (nFits == 1) {
754 // for first fit add any missing float pars
755 std::unique_ptr<RooAbsCollection> snap(cachedFit->floatParsFinal().snapshot());
756 snap->remove(*fPars, true, true);
757 fPars->addClone(*snap);
758 // add also the non-string const pars
759 for (auto &p : cachedFit->constPars()) {
760 if (p->getAttribute("global"))
761 continue; // don't consider globals
762 auto v = dynamic_cast<RooAbsReal *>(p);
763 if (!v) {
764 continue;
765 };
766 if (!fPars->contains(*v))
767 fPars->addClone(*v);
768 }
769 }
770 // get names of all the floats
771 std::set<std::string> floatPars;
772 for (auto &p : cachedFit->floatParsFinal())
773 floatPars.insert(p->GetName());
774 // see if
775
776 // build a set of the const par values
777 std::set<std::tuple<std::string, double, std::string>> constPars;
778 for (auto &p : cachedFit->constPars()) {
779 if (p->getAttribute("global"))
780 continue; // don't consider globals when looking for cfits
781 auto v = dynamic_cast<RooAbsReal *>(p);
782 if (!v) {
783 continue;
784 };
785 constPars.insert(
786 std::make_tuple(v->GetName(), v->getVal(),
787 v->getStringAttribute("altVal") ? v->getStringAttribute("altVal") : ""));
788 }
789 // now see if this is a subset of any existing cfit ...
790 for (auto &&[key, value] : cfits) {
791 if (constPars == key)
792 continue; // ignore cases where we already recorded this list of constPars
793 if (std::includes(constPars.begin(), constPars.end(), key.begin(), key.end())) {
794 // usual case ... cachedFit has more constPars than one of the fits we have already encountered
795 // (the ufit)
796 // => cachedFit is a cfit of key fr ...
797 std::set<std::string> pois;
798 for (auto &&par : constPars) {
799 if (key.find(par) == key.end()) {
800 pois.insert(std::get<0>(par));
801 allpois.insert(std::get<0>(par));
802 }
803 }
804 if (!pois.empty()) {
805 cfits[constPars].insert(pois);
806 // std::cout << cachedFit->GetName() << " ";
807 // for(auto ff: constPars) std::cout << ff.first << "=" <<
808 // ff.second << " "; std::cout << std::endl;
809 }
810 }
811 /* FOR NOW we will skip cases where we encounter the cfit before the ufit - usually should eval the
812 ufit first
813 * else if (std::includes(key.begin(), key.end(), constPars.begin(), constPars.end())) {
814 // constPars are subset of key
815 // => key is a ufit of the cachedFit
816 // add all par names of key that aren't in constPars ... these are the poi
817 std::set<std::string> pois;
818 for (auto &&par: key) {
819 if (constPars.find(par) == constPars.end()) {
820 pois.insert(std::get<0>(par));
821 allpois.insert(std::get<0>(par));
822 }
823 }
824 if (!pois.empty()) {
825 std::cout << "found cfit BEFORE ufit??" << std::endl;
826 value.insert(pois);
827 }
828 } */
829 }
830 // ensure that this combination of constPars has entry in map,
831 // even if it doesn't end up with any poi identified from cfits to it
832 cfits[constPars];
833 delete cachedFit;
834 }
835 }
836 }
837 }
838 };
839 processDir(dir);
840 ::Info("xRooHypoSpace::xRooHypoSpace", "%s - Loaded %d fits", apath, nFits);
841
842 if (allpois.size() == 1) {
843 ::Info("xRooHypoSpace::xRooHypoSpace", "Detected POI: %s", allpois.begin()->c_str());
844
845 auto nll = std::make_shared<xRooNLLVar>(nullptr, nullptr);
846 auto dummyNll = std::make_shared<RooRealVar>(apath, "Dummy NLL", std::numeric_limits<double>::quiet_NaN());
847 nll->std::shared_ptr<RooAbsReal>::operator=(dummyNll);
848 dummyNll->setAttribute("readOnly");
849 // add pars as 'servers' on the dummy NLL
850 if (fPars) {
851 for (auto &&p : *fPars) {
852 dummyNll->addServer(
853 *p); // this is ok provided fPars (i.e. hypoSpace) stays alive as long as the hypoPoint ...
854 }
855 // flag poi
856 for (auto &p : allpois) {
857 fPars->find(p.c_str())->setAttribute("poi", true);
858 }
859 }
860 nll->reinitialize(); // triggers filling of par lists etc
861
862 for (auto &&[key, value] : cfits) {
863 if (value.find(allpois) != value.end()) {
864 // get the value of the poi in the key set
865 auto _coords = std::make_shared<RooArgSet>();
866 for (auto &k : key) {
867 auto v = _coords->addClone(RooRealVar(std::get<0>(k).c_str(), std::get<0>(k).c_str(), std::get<1>(k)));
868 v->setAttribute("poi", allpois.find(std::get<0>(k)) != allpois.end());
869 if (!std::get<2>(k).empty()) {
870 v->setStringAttribute("altVal", std::get<2>(k).c_str());
871 }
872 }
874 // hp.fPOIName = allpois.begin()->c_str();
875 // hp.fNullVal = _coords->getRealValue(hp.fPOIName.c_str());
876 hp.coords = _coords;
877 hp.nllVar = nll;
878
879 // auto altVal =
880 // hp.null_cfit()->constPars().find(hp.fPOIName.c_str())->getStringAttribute("altVal");
881 // if(altVal) hp.fAltVal = TString(altVal).Atof();
882 // else hp.fAltVal = std::numeric_limits<double>::quiet_NaN();
883
884 // decide based on values
885 if (std::isnan(hp.fAltVal())) {
887 } else if (hp.fNullVal() >= hp.fAltVal()) {
889 } else {
891 }
892
893 emplace_back(hp);
894 }
895 }
896 } else if (nFits > 0) {
897 std::cout << "possible POI: ";
898 for (auto p : allpois)
899 std::cout << p << ",";
900 std::cout << std::endl;
901 }
902}
903
905{
906
907 auto _axes = axes();
908
909 size_t badFits = 0;
910
911 for (size_t i = 0; i < size(); i++) {
912 std::cout << i << ") ";
913 for (auto a : _axes) {
914 if (a != _axes.first())
915 std::cout << ",";
916 std::cout << a->GetName() << "="
917 << at(i).coords->getRealValue(a->GetName(), std::numeric_limits<double>::quiet_NaN());
918 }
919 std::cout << " status=[ufit:";
920 auto ufit = const_cast<xRooHypoPoint &>(at(i)).ufit(true);
921 if (!ufit) {
922 std::cout << "-";
923 } else {
924 std::cout << ufit->status();
925 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(ufit->status()) == 0);
926 }
927 std::cout << ",cfit_null:";
928 auto cfit = const_cast<xRooHypoPoint &>(at(i)).cfit_null(true);
929 if (!cfit) {
930 std::cout << "-";
931 } else {
932 std::cout << cfit->status();
933 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(cfit->status()) == 0);
934 }
935 std::cout << ",cfit_alt:";
936 auto afit = const_cast<xRooHypoPoint &>(at(i)).cfit_alt(true);
937 if (!afit) {
938 std::cout << "-";
939 } else {
940 std::cout << afit->status();
941 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(afit->status()) == 0);
942 }
943 if (auto asiPoint = const_cast<xRooHypoPoint &>(at(i)).asimov(true)) {
944 std::cout << ",asimov.ufit:";
945 auto asi_ufit = asiPoint->ufit(true);
946 if (!asi_ufit) {
947 std::cout << "-";
948 } else {
949 std::cout << asi_ufit->status();
950 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(asi_ufit->status()) == 0);
951 }
952 std::cout << ",asimov.cfit_null:";
953 auto asi_cfit = asiPoint->cfit_null(true);
954 if (!asi_cfit) {
955 std::cout << "-";
956 } else {
957 std::cout << asi_cfit->status();
958 badFits += (xRooNLLVar::xRooHypoPoint::allowedStatusCodes.count(asi_cfit->status()) == 0);
959 }
960 }
961 std::cout << "]";
962 auto sigma_mu = const_cast<xRooHypoPoint &>(at(i)).sigma_mu(true);
963 if (!std::isnan(sigma_mu.first)) {
964 std::cout << " sigma_mu=" << sigma_mu.first;
965 if (sigma_mu.second)
966 std::cout << " +/- " << sigma_mu.second;
967 }
968 std::cout << std::endl;
969 }
970 std::cout << "--------------------------" << std::endl;
971 std::cout << "Number of bad fits: " << badFits << std::endl;
972}
973
974std::shared_ptr<TGraphErrors> xRooNLLVar::xRooHypoSpace::graph(
975 const char *opt /*, const std::function<void(xRooNLLVar::xRooHypoSpace*)>& progress*/) const
976{
977
978 TString sOpt(opt);
979 sOpt.ToLower();
980
981 bool doCLs = sOpt.Contains("cls");
982 bool readOnly = sOpt.Contains("readonly");
983 bool visualize = sOpt.Contains("visualize") && !readOnly;
984
985 double nSigma =
986 (sOpt.Contains("exp"))
987 ? (TString(sOpt(sOpt.Index("exp") + 3,
988 sOpt.Index(" ", sOpt.Index("exp")) == -1 ? sOpt.Length() : sOpt.Index(" ", sOpt.Index("exp"))))
989 .Atof())
990 : std::numeric_limits<double>::quiet_NaN();
991 bool expBand =
992 !std::isnan(nSigma) && nSigma && !(sOpt(sOpt.Index("exp") + 3) == '+' || sOpt(sOpt.Index("exp") + 3) == '-');
993
994 auto _axes = axes();
995 if (_axes.size() != 1)
996 return nullptr;
997
998 auto out = std::make_shared<TGraphErrors>();
999 out->SetName(GetName());
1000 out->SetEditable(false);
1001 const char *sCL = (doCLs) ? "CLs" : "null";
1002
1003 TString title =
1004 TString::Format("%s;%s;p_{%s}",
1005 (std::isnan(nSigma))
1006 ? "Observed"
1007 : (!nSigma ? "Expected"
1008 : TString::Format("%s%d#sigma Expected",
1009 expBand || !nSigma ? "" : ((nSigma < 0) ? "-" : "+"), int(nSigma))
1010 .Data()),
1011 _axes.at(0)->GetTitle(), sCL);
1012
1013 if (std::isnan(nSigma)) {
1014 out->SetNameTitle(TString::Format("obs_p%s", sCL), title);
1015 out->SetMarkerStyle(20);
1016 out->SetMarkerSize(0.5);
1017 if (sOpt.Contains("ts")) {
1018 out->SetNameTitle("obs_ts", TString::Format("Observed;%s;%s", _axes.at(0)->GetTitle(),
1019 (empty() ? "" : front().tsTitle(true).Data())));
1020 }
1021 } else {
1022 out->SetNameTitle(TString::Format("exp%d_p%s", int(nSigma), sCL), title);
1023 out->SetMarkerStyle(0);
1024 out->SetMarkerSize(0);
1025 out->SetLineStyle(2 + int(nSigma));
1026 if (expBand && nSigma) {
1027 out->SetFillColor((nSigma == 2) ? kYellow : kGreen);
1028 // out->SetFillStyle(3005);
1029 out->SetLineStyle(0);
1030 out->SetLineWidth(0);
1031 auto x = out->Clone("up");
1032 x->SetBit(kCanDelete);
1033 dynamic_cast<TAttFill *>(x)->SetFillStyle(nSigma == 2 ? 3005 : 3004);
1034 dynamic_cast<TAttFill *>(x)->SetFillColor(kBlack);
1035 out->GetListOfFunctions()->Add(x, "F");
1036 x = out->Clone("down");
1037 x->SetBit(kCanDelete);
1038 // dynamic_cast<TAttFill*>(x)->SetFillColor((nSigma==2) ? kYellow : kGreen);
1039 // dynamic_cast<TAttFill*>(x)->SetFillStyle(1001);
1040 out->GetListOfFunctions()->Add(x, "F");
1041 }
1042 if (sOpt.Contains("ts")) {
1043 out->SetNameTitle(TString::Format("exp_ts%d", int(nSigma)),
1044 TString::Format("Expected;%s;%s", _axes.at(0)->GetTitle(), front().tsTitle(true).Data()));
1045 }
1046 }
1047
1048 auto badPoints = [&]() {
1049 auto badPoints2 = dynamic_cast<TGraph *>(out->GetListOfFunctions()->FindObject("badPoints"));
1050 if (!badPoints2) {
1051 badPoints2 = new TGraph;
1052 badPoints2->SetBit(kCanDelete);
1053 badPoints2->SetName("badPoints");
1054 badPoints2->SetMarkerStyle(5);
1055 badPoints2->SetMarkerColor(std::isnan(nSigma) ? kRed : kBlue);
1056 badPoints2->SetMarkerSize(1);
1057 out->GetListOfFunctions()->Add(badPoints2, "P");
1058 }
1059 return badPoints2;
1060 };
1061 int nPointsDown = 0;
1062 bool above = true;
1063 TStopwatch s;
1064 s.Start();
1065 size_t nDone = 0;
1066 for (auto &p : *this) {
1067 if (s.RealTime() > 5) {
1068 if (visualize) {
1069 // draw readonly version of the graph
1070 auto gra = graph(sOpt + " readOnly");
1071 if (gra && gra->GetN()) {
1072 if (!gPad && gROOT->GetSelectedPad())
1073 gROOT->GetSelectedPad()->cd();
1074 if (gPad)
1075 gPad->Clear();
1076 gra->DrawClone(expBand ? "AF" : "ALP")->SetBit(kCanDelete);
1077#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1078 if (auto pad = gROOT->GetSelectedPad()) {
1079 pad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1080 }
1081#endif
1083 }
1084 } else {
1085 ::Info("xRooHypoSpace::graph", "Completed %d/%d points for %s", int(nDone), int(size()), sOpt.Data());
1086 }
1087 s.Start();
1088 } else {
1089 s.Continue();
1090 }
1091 double _x = p.coords->getRealValue(_axes.at(0)->GetName(), std::numeric_limits<double>::quiet_NaN());
1092 auto pval = const_cast<xRooHypoPoint &>(p).getVal(sOpt);
1093 auto idx = out->GetN() - nPointsDown;
1094
1095 if (std::isnan(pval.first)) {
1096 if (p.status() != 0) { // if status is 0 then bad pval is really just absence of fits, not bad fits
1097 badPoints()->SetPoint(badPoints()->GetN(), _x, 0);
1098 }
1099 } else {
1100 out->InsertPointBefore(idx, _x, pval.first);
1101 out->SetPointError(idx, 0, pval.second);
1102 }
1103
1104 if (expBand && nSigma) {
1105 TString sOpt2 = sOpt;
1106 sOpt2.ReplaceAll("exp", "exp-");
1107 pval = const_cast<xRooHypoPoint &>(p).getVal(sOpt2);
1108 if (std::isnan(pval.first)) {
1109 if (p.status() != 0) { // if status is 0 then bad pval is really just absence of fits, not bad fits
1110 badPoints()->SetPoint(badPoints()->GetN(), _x, 0);
1111 }
1112 } else {
1113 out->InsertPointBefore(idx + 1, _x, pval.first);
1114 out->SetPointError(idx + 1, 0, pval.second);
1115 nPointsDown++;
1116 if (out->GetPointY(idx) < pval.first)
1117 above = false; // the -sigma points are actually above +sigma
1118 }
1119 }
1120 nDone++;
1121 }
1122
1123 if (out->GetN() == 0)
1124 return out;
1125
1126 if (!expBand) {
1127 out->Sort();
1128 if (out->GetListOfFunctions()->FindObject("badPoints")) {
1129 // try to interpolate the points
1130 for (int i = 0; i < badPoints()->GetN(); i++) {
1131 badPoints()->SetPointY(i, out->Eval(badPoints()->GetPointX(i)));
1132 }
1133 }
1134 } else {
1135 out->Sort(&TGraph::CompareX, true, 0, out->GetN() - nPointsDown - 1); // sort first half
1136 out->Sort(&TGraph::CompareX, false, out->GetN() - nPointsDown, out->GetN() - 1); // reverse sort second half
1137
1138 // now populate the up and down values
1139 auto up = dynamic_cast<TGraph *>(out->GetListOfFunctions()->FindObject("up"));
1140 auto down = dynamic_cast<TGraph *>(out->GetListOfFunctions()->FindObject("down"));
1141
1142 for (int i = 0; i < out->GetN(); i++) {
1143 if (i < out->GetN() - nPointsDown) {
1144 up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) + out->GetErrorY(i) * (above ? 1. : -1.));
1145 down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.));
1146 } else {
1147 up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.));
1148 down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) + out->GetErrorY(i) * (above ? 1. : -1.));
1149 }
1150 }
1151 }
1152
1153 if (visualize) {
1154 // draw result
1155 if (!gPad && gROOT->GetSelectedPad())
1156 gROOT->GetSelectedPad()->cd();
1157 if (gPad)
1158 gPad->Clear();
1159 out->DrawClone(expBand ? "AF" : "ALP")->SetBit(kCanDelete);
1160#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1161 if (auto pad = gROOT->GetSelectedPad()) {
1162 pad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1163 }
1164#endif
1166 }
1167
1168 return out;
1169}
1170
1171std::shared_ptr<TMultiGraph> xRooNLLVar::xRooHypoSpace::graphs(const char *opt)
1172{
1173 TString sOpt(opt);
1174 sOpt.ToLower();
1175 std::shared_ptr<TMultiGraph> out;
1176 if (sOpt.Contains("pcls") || sOpt.Contains("pnull") || sOpt.Contains("ts")) {
1177
1178 bool visualize = sOpt.Contains("visualize");
1179 sOpt.ReplaceAll("visualize", "");
1180
1181 auto exp2 = graph(sOpt + " exp2");
1182 auto exp1 = graph(sOpt + " exp1");
1183 auto exp = graph(sOpt + " exp");
1184 bool doObs = true;
1185 // for(auto& hp : *this) { if(hp.isExpected) {doObs=false; break;} }
1186 auto obs = (doObs) ? graph(sOpt) : nullptr;
1187
1188 out = std::make_shared<TMultiGraph>(GetName(), GetTitle());
1189 if (exp2 && exp2->GetN() > 1)
1190 out->Add(static_cast<TGraph *>(exp2->Clone()), "FP");
1191 if (exp1 && exp1->GetN() > 1)
1192 out->Add(static_cast<TGraph *>(exp1->Clone()), "FP");
1193 if (exp && exp->GetN() > 1)
1194 out->Add(static_cast<TGraph *>(exp->Clone()), "LP");
1195 if (obs && obs->GetN() > 1)
1196 out->Add(static_cast<TGraph *>(obs->Clone()), "LP");
1197
1198 if (!out->GetListOfGraphs()) {
1199 return nullptr;
1200 }
1201
1202 TGraph *testedPoints = nullptr;
1203 if (sOpt.Contains("pcls")) {
1204 TGraph *line = new TGraph;
1205 line->SetName("alpha");
1206 line->SetLineStyle(2);
1207 line->SetEditable(false);
1208 line->SetPoint(line->GetN(), out->GetHistogram()->GetXaxis()->GetXmin() - 10, 0.05);
1209 testedPoints = new TGraph;
1210 testedPoints->SetName("hypoPoints");
1211 testedPoints->SetEditable(false);
1212 testedPoints->SetMarkerStyle(24);
1213 testedPoints->SetMarkerSize(0.4); // use line to indicate tested points
1214 if (exp) {
1215 for (int i = 0; i < exp->GetN(); i++) {
1216 testedPoints->SetPoint(testedPoints->GetN(), exp->GetPointX(i), 0.05);
1217 }
1218 }
1219 line->SetPoint(line->GetN(), out->GetHistogram()->GetXaxis()->GetXmax() + 10, 0.05);
1221 out->GetListOfFunctions()->Add(line, "L");
1222 }
1223 if (exp) {
1224 out->GetHistogram()->GetXaxis()->SetTitle(exp->GetHistogram()->GetXaxis()->GetTitle());
1225 out->GetHistogram()->GetYaxis()->SetTitle(exp->GetHistogram()->GetYaxis()->GetTitle());
1226 }
1227 auto leg = new TLegend(1. - gStyle->GetPadRightMargin() - 0.3, 1. - gStyle->GetPadTopMargin() - 0.35,
1228 1. - gStyle->GetPadRightMargin() - 0.05, 1. - gStyle->GetPadTopMargin() - 0.05);
1229 leg->SetName("legend");
1230 leg->SetBit(kCanDelete);
1231
1232 out->GetListOfFunctions()->Add(leg);
1233 // out->GetListOfFunctions()->Add(out->GetHistogram()->Clone(".axis"),"sameaxis"); // redraw axis
1234
1235 for (auto g : *out->GetListOfGraphs()) {
1236 if (auto o = dynamic_cast<TGraph *>(g)->GetListOfFunctions()->FindObject("down")) {
1237 leg->AddEntry(o, "", "F");
1238 } else {
1239 leg->AddEntry(g, "", "LPE");
1240 }
1241 }
1242
1243 if (sOpt.Contains("pcls")) {
1244 // add current limit estimates to legend
1245 if (exp2 && exp2->GetN() > 1) {
1246 auto l = xRooFit::matchPrecision(GetLimit(*graph(sOpt + "exp-2")));
1247 leg->AddEntry((TObject *)nullptr, TString::Format("-2#sigma: %g +/- %g", l.first, l.second), "");
1248 }
1249 if (exp1 && exp1->GetN() > 1) {
1250 auto l = xRooFit::matchPrecision(GetLimit(*graph(sOpt + "exp-1")));
1251 leg->AddEntry((TObject *)nullptr, TString::Format("-1#sigma: %g +/- %g", l.first, l.second), "");
1252 }
1253 if (exp && exp->GetN() > 1) {
1254 auto l = xRooFit::matchPrecision(GetLimit(*exp));
1255 leg->AddEntry((TObject *)nullptr, TString::Format("0#sigma: %g +/- %g", l.first, l.second), "");
1256 }
1257 if (exp1 && exp1->GetN() > 1) {
1258 auto l = xRooFit::matchPrecision(GetLimit(*graph(sOpt + "exp+1")));
1259 leg->AddEntry((TObject *)nullptr, TString::Format("+1#sigma: %g +/- %g", l.first, l.second), "");
1260 }
1261 if (exp2 && exp2->GetN() > 1) {
1262 auto l = xRooFit::matchPrecision(GetLimit(*graph(sOpt + "exp+2")));
1263 leg->AddEntry((TObject *)nullptr, TString::Format("+2#sigma: %g +/- %g", l.first, l.second), "");
1264 }
1265 if (obs && obs->GetN() > 1) {
1266 auto l = xRooFit::matchPrecision(GetLimit(*obs));
1267 leg->AddEntry((TObject *)nullptr, TString::Format("Observed: %g +/- %g", l.first, l.second), "");
1268 }
1269 }
1270 if (testedPoints)
1271 out->Add(testedPoints, "P");
1272
1273 if (visualize) {
1274 if (!gPad && gROOT->GetSelectedPad())
1275 gROOT->GetSelectedPad()->cd();
1276 if (gPad)
1277 gPad->Clear();
1278 auto gra2 = static_cast<TMultiGraph *>(out->DrawClone("A"));
1279 gra2->SetBit(kCanDelete);
1280 if (sOpt.Contains("pcls") || sOpt.Contains("pnull")) {
1281 gra2->GetHistogram()->SetMinimum(1e-6);
1282 }
1283 if (gPad) {
1284 gPad->RedrawAxis();
1285 gPad->GetCanvas()->Paint();
1286 gPad->GetCanvas()->Update();
1287#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1288 gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1289#endif
1290 }
1292 }
1293 }
1294
1295 return out;
1296}
1297
1299{
1300
1301 if (std::isnan(target)) {
1302 target = 1. - gEnv->GetValue("xRooHypoSpace.CL", 95.) / 100.;
1303 }
1304
1305 auto gr = std::make_shared<TGraph>(pValues);
1306 // remove any nan points and duplicates
1307 int i = 0;
1308 std::set<double> existingX;
1309 while (i < gr->GetN()) {
1310 if (std::isnan(gr->GetPointY(i))) {
1311 gr->RemovePoint(i);
1312 } else if (existingX.find(gr->GetPointX(i)) != existingX.end()) {
1313 gr->RemovePoint(i);
1314 } else {
1315 existingX.insert(gr->GetPointX(i));
1316 // convert to log ....
1317 gr->SetPointY(i, log(std::max(gr->GetPointY(i), 1e-10)));
1318 i++;
1319 }
1320 }
1321
1322 gr->Sort();
1323
1324 // simple linear extrapolation to critical value ... return nan if problem
1325 if (gr->GetN() < 2) {
1326 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1327 }
1328
1329 double alpha = log(target);
1330
1331 bool above = gr->GetPointY(0) > alpha;
1332 for (int ii = 1; ii < gr->GetN(); ii++) {
1333 if ((above && (gr->GetPointY(ii) <= alpha)) || (!above && (gr->GetPointY(ii) >= alpha))) {
1334 // found the limit ... return linearly extrapolated point
1335 double lim = gr->GetPointX(ii - 1) + (gr->GetPointX(ii) - gr->GetPointX(ii - 1)) *
1336 (alpha - gr->GetPointY(ii - 1)) /
1337 (gr->GetPointY(ii) - gr->GetPointY(ii - 1));
1338 // use points either side as error
1339 double err = std::max(lim - gr->GetPointX(ii - 1), gr->GetPointX(ii) - lim);
1340 // give err negative sign to indicate if error due to negative side
1341 if ((lim - gr->GetPointX(ii - 1)) > (gr->GetPointX(ii) - lim))
1342 err *= -1;
1343 return std::pair(lim, err);
1344 }
1345 }
1346 // if reach here need to extrapolate ...
1347 if ((above && gr->GetPointY(gr->GetN() - 1) <= gr->GetPointY(0)) ||
1348 (!above && gr->GetPointY(gr->GetN() - 1) >= gr->GetPointY(0))) {
1349 // extrapolating above based on last two points
1350 // in fact, if 2nd last point is a p=1 (log(p)=0) then go back
1351 int offset = 2;
1352 while (offset < gr->GetN() && gr->GetPointY(gr->GetN() - offset) == 0)
1353 offset++;
1354 double x1 = gr->GetPointX(gr->GetN() - offset);
1355 double y1 = gr->GetPointY(gr->GetN() - offset);
1356 double m = (gr->GetPointY(gr->GetN() - 1) - y1) / (gr->GetPointX(gr->GetN() - 1) - x1);
1357 if (m == 0.)
1358 return std::pair(std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity());
1359 return std::pair((alpha - y1) / m + x1, std::numeric_limits<double>::infinity());
1360 } else {
1361 // extrapolating below based on first two points
1362 double x1 = gr->GetPointX(0);
1363 double y1 = gr->GetPointY(0);
1364 double m = (gr->GetPointY(1) - y1) / (gr->GetPointX(1) - x1);
1365 if (m == 0.)
1366 return std::pair(-std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity());
1367 return std::pair((alpha - y1) / m + x1, -std::numeric_limits<double>::infinity());
1368 }
1369}
1370
1372{
1373 TString sOpt = TString::Format("p%s", type);
1374 if (std::isnan(nSigma)) {
1375 sOpt += "obs";
1376 } else {
1377 sOpt += TString::Format("exp%s%d", nSigma > 0 ? "+" : "", int(nSigma));
1378 }
1379 return GetLimit(*graph(sOpt + " readonly"));
1380}
1381
1383xRooNLLVar::xRooHypoSpace::findlimit(const char *opt, double relUncert, unsigned int maxTries)
1384{
1385 TString sOpt(opt);
1386 bool visualize = sOpt.Contains("visualize");
1387 sOpt.ReplaceAll("visualize", "");
1388 std::shared_ptr<TGraphErrors> gr = graph(sOpt + " readonly");
1389 if (visualize) {
1390 auto gra = graphs(sOpt.Contains("toys") ? "pcls readonly toys" : "pcls readonly");
1391 if (gra) {
1392 if (!gPad)
1393 gra->Draw(); // in 6.28 DrawClone wont make the gPad defined :( ... so Draw then clear and Draw Clone
1394 gPad->Clear();
1395 gra->DrawClone("A")->SetBit(kCanDelete);
1396 gPad->RedrawAxis();
1397 gra->GetHistogram()->SetMinimum(1e-9);
1398 gra->GetHistogram()->GetYaxis()->SetRangeUser(1e-9, 1);
1399 gPad->Modified();
1400#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1401 gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1402#endif
1404 }
1405 }
1406
1407 // resync parameter boundaries from nlls (may have been modified by fits)
1408 for (auto p : axes()) {
1409 for (auto &[pdf, nll] : fNlls) {
1410 if (auto _v = dynamic_cast<RooRealVar *>(nll->pars()->find(*p))) {
1411 dynamic_cast<RooRealVar *>(p)->setRange(_v->getMin(), _v->getMax());
1412 }
1413 }
1414 }
1415
1416 if (!gr || gr->GetN() < 2) {
1417 auto v = (axes().empty()) ? nullptr : dynamic_cast<RooRealVar *>(*axes().rbegin());
1418 if (!v)
1419 return std::pair(std::numeric_limits<double>::quiet_NaN(), 0.);
1420 double muMax = std::min(std::min(v->getMax("physical"), v->getMax()), v->getMax("scan"));
1421 double muMin = std::max(std::max(v->getMin("physical"), v->getMin()), v->getMin("scan"));
1422 if (!gr || gr->GetN() < 1) {
1423 if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), muMin)).getVal(sOpt).first)) {
1424 // first point failed ... give up
1425 Error("findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), muMin);
1426 return std::pair(std::numeric_limits<double>::quiet_NaN(), 0.);
1427 }
1428 gr.reset();
1429 return findlimit(opt, relUncert, maxTries - 1); // do this to resync parameter limits
1430 }
1431
1432 // can approximate expected limit using
1433 // mu_hat + sigma_mu*ROOT::Math::gaussian_quantile(1.-alpha/2.,1) for cls
1434 // or mu_hat + sigma_mu*ROOT::Math::gaussian_quantile((1.-alpha),1) for cls+b
1435 // get a very first estimate of sigma_mu from ufit to expected data, take error on mu as sigma_mu
1436 double nextPoint = muMin + (muMax - muMin) / 50;
1437
1438 // if done an expected limit, assume data is like expected and choose expected limit point as first test point
1439 if (sOpt.Contains("obs")) {
1440 TString sOpt2 = sOpt;
1441 sOpt2.ReplaceAll("obs", "exp");
1442 auto expLim = findlimit(sOpt2, std::numeric_limits<double>::infinity(), 0);
1443 if (!std::isnan(expLim.first) && expLim.first < nextPoint)
1444 nextPoint = expLim.first;
1445 }
1446
1447 auto point =
1448 (sOpt.Contains("exp")) ? back().asimov() : std::shared_ptr<xRooHypoPoint>(&back(), [](xRooHypoPoint *) {});
1449 point = nullptr;
1450 if (point && point->ufit()) {
1451 double rough_sigma_mu = point->mu_hat().getError();
1452 double another_estimate = point->mu_hat().getVal() + rough_sigma_mu * ROOT::Math::gaussian_quantile(0.95, 1);
1453 // if (another_estimate < nextPoint) {
1454 nextPoint = another_estimate;
1455 ::Info("xRooHypoSpace::findlimit", "Guessing %g based on rough sigma_mu = %g", nextPoint, rough_sigma_mu);
1456 //}
1457 }
1458
1459 if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), nextPoint)).getVal(sOpt).first)) {
1460 // second point failed ... give up
1461 Error("findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), nextPoint);
1462 return std::pair(std::numeric_limits<double>::quiet_NaN(), 0.);
1463 }
1464 gr.reset();
1465 return findlimit(opt, relUncert, maxTries - 1);
1466 }
1467
1468 auto lim = GetLimit(*gr);
1469
1470 if (std::isnan(lim.first)) {
1471 return lim;
1472 }
1473
1474 auto v = dynamic_cast<RooRealVar *>(*axes().rbegin());
1475 double maxMu = std::min(std::min(v->getMax("physical"), v->getMax()), v->getMax("scan"));
1476 double minMu = std::max(std::max(v->getMin("physical"), v->getMin()), v->getMin("scan"));
1477
1478 // static double MIN_LIMIT_UNCERT = 1e-4; // stop iterating once uncert gets this small
1479 if (lim.first > -std::numeric_limits<double>::infinity() && lim.first < std::numeric_limits<double>::infinity() &&
1480 (std::abs(lim.second) <= relUncert * std::abs(lim.first) /* || std::abs(lim.second)<MIN_LIMIT_UNCERT*/))
1481 return lim;
1482
1483 double nextPoint;
1484
1485 if (lim.second == std::numeric_limits<double>::infinity()) {
1486 // limit was found by extrapolating to right
1487 nextPoint = lim.first;
1488 if (nextPoint == std::numeric_limits<double>::infinity() || nextPoint > maxMu) {
1489 nextPoint = gr->GetPointX(gr->GetN() - 1) + (maxMu - minMu) / 50;
1490 }
1491
1492 // prefer extrapolation with sigma_mu, if available, if it takes us further
1493 // as shape of p-value curve is usually
1494 auto point =
1495 (sOpt.Contains("exp")) ? back().asimov() : std::shared_ptr<xRooHypoPoint>(&back(), [](xRooHypoPoint *) {});
1496 point = nullptr;
1497 if (point && point->ufit()) {
1498 double rough_sigma_mu = point->mu_hat().getError();
1499 double another_estimate = point->mu_hat().getVal() + rough_sigma_mu * ROOT::Math::gaussian_quantile(0.95, 1);
1500 // if (another_estimate < nextPoint) {
1501 nextPoint = std::max(nextPoint, another_estimate);
1502 ::Info("xRooHypoSpace::findlimit", "Guessing %g based on rough sigma_mu = %g", nextPoint, rough_sigma_mu);
1503 //}
1504 }
1505 nextPoint = std::min(nextPoint + nextPoint * relUncert * 0.99, maxMu); // ensure we step over location if possible
1506
1507 if (nextPoint > maxMu)
1508 return lim;
1509 } else if (lim.second == -std::numeric_limits<double>::infinity()) {
1510 // limit from extrapolating to left
1511 nextPoint = lim.first;
1512 if (nextPoint < minMu)
1513 nextPoint = gr->GetPointX(0) - (maxMu - minMu) / 50;
1514 if (nextPoint < minMu)
1515 return lim;
1516 } else {
1517 nextPoint = lim.first + lim.second * relUncert * 0.99;
1518 }
1519
1520 // got here need a new point .... evaluate the estimated lim location +/- the relUncert (signed error takes care of
1521 // direction)
1522
1523 ::Info("xRooHypoSpace::findlimit", "%s -- Testing new point @ %s=%g (delta=%g)", sOpt.Data(), v->GetName(),
1524 nextPoint, lim.second);
1525 if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), nextPoint)).getVal(sOpt).first)) {
1526 if (maxTries == 0) {
1527 Warning("findlimit", "Reached max number of point evaluations");
1528 } else {
1529 Error("findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), nextPoint);
1530 }
1531 return lim;
1532 }
1533 gr.reset();
1534 return findlimit(opt, relUncert, maxTries - 1);
1535}
1536
1538{
1539
1540 TString sOpt(opt);
1541 sOpt.ToLower();
1542
1543 if ((sOpt == "" || sOpt == "same") && !empty()) {
1544 if (front().fPllType == xRooFit::Asymptotics::OneSidedPositive) {
1545 sOpt += "pcls"; // default to showing cls p-value scan if drawing a limit
1546 for (auto &hp : *this) {
1547 if (!hp.nullToys.empty() || !hp.altToys.empty()) {
1548 sOpt += " toys";
1549 break; // default to toys if done toys
1550 }
1551 }
1552 } else if (front().fPllType == xRooFit::Asymptotics::TwoSided) {
1553 sOpt += "ts";
1554 }
1555 }
1556
1557 // split up by ; and call Draw for each (with 'same' appended)
1558 auto _axes = axes();
1559 if (_axes.empty())
1560 return;
1561
1562 if (sOpt == "status") {
1563 // draw the points in the space
1564 if (_axes.size() <= 2) {
1565 TGraphErrors *out = new TGraphErrors;
1566 out->SetBit(kCanDelete);
1567 out->SetName("points");
1568 out->SetMarkerSize(0.5);
1569 TGraph *tsAvail = new TGraph;
1570 tsAvail->SetName("ts");
1571 tsAvail->SetBit(kCanDelete);
1572 tsAvail->SetMarkerStyle(20);
1573 TGraph *expAvail = new TGraph;
1574 expAvail->SetName("exp");
1575 expAvail->SetBit(kCanDelete);
1576 expAvail->SetMarkerStyle(25);
1577 expAvail->SetMarkerSize(out->GetMarkerSize() * 1.5);
1578 TGraph *badPoints = new TGraph;
1579 badPoints->SetName("bad_ufit");
1580 badPoints->SetBit(kCanDelete);
1581 badPoints->SetMarkerStyle(5);
1582 badPoints->SetMarkerColor(kRed);
1583 badPoints->SetMarkerSize(out->GetMarkerSize());
1584 TGraph *badPoints2 = new TGraph;
1585 badPoints2->SetName("bad_cfit_null");
1586 badPoints2->SetBit(kCanDelete);
1587 badPoints2->SetMarkerStyle(2);
1588 badPoints2->SetMarkerColor(kRed);
1589 badPoints2->SetMarkerSize(out->GetMarkerSize());
1590
1591 out->SetTitle(TString::Format("%s;%s;%s", GetTitle(), _axes.at(0)->GetTitle(),
1592 (_axes.size() == 1) ? "" : _axes.at(1)->GetTitle()));
1593 for (auto &p : *this) {
1594 bool _readOnly = p.nllVar ? p.nllVar->get()->getAttribute("readOnly") : false;
1595 if (p.nllVar)
1596 p.nllVar->get()->setAttribute("readOnly", true);
1597 double x = p.coords->getRealValue(_axes.at(0)->GetName());
1598 double y = _axes.size() == 1 ? p.ts_asymp().first : p.coords->getRealValue(_axes.at(1)->GetName());
1599 out->SetPoint(out->GetN(), x, y);
1600 if (!std::isnan(p.ts_asymp().first)) {
1601 if (_axes.size() == 1)
1602 out->SetPointError(out->GetN() - 1, 0, p.ts_asymp().second);
1603 tsAvail->SetPoint(tsAvail->GetN(), x, y);
1604 } else if (p.fUfit && (std::isnan(p.fUfit->minNll()) ||
1605 xRooNLLVar::xRooHypoPoint::allowedStatusCodes.find(p.fUfit->status()) ==
1607 badPoints->SetPoint(badPoints->GetN(), x, y);
1608 } else if (p.fNull_cfit && (std::isnan(p.fNull_cfit->minNll()) ||
1609 xRooNLLVar::xRooHypoPoint::allowedStatusCodes.find(p.fNull_cfit->status()) ==
1611 badPoints2->SetPoint(badPoints2->GetN(), x, y);
1612 }
1613 if (!std::isnan(p.ts_asymp(0).first)) {
1614 expAvail->SetPoint(expAvail->GetN(), x, y);
1615 } else if (p.asimov() && p.asimov()->fUfit &&
1616 (std::isnan(p.asimov()->fUfit->minNll()) ||
1617 xRooNLLVar::xRooHypoPoint::allowedStatusCodes.find(p.asimov()->fUfit->status()) ==
1619
1620 } else if (p.asimov() && p.asimov()->fNull_cfit &&
1621 (std::isnan(p.asimov()->fNull_cfit->minNll()) ||
1622 xRooNLLVar::xRooHypoPoint::allowedStatusCodes.find(p.asimov()->fNull_cfit->status()) ==
1624 }
1625 if (p.nllVar)
1626 p.nllVar->get()->setAttribute("readOnly", _readOnly);
1627 }
1628
1629 if (_axes.size() == 1) {
1630 TGraph tmp;
1631 for (int i = 0; i < out->GetN(); i++) {
1632 if (!std::isnan(out->GetPointY(i)))
1633 tmp.SetPoint(tmp.GetN(), out->GetPointX(i), out->GetPointY(i));
1634 }
1635 auto fixPoints = [&](TGraph *g) {
1636 for (int i = 0; i < g->GetN(); i++) {
1637 if (std::isnan(g->GetPointY(i)))
1638 g->SetPointY(i, std::isnan(tmp.Eval(g->GetPointX(i))) ? 0. : tmp.Eval(g->GetPointX(i)));
1639 }
1640 };
1641 fixPoints(out);
1642 fixPoints(tsAvail);
1643 fixPoints(expAvail);
1644 fixPoints(badPoints);
1645 fixPoints(badPoints2);
1646 }
1647
1648 out->SetMarkerStyle(4);
1649 out->Draw("AP");
1650 auto leg = new TLegend(1. - gPad->GetRightMargin() - 0.3, 1. - gPad->GetTopMargin() - 0.35,
1651 1. - gPad->GetRightMargin() - 0.05, 1. - gPad->GetTopMargin() - 0.05);
1652 leg->SetName("legend");
1653 leg->AddEntry(out, "Uncomputed", "P");
1654
1655 if (tsAvail->GetN()) {
1656 out->GetListOfFunctions()->Add(tsAvail, "P");
1657 leg->AddEntry(tsAvail, "Computed", "P");
1658 } else {
1659 delete tsAvail;
1660 }
1661 if (expAvail->GetN()) {
1662 out->GetListOfFunctions()->Add(expAvail, "P");
1663 leg->AddEntry(expAvail, "Expected computed", "P");
1664 } else {
1665 delete expAvail;
1666 }
1667 if (badPoints->GetN()) {
1668 out->GetListOfFunctions()->Add(badPoints, "P");
1669 leg->AddEntry(badPoints, "Bad ufit", "P");
1670 } else {
1671 delete badPoints;
1672 }
1673 if (badPoints2->GetN()) {
1674 out->GetListOfFunctions()->Add(badPoints2, "P");
1675 leg->AddEntry(badPoints2, "Bad null cfit", "P");
1676 } else {
1677 delete badPoints2;
1678 }
1680 leg->Draw();
1681 // if(_axes.size()==1) out->GetHistogram()->GetYaxis()->SetRangeUser(0,1);
1682 gPad->SetGrid(true, _axes.size() > 1);
1683 if (_axes.size() == 1)
1684 gPad->SetLogy(false);
1685 }
1686
1688
1689 return;
1690 }
1691
1692 if (sOpt.Contains("pcls") || sOpt.Contains("pnull") || sOpt.Contains("ts")) {
1693 auto gra = graphs(sOpt + " readonly");
1694 if (!gPad && gROOT->GetSelectedPad())
1695 gROOT->GetSelectedPad()->cd();
1696 if (!sOpt.Contains("same") && gPad) {
1697 gPad->Clear();
1698 }
1699 if (gra) {
1700 auto gra2 = static_cast<TMultiGraph *>(gra->DrawClone(sOpt.Contains("same") ? "" : "A"));
1701 gra2->SetBit(kCanDelete);
1702 if (sOpt.Contains("pcls") || sOpt.Contains("pnull")) {
1703 gra2->GetHistogram()->SetMinimum(1e-6);
1704 }
1705 if (gPad) {
1706 gPad->RedrawAxis();
1707 gPad->GetCanvas()->Paint();
1708 gPad->GetCanvas()->Update();
1709#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1710 gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1711#endif
1713 }
1714 }
1715 if (!sOpt.Contains("same") && gPad) {
1716 // auto mg = static_cast<TMultiGraph*>(gPad->GetPrimitive(gra->GetName()));
1717 // mg->GetHistogram()->SetMinimum(1e-9);
1718 // mg->GetHistogram()->GetYaxis()->SetRangeUser(1e-9,1);
1719 // gPad->SetGrid(0, 0);
1720 // gPad->SetLogy(1);
1721 }
1722
1724
1725 return;
1726 }
1727
1728 // graphs("ts visualize");return;
1729
1730 TGraphErrors *out = new TGraphErrors;
1731 out->SetName(GetName());
1732
1733 TString title = (!axes().empty()) ? TString::Format(";%s", axes().first()->GetTitle()) : "";
1734
1735 auto pllType = xRooFit::Asymptotics::TwoSided;
1736 if (!empty() && axes().size() == 1) {
1737 for (auto &p : *this) {
1738 if (p.fPllType != xRooFit::Asymptotics::TwoSided) {
1739 pllType = p.fPllType;
1740 }
1741 }
1742 title += ";";
1743 title += front().tsTitle(true);
1744 }
1745
1746 out->SetTitle(title);
1747 *dynamic_cast<TAttFill *>(out) = *this;
1748 *dynamic_cast<TAttLine *>(out) = *this;
1749 *dynamic_cast<TAttMarker *>(out) = *this;
1750 out->SetBit(kCanDelete);
1751
1752 if (!gPad)
1754 TVirtualPad *basePad = gPad;
1755 if (!sOpt.Contains("same"))
1756 basePad->Clear();
1757
1758 bool doFits = false;
1759 if (sOpt.Contains("fits")) {
1760 doFits = true;
1761 sOpt.ReplaceAll("fits", "");
1762 }
1763
1764 auto mainPad = gPad;
1765
1766 out->SetEditable(false);
1767
1768 if (doFits) {
1769 gPad->Divide(1, 2);
1770 gPad->cd(1);
1771 gPad->SetBottomMargin(gPad->GetBottomMargin() * 2.); // increase margin to be same as before
1772 gPad->SetGrid();
1773 out->Draw(sOpt);
1774 basePad->cd(2);
1775 mainPad = basePad->GetPad(1);
1776 } else {
1777 gPad->SetGrid();
1778 out->Draw("ALP");
1779 }
1780
1781 std::pair<double, double> minMax(std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity());
1782 for (auto &p : *this) {
1783 if (p.fPllType != pllType)
1784 continue; // must all have same pll type
1785 auto val = p.pll(true).first;
1786 if (std::isnan(val))
1787 continue;
1788 minMax.first = std::min(minMax.first, val);
1789 minMax.second = std::max(minMax.second, val);
1790 }
1791 if (minMax.first < std::numeric_limits<double>::infinity())
1792 out->GetHistogram()->SetMinimum(minMax.first);
1793 if (minMax.second > -std::numeric_limits<double>::infinity())
1794 out->GetHistogram()->SetMaximum(minMax.second);
1795
1796 TGraph *badPoints = nullptr;
1797
1798 TStopwatch s;
1799 s.Start();
1800 std::shared_ptr<const RooFitResult> ufr;
1801 for (auto &p : *this) {
1802 if (p.fPllType != pllType)
1803 continue; // must all have same pll type
1804 auto val = p.pll().first;
1805 if (!ufr)
1806 ufr = p.ufit();
1807 if (out->GetN() == 0 && ufr && ufr->status() == 0) {
1808 out->SetPoint(out->GetN(),
1809 ufr->floatParsFinal().getRealValue(axes().first()->GetName(),
1810 ufr->constPars().getRealValue(axes().first()->GetName())),
1811 0.);
1812 out->SetPointError(out->GetN() - 1, 0, ufr->edm());
1813 }
1814 if (auto fr = p.fNull_cfit;
1815 fr && doFits) { // access member to avoid unnecessarily creating fit result if wasnt needed
1816 // create a new subpad and draw fitResult on it
1817 auto _pad = gPad;
1818 auto pad =
1819 new TPad(fr->GetName(), TString::Format("%s = %g", poi().first()->GetTitle(), p.fNullVal()), 0, 0, 1., 1);
1820 pad->SetNumber(out->GetN() + 1); // can't use "0" for a subpad
1821 pad->cd();
1822 xRooNode(fr).Draw("goff");
1823 _pad->cd();
1824 //_pad->GetListOfPrimitives()->AddFirst(pad);
1825 pad->AppendPad();
1826 }
1827 if (std::isnan(val) && p.status() != 0) {
1828 if (!badPoints) {
1829 badPoints = new TGraph;
1830 badPoints->SetBit(kCanDelete);
1831 badPoints->SetName("badPoints");
1832 badPoints->SetMarkerStyle(5);
1833 badPoints->SetMarkerColor(kRed);
1834 badPoints->SetMarkerSize(1);
1835 out->GetListOfFunctions()->Add(badPoints, "P");
1836 }
1837 badPoints->SetPoint(badPoints->GetN(), p.fNullVal(), out->Eval(p.fNullVal()));
1838 mainPad->Modified();
1839 } else if (!std::isnan(val)) {
1840 out->SetPoint(out->GetN(), p.coords->getRealValue(axes().first()->GetName()), p.pll().first);
1841 out->SetPointError(out->GetN() - 1, 0, p.pll().second);
1842 out->Sort();
1843
1844 // reposition bad points
1845 if (badPoints) {
1846 for (int i = 0; i < badPoints->GetN(); i++)
1847 badPoints->SetPointY(i, out->Eval(badPoints->GetPointX(i)));
1848 }
1849
1850 mainPad->Modified();
1851 }
1852 if (s.RealTime() > 3) { // stops the clock
1853 basePad->GetCanvas()->Paint();
1854 basePad->GetCanvas()->Update();
1856 s.Reset();
1857 s.Start();
1858 }
1859 s.Continue();
1860 }
1861 basePad->GetCanvas()->Paint();
1862 basePad->GetCanvas()->Update();
1863#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1864 basePad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1865#endif
1867
1868 // finish by overlaying ufit
1869 if (ufr && doFits) {
1870 auto _pad = gPad;
1871 auto pad = new TPad(ufr->GetName(), "unconditional fit", 0, 0, 1., 1.);
1872 pad->SetNumber(-1);
1873 pad->cd();
1874 xRooNode(ufr).Draw("goff");
1875 _pad->cd();
1876 pad->AppendPad();
1877
1878 // draw one more pad to represent the selected, and draw the ufit pad onto that pad
1879 pad = new TPad("selected", "selected", 0, 0, 1, 1);
1880 pad->Draw();
1881 pad->cd();
1882 basePad->GetPad(2)->GetPad(-1)->AppendPad();
1883 pad->Modified();
1884 pad->Update();
1886
1887 basePad->cd();
1888 }
1889
1890 if (doFits) {
1891 if (!xRooNode::gIntObj) {
1893 }
1894 gPad->GetCanvas()->Connect("Highlighted(TVirtualPad*,TObject*,Int_t,Int_t)", "xRooNode::InteractiveObject",
1895 xRooNode::gIntObj, "Interactive_PLLPlot(TVirtualPad*,TObject*,Int_t,Int_t)");
1896 }
1897
1898 return;
1899}
1900
1902{
1903
1904 RooStats::HypoTestInverterResult *out = nullptr;
1905
1906 auto _axes = axes();
1907 if (_axes.empty())
1908 return out;
1909
1910 out = new RooStats::HypoTestInverterResult(GetName(), *dynamic_cast<RooRealVar *>(_axes.at(0)), 0.95);
1911 out->SetTitle(GetTitle());
1912
1913 for (auto &p : *this) {
1914 double _x = p.coords->getRealValue(_axes.at(0)->GetName(), std::numeric_limits<double>::quiet_NaN());
1915 out->Add(_x, p.result());
1916 }
1917
1918 return out;
1919}
1920
#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
double _val
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:367
#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:1983
std::shared_ptr< const RooAbsCollection > coords
Definition xRooNLLVar.h:242
xValueWithError limit(const char *type="cls", double nSigma=std::numeric_limits< double >::quiet_NaN()) const
static xValueWithError GetLimit(const TGraph &pValues, double target=std::numeric_limits< double >::quiet_NaN())
std::map< std::string, xValueWithError > 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())
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)
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::shared_ptr< RooAbsPdf > pdf() const
Definition xRooNLLVar.h:423
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:52
void Draw(Option_t *opt="") override
Default Draw method for all objects.
static InteractiveObject * gIntObj
Definition xRooNode.h:470
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.
A space to attach TBranches.
virtual value_type getCurrentIndex() const
Return index number of current state.
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:2476
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:4089
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:1527
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2325
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:688
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:2469
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:2020
void SetName(const char *name="") override
Set graph name.
Definition TGraph.cxx:2364
virtual Double_t GetPointY(Int_t i) const
Get y value for point i.
Definition TGraph.cxx:1538
virtual void SetPointY(Int_t i, Double_t y)
Set y value for point i.
Definition TGraph.cxx:2357
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:2284
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.
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
const char * Data() const
Definition TString.h:376
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition TString.h:704
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 BEGIN_XROOFIT_NAMESPACE
Definition Config.h:24
#define END_XROOFIT_NAMESPACE
Definition Config.h:25
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4