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