Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
xRooNLLVar.cxx
Go to the documentation of this file.
1/*
2 * Project: xRooFit
3 * Author:
4 * Will Buttinger, RAL 2022
5 *
6 * Copyright (c) 2022, CERN
7 *
8 * Redistribution and use in source and binary forms,
9 * with or without modification, are permitted according to the terms
10 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
11 */
12
13/** \class ROOT::Experimental::XRooFit::xRooNLLVar
14\ingroup xroofit
15
16This xRooNLLVar object has several special methods, e.g. for fitting and toy dataset generation.
17
18 */
19
20#include "RVersion.h"
21
22#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
23#define protected public
24#endif
25#include "RooFitResult.h"
26#include "RooNLLVar.h"
27#ifdef protected
28#undef protected
29#endif
30
31#include "xRooFit/xRooFit.h"
32
33#include "RooCmdArg.h"
34#include "RooAbsPdf.h"
35#include "RooAbsData.h"
36
37#include "RooConstraintSum.h"
38#include "RooSimultaneous.h"
40#include "TPRegexp.h"
41#include "TEfficiency.h"
42
43#include "RooRealVar.h"
44#include "Math/ProbFunc.h"
45#include "RooRandom.h"
46
47#include "TPad.h"
48#include "TSystem.h"
49
50#include "coutCapture.h"
51
52#include <chrono>
53
54#include "Math/GenAlgoOptions.h"
55
56#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
57#define private public
58#define GETWS(a) a->_myws
59#define GETWSSETS(w) w->_namedSets
60#else
61#define GETWS(a) a->workspace()
62#define GETWSSETS(w) w->sets()
63#endif
64#include "RooWorkspace.h"
65#ifdef private
66#undef private
67#endif
68
69#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
70#define protected public
71#endif
73#ifdef protected
74#undef protected
75#endif
76
77#include "TMultiGraph.h"
78#include "TCanvas.h"
79#include "TArrow.h"
80#include "RooStringVar.h"
81#include "TDirectory.h"
82#include "TStyle.h"
83#include "TH1D.h"
84#include "TLegend.h"
85#include "RooCategory.h"
86#include "TTree.h"
87#include "TGraph2D.h"
88
89#include "RooGaussian.h"
90#include "RooPoisson.h"
91
92#include "TROOT.h"
93#include "TKey.h"
94#include "TRegexp.h"
95
97
98std::set<int> xRooNLLVar::xRooHypoPoint::allowedStatusCodes = {0};
99
101public:
102 AutoRestorer(const RooAbsCollection &s, xRooNLLVar *nll = nullptr) : fSnap(s.snapshot()), fNll(nll)
103 {
104 fPars.add(s);
105 if (fNll) {
106 // if (!fNll->kReuseNLL) fOldNll = *fNll;
107 fOldData = fNll->getData();
108 fOldName = fNll->get()->GetName();
109 fOldTitle = fNll->get()->getStringAttribute("fitresultTitle");
110 }
111 }
113 {
115 if (fNll) {
116 // commented out code was attempt to speed up things avoid unnecessarily reinitializing things over and over
117 // if (!fNll->kReuseNLL) {
118 // // can be faster just by putting back in old nll
119 // fNll->std::shared_ptr<RooAbsReal>::operator=(fOldNll);
120 // fNll->fData = fOldData.first;
121 // fNll->fGlobs = fOldData.second;
122 // } else {
123 // fNll->setData(fOldData);
124 // fNll->get()->SetName(fOldName);
125 // fNll->get()->setStringAttribute("fitresultTitle", (fOldTitle == "") ? nullptr : fOldTitle);
126 // }
127 fNll->fGlobs = fOldData.second; // will mean globs matching checks are skipped in setData
128 fNll->setData(fOldData);
129 fNll->get()->SetName(fOldName);
130 fNll->get()->setStringAttribute("fitresultTitle", (fOldTitle == "") ? nullptr : fOldTitle);
131 }
132 }
134 std::unique_ptr<RooAbsCollection> fSnap;
135 xRooNLLVar *fNll = nullptr;
136 // std::shared_ptr<RooAbsReal> fOldNll;
137 std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> fOldData;
139};
140
141xRooNLLVar::~xRooNLLVar() {}
142
143xRooNLLVar::xRooNLLVar(RooAbsPdf &pdf, const std::pair<RooAbsData *, const RooAbsCollection *> &data,
144 const RooLinkedList &nllOpts)
145 : xRooNLLVar(std::shared_ptr<RooAbsPdf>(&pdf, [](RooAbsPdf *) {}),
146 std::make_pair(std::shared_ptr<RooAbsData>(data.first, [](RooAbsData *) {}),
147 std::shared_ptr<const RooAbsCollection>(data.second, [](const RooAbsCollection *) {})),
148 nllOpts)
149{
150}
151
152xRooNLLVar::xRooNLLVar(const std::shared_ptr<RooAbsPdf> &pdf,
153 const std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> &data,
154 const RooLinkedList &opts)
155 : fPdf(pdf), fData(data.first), fGlobs(data.second)
156{
157
159
160 fOpts = std::shared_ptr<RooLinkedList>(new RooLinkedList, [](RooLinkedList *l) {
161 if (l)
162 l->Delete();
163 delete l;
164 });
165 fOpts->SetName("");
166
167 // we *must* take global observables from the model even if they are included in the dataset
168 // this is because the way xRooNLLVar is coded up it assumes the globs in the funcVars *ARE*
169 // part of the model
170 fOpts->Add(RooFit::GlobalObservablesSource("model").Clone(nullptr));
171
172 for (int i = 0; i < opts.GetSize(); i++) {
173 if (strlen(opts.At(i)->GetName()) == 0)
174 continue; // skipping "none" cmds
175 if (strcmp(opts.At(i)->GetName(), "GlobalObservables") == 0) {
176 // will skip here to add with the obs from the function below
177 // must match global observables
178 auto gl = dynamic_cast<RooCmdArg *>(opts.At(i))->getSet(0);
179 if (!fGlobs || !fGlobs->equals(*gl)) {
180 throw std::runtime_error("GlobalObservables mismatch");
181 }
182 } else if (strcmp(opts.At(i)->GetName(), "Hesse") == 0) {
183 fitConfig()->SetParabErrors(dynamic_cast<RooCmdArg *>(opts.At(i))->getInt(0)); // controls hesse
184 } else if (strcmp(opts.At(i)->GetName(), "Strategy") == 0) {
185 fitConfig()->MinimizerOptions().SetStrategy(dynamic_cast<RooCmdArg *>(opts.At(i))->getInt(0));
186 } else if (strcmp(opts.At(i)->GetName(), "StrategySequence") == 0) {
187 fitConfigOptions()->SetNamedValue("StrategySequence", dynamic_cast<RooCmdArg *>(opts.At(i))->getString(0));
188 } else if (strcmp(opts.At(i)->GetName(), "Tolerance") == 0) {
189 fitConfig()->MinimizerOptions().SetTolerance(dynamic_cast<RooCmdArg *>(opts.At(i))->getInt(0));
190 } else if (strcmp(opts.At(i)->GetName(), "PrintLevel") == 0) {
191 fitConfig()->MinimizerOptions().SetPrintLevel(dynamic_cast<RooCmdArg *>(opts.At(i))->getInt(0));
192 } else {
193 if (strcmp(opts.At(i)->GetName(), "Optimize") == 0) {
194 // this flag will trigger constOptimizeTestStatistic to be called on the nll in createNLL method
195 // we should ensure that the fitconfig setting is consistent with it ...
196 fitConfigOptions()->SetValue("OptimizeConst", dynamic_cast<RooCmdArg *>(opts.At(i))->getInt(0));
197 }
198 fOpts->Add(opts.At(i)->Clone(nullptr)); // nullptr needed because accessing Clone via TObject base class puts
199 // "" instead, so doesnt copy names
200 }
201 }
202 if (fGlobs) {
203 // add global observables opt with function obs
204 auto _vars = std::unique_ptr<RooArgSet>(fPdf->getVariables());
205 if (auto extCon = dynamic_cast<RooCmdArg *>(fOpts->find("ExternalConstraints"))) {
206 for (auto con : *extCon->getSet(0)) {
207 _vars->add(*std::unique_ptr<RooArgSet>(con->getVariables()));
208 }
209 }
210 auto _funcGlobs = std::unique_ptr<RooArgSet>(dynamic_cast<RooArgSet *>(_vars->selectCommon(*fGlobs)));
211 fOpts->Add(RooFit::GlobalObservables(*_funcGlobs).Clone());
212 }
213
214 if (auto flag = dynamic_cast<RooCmdArg *>(fOpts->find("ReuseNLL"))) {
215 kReuseNLL = flag->getInt(0);
216 }
217
218 // if fit range specified, and pdf is a RooSimultaneous, may need to 'reduce' the model if some of the pdfs are in
219 // range and others are not
220 if (auto range = dynamic_cast<RooCmdArg *>(fOpts->find("RangeWithName"))) {
221 TString rangeName = range->getString(0);
222
223 // reduce the data here for convenience, not really necessary because will happen inside RooNLLVar but still
224 // fData.reset( fData->reduce(RooFit::SelectVars(*fData->get()),RooFit::CutRange(rangeName)) );
225
226 if (auto s = dynamic_cast<RooSimultaneous *>(fPdf.get()); s) {
227 auto &_cat = const_cast<RooAbsCategoryLValue &>(s->indexCat());
228 std::vector<TString> chanPatterns;
229 TStringToken pattern(rangeName, ",");
230 bool hasRange(false);
231 std::string noneCatRanges;
232 while (pattern.NextToken()) {
233 chanPatterns.emplace_back(pattern);
234 if (_cat.hasRange(chanPatterns.back())) {
235 hasRange = true;
236 } else {
237 if (!noneCatRanges.empty())
238 noneCatRanges += ",";
239 noneCatRanges += chanPatterns.back();
240 }
241 }
242 if (hasRange) {
243 // must remove the ranges that referred to selections on channel category
244 // otherwise RooFit will incorrectly evaluate the NLL (it creates a partition for each range given in the
245 // list, which all end up being equal) the NLL would become scaled by the number of ranges given
246 if (noneCatRanges.empty()) {
247 fOpts->Remove(range);
248 SafeDelete(range);
249 } else {
250 range->setString(0, noneCatRanges.c_str());
251 }
252 // must reduce because category var has one of the ranges
253 auto newPdf =
254 std::make_shared<RooSimultaneous>(TString::Format("%s_reduced", s->GetName()), "Reduced model", _cat);
255 for (auto &c : _cat) {
256 auto _pdf = s->getPdf(c.first.c_str());
257 if (!_pdf)
258 continue;
259 _cat.setIndex(c.second);
260 bool matchAny = false;
261 for (auto &p : chanPatterns) {
262 if (_cat.hasRange(p) && _cat.inRange(p)) {
263 matchAny = true;
264 break;
265 }
266 }
267 if (matchAny) {
268 newPdf->addPdf(*_pdf, c.first.c_str());
269 }
270 }
271 fPdf = newPdf;
272 }
273 }
274 }
275
276 // if (fGlobs) {
277 // // must check GlobalObservables is in the list
278 // }
279 //
280 // if (auto globs = dynamic_cast<RooCmdArg*>(fOpts->find("GlobalObservables"))) {
281 // // first remove any obs the pdf doesnt depend on
282 // auto _vars = std::unique_ptr<RooAbsCollection>( fPdf->getVariables() );
283 // auto _funcGlobs = std::unique_ptr<RooAbsCollection>(_vars->selectCommon(*globs->getSet(0)));
284 // fGlobs.reset( std::unique_ptr<RooAbsCollection>(globs->getSet(0)->selectCommon(*_funcGlobs))->snapshot() );
285 // globs->setSet(0,dynamic_cast<const RooArgSet&>(*_funcGlobs)); // globs in linked list has its own argset
286 // but args need to live as long as the func
287 // /*RooArgSet toRemove;
288 // for(auto a : *globs->getSet(0)) {
289 // if (!_vars->find(*a)) toRemove.add(*a);
290 // }
291 // const_cast<RooArgSet*>(globs->getSet(0))->remove(toRemove);
292 // fGlobs.reset( globs->getSet(0)->snapshot() );
293 // fGlobs->setAttribAll("Constant",true);
294 // const_cast<RooArgSet*>(globs->getSet(0))->replace(*fGlobs);*/
295 // }
296}
297
298xRooNLLVar::xRooNLLVar(const std::shared_ptr<RooAbsPdf> &pdf, const std::shared_ptr<RooAbsData> &data,
299 const RooLinkedList &opts)
300 : xRooNLLVar(
301 pdf,
302 std::make_pair(data, std::shared_ptr<const RooAbsCollection>(
303 (opts.find("GlobalObservables"))
304 ? dynamic_cast<RooCmdArg *>(opts.find("GlobalObservables"))->getSet(0)->snapshot()
305 : nullptr)),
306 opts)
307{
308}
309
311{
312 std::cout << "PDF: ";
313 if (fPdf) {
314 fPdf->Print();
315 } else {
316 std::cout << "<null>" << std::endl;
317 }
318 std::cout << "Data: ";
319 if (fData) {
320 fData->Print();
321 } else {
322 std::cout << "<null>" << std::endl;
323 }
324 std::cout << "NLL Options: " << std::endl;
325 for (int i = 0; i < fOpts->GetSize(); i++) {
326 auto c = dynamic_cast<RooCmdArg *>(fOpts->At(i));
327 if (!c)
328 continue;
329 std::cout << " " << c->GetName() << " : ";
330 if (c->getString(0)) {
331 std::cout << c->getString(0);
332 } else if (c->getSet(0) && !c->getSet(0)->empty()) {
333 std::cout << (c->getSet(0)->contentsString());
334 } else {
335 std::cout << c->getInt(0);
336 }
337 std::cout << std::endl;
338 }
339 if (fFitConfig) {
340 std::cout << "Fit Config: " << std::endl;
341 std::cout << " UseParabErrors: " << (fFitConfig->ParabErrors() ? "True" : "False")
342 << " [toggles HESSE algorithm]" << std::endl;
343 std::cout << " MinimizerOptions: " << std::endl;
344 fFitConfig->MinimizerOptions().Print();
345 }
346 std::cout << "Last Rebuild Log Output: " << fFuncCreationLog << std::endl;
347}
348
350{
351 TString oldName = "";
352 if (std::shared_ptr<RooAbsReal>::get())
353 oldName = std::shared_ptr<RooAbsReal>::get()->GetName();
354 if (fPdf) {
356 // need to find all RooRealSumPdf nodes and mark them binned or unbinned as required
357 RooArgSet s;
358 fPdf->treeNodeServerList(&s, nullptr, true, false);
359 s.add(*fPdf); // ensure include self in case fitting a RooRealSumPdf
360 bool isBinned = false;
361 bool hasBinned = false; // if no binned option then 'auto bin' ...
362 if (auto a = dynamic_cast<RooCmdArg *>(fOpts->find("Binned")); a) {
363 hasBinned = true;
364 isBinned = a->getInt(0);
365 }
366 std::map<RooAbsArg *, bool> origValues;
367 if (hasBinned) {
368 for (auto a : s) {
369 if (a->InheritsFrom("RooRealSumPdf")) {
370 // since RooNLLVar will assume binBoundaries available (not null), we should check bin boundaries
371 // available
372 bool setBinned = false;
373 if (isBinned) {
374 std::unique_ptr<RooArgSet> obs(a->getObservables(fData->get()));
375 if (obs->size() == 1) { // RooNLLVar requires exactly 1 obs
376 auto *var = static_cast<RooRealVar *>(obs->first());
377 std::unique_ptr<std::list<double>> boundaries{dynamic_cast<RooAbsReal *>(a)->binBoundaries(
378 *var, -std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity())};
379 if (boundaries) {
380 if (!std::shared_ptr<RooAbsReal>::get()) {
381 Info("xRooNLLVar", "%s will be evaluated as a Binned PDF (%d bins)", a->GetName(),
382 int(boundaries->size() - 1));
383 }
384 setBinned = true;
385 }
386 }
387 }
388 origValues[a] = a->getAttribute("BinnedLikelihood");
389 a->setAttribute("BinnedLikelihood", setBinned);
390 }
391 }
392 }
393 // before creating, clear away caches if any if pdf is in ws
394 if (GETWS(fPdf)) {
395 std::set<std::string> setNames;
396 for (auto &a : GETWSSETS(GETWS(fPdf))) {
397 if (TString(a.first.c_str()).BeginsWith("CACHE_")) {
398 setNames.insert(a.first);
399 }
400 }
401 for (auto &a : setNames) {
402 GETWS(fPdf)->removeSet(a.c_str());
403 }
404 }
405 std::set<std::string> attribs;
406 if (std::shared_ptr<RooAbsReal>::get())
407 attribs = std::shared_ptr<RooAbsReal>::get()->attributes();
408 this->reset(std::unique_ptr<RooAbsReal>{fPdf->createNLL(*fData, *fOpts)}.release());
409 std::shared_ptr<RooAbsReal>::get()->SetName(TString::Format("nll_%s/%s", fPdf->GetName(), fData->GetName()));
410 // RooFit only swaps in what it calls parameters, this misses out the RooConstVars which we treat as pars as well
411 // so swap those in ... question: is recursiveRedirectServers usage in RooAbsOptTestStatic (and here) a memory
412 // leak?? where do the replaced servers get deleted??
413
414 for (auto &a : attribs)
415 std::shared_ptr<RooAbsReal>::get()->setAttribute(a.c_str());
416 // create parent on next line to avoid triggering workspace initialization code in constructor of xRooNode
417 if (GETWS(fPdf)) {
418 xRooNode(*GETWS(fPdf), std::make_shared<xRooNode>()).sterilize();
419 } // there seems to be a nasty bug somewhere that can make the cache become invalid, so clear it here
420 if (oldName != "")
421 std::shared_ptr<RooAbsReal>::get()->SetName(oldName);
422 if (!origValues.empty()) {
423 // need to evaluate NOW so that slaves are created while the BinnedLikelihood settings are in place
424 std::shared_ptr<RooAbsReal>::get()->getVal();
425 for (auto &[o, v] : origValues)
426 o->setAttribute("BinnedLikelihood", v);
427 }
428 }
429
430 fFuncVars = std::unique_ptr<RooArgSet>{std::shared_ptr<RooAbsReal>::get()->getVariables()};
431 if (fGlobs) {
432 fFuncGlobs.reset(fFuncVars->selectCommon(*fGlobs));
433 fFuncGlobs->setAttribAll("Constant", true);
434 }
435 fConstVars.reset(fFuncVars->selectByAttrib("Constant", true)); // will check if any of these have floated
436}
437
438std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>>
439xRooNLLVar::generate(bool expected, int seed)
440{
441 if (!fPdf)
442 return std::pair(nullptr, nullptr);
443 auto fr = std::make_shared<RooFitResult>(TUUID().AsString());
444 fr->setFinalParList(RooArgList());
446 l.add((fFuncVars) ? *fFuncVars : *std::unique_ptr<RooAbsCollection>(fPdf->getParameters(*fData)));
447 fr->setConstParList(l);
448 const_cast<RooArgList &>(fr->constPars()).setAttribAll("global", false);
449 if (fGlobs)
450 std::unique_ptr<RooAbsCollection>(fr->constPars().selectCommon(*fGlobs))->setAttribAll("global", true);
451 return xRooFit::generateFrom(*fPdf, *fr, expected, seed);
452}
453
454xRooNLLVar::xRooFitResult::xRooFitResult(const std::shared_ptr<xRooNode> &in, const std::shared_ptr<xRooNLLVar> &nll)
455 : std::shared_ptr<const RooFitResult>(std::dynamic_pointer_cast<const RooFitResult>(in->fComp)),
456 fNode(in),
457 fNll(nll),
458 fCfits(std::make_shared<std::map<std::string, xRooFitResult>>())
459{
460}
462{
463 return fNode->get<RooFitResult>();
464}
465// xRooNLLVar::xRooFitResult::operator std::shared_ptr<const RooFitResult>() const { return
466// std::dynamic_pointer_cast<const RooFitResult>(fNode->fComp); }
467xRooNLLVar::xRooFitResult::operator const RooFitResult *() const
468{
469 return fNode->get<const RooFitResult>();
470}
472{
473 fNode->Draw(opt);
474}
475
476xRooNLLVar::xRooFitResult xRooNLLVar::xRooFitResult::cfit(const char *poiValues, const char *alias)
477{
478
479 // create a hypoPoint with ufit equal to this fit
480 // and poi equal to given poi
481 if (!fNll)
482 throw std::runtime_error("xRooFitResult::cfit: Cannot create cfit without nll");
483
484 // see if fit already done
485 if (alias) {
486 if (auto res = fCfits->find(alias); res != fCfits->end()) {
487 return res->second;
488 }
489 }
490 if (auto res = fCfits->find(poiValues); res != fCfits->end()) {
491 return res->second;
492 }
493
494 AutoRestorer s(*fNll->fFuncVars);
495 *fNll->fFuncVars = get()->floatParsFinal();
496 fNll->fFuncVars->assignValueOnly(get()->constPars());
497 std::unique_ptr<RooAbsCollection>(fNll->fFuncVars->selectCommon(get()->floatParsFinal()))
498 ->setAttribAll("Constant", false);
499 std::unique_ptr<RooAbsCollection>(fNll->fFuncVars->selectCommon(get()->constPars()))->setAttribAll("Constant", true);
500
501 auto hp = fNll->hypoPoint(poiValues, std::numeric_limits<double>::quiet_NaN(), xRooFit::Asymptotics::Unknown);
502 hp.fUfit = *this;
503 auto out = xRooNLLVar::xRooFitResult(std::make_shared<xRooNode>(hp.cfit_null(), fNode->fParent), fNll);
504 fCfits->insert(std::pair((alias) ? alias : poiValues, out));
505 return out;
506}
508{
509 RooRealVar *npVar = dynamic_cast<RooRealVar *>((prefit ? get()->floatParsInit() : get()->floatParsFinal()).find(np));
510 if (!npVar)
511 throw std::runtime_error("xRooFitResult::ifit: par not found");
512 return cfit(TString::Format("%s=%f", np, npVar->getVal() + (up ? npVar->getErrorHi() : npVar->getErrorLo())));
513}
514double xRooNLLVar::xRooFitResult::impact(const char *poi, const char *np, bool up, bool prefit, bool covApprox)
515{
516 if (!covApprox) {
517 // get the ifit and get the difference between the postFit poi values
518 RooRealVar *poiHat = dynamic_cast<RooRealVar *>((get()->floatParsFinal()).find(poi));
519 if (!poiHat)
520 throw std::runtime_error("xRooFitResult::impact: poi not found");
521 auto _ifit = ifit(np, up, prefit);
522 if (!_ifit)
523 throw std::runtime_error("xRooFitResult::impact: null ifit");
524 if (_ifit->status() != 0)
525 fNode->Warning("impact", "ifit status code is %d", _ifit->status());
526 return _ifit->floatParsFinal().getRealValue(poi) - poiHat->getVal();
527 } else {
528 // estimate impact from the covariance matrix ....
529 int iPoi = get()->floatParsFinal().index(poi);
530 int iNp = get()->floatParsFinal().index(np);
531 if (iPoi == -1)
532 throw std::runtime_error("xRooFitResult::impact: poi not found");
533 if (iNp == -1)
534 throw std::runtime_error("xRooFitResult::impact: np not found");
535 RooRealVar *npVar =
536 dynamic_cast<RooRealVar *>((prefit ? get()->floatParsInit() : get()->floatParsFinal()).find(np));
537 return get()->covarianceMatrix()(iPoi, iNp) / (up ? npVar->getErrorHi() : npVar->getErrorLo());
538 }
539 return std::numeric_limits<double>::quiet_NaN();
540}
541
542double xRooNLLVar::xRooFitResult::conditionalError(const char *poi, const char *nps, bool up, bool covApprox)
543{
544 // run a fit with given NPs held constant, return quadrature difference
545
546 TString npNames;
547 RooArgList vars;
548 RooAbsArg *poiVar = nullptr;
549 for (auto p : get()->floatParsFinal()) {
550 if (strcmp(p->GetName(), poi) == 0) {
551 vars.add(*p);
552 poiVar = p;
553 continue;
554 }
555 TStringToken pattern(nps, ",");
556 bool matches = false;
557 while (pattern.NextToken()) {
558 TString s(pattern);
559 if ((p->getStringAttribute("group") && s == p->getStringAttribute("group")) ||
560 TString(p->GetName()).Contains(TRegexp(s, true)) || p->getAttribute(s)) {
561 matches = true;
562 break;
563 }
564 }
565 if (matches) {
566 if (npNames.Length())
567 npNames += ",";
568 npNames += p->GetName();
569 } else {
570 vars.add(*p); // keeping in reduced cov matrix
571 }
572 }
573 if (!poiVar) {
574 throw std::runtime_error(TString::Format("Could not find poi: %s", poi));
575 }
576 if (npNames == "") {
577 fNode->Warning("conditionalError", "No parameters selected by: %s", nps);
578 return (up) ? static_cast<RooRealVar *>(poiVar)->getErrorHi() : static_cast<RooRealVar *>(poiVar)->getErrorLo();
579 }
580
581 if (covApprox) {
582 int idx = vars.index(poi);
583 return sqrt(get()->conditionalCovarianceMatrix(vars)(idx, idx));
584 }
585
586 auto _cfit = cfit(npNames.Data(), nps);
587
588 auto _poi = _cfit->floatParsFinal().find(poi);
589
590 return (up) ? static_cast<RooRealVar *>(_poi)->getErrorHi() : static_cast<RooRealVar *>(_poi)->getErrorLo();
591}
592
593RooArgList xRooNLLVar::xRooFitResult::ranknp(const char *poi, bool up, bool prefit, double approxThreshold)
594{
595
596 RooRealVar *poiHat = dynamic_cast<RooRealVar *>((get()->floatParsFinal()).find(poi));
597 if (!poiHat)
598 throw std::runtime_error("xRooFitResult::ranknp: poi not found");
599
600 std::vector<std::pair<std::string, double>> ranks;
601 // first do with the covariance approximation, since that's always available
602 for (auto par : get()->floatParsFinal()) {
603 if (par == poiHat)
604 continue;
605 ranks.emplace_back(std::pair(par->GetName(), impact(poi, par->GetName(), up, prefit, true)));
606 }
607
608 std::sort(ranks.begin(), ranks.end(), [](auto &left, auto &right) {
609 if (std::isnan(left.second) && !std::isnan(right.second))
610 return false;
611 if (!std::isnan(left.second) && std::isnan(right.second))
612 return true;
613 return fabs(left.second) > fabs(right.second);
614 });
615
616 // now redo the ones above the threshold
617 for (auto &[n, v] : ranks) {
618 if (v >= approxThreshold) {
619 try {
620 v = impact(poi, n.c_str(), up, prefit);
621 } catch (...) {
622 v = std::numeric_limits<double>::quiet_NaN();
623 };
624 }
625 }
626
627 // resort
628 std::sort(ranks.begin(), ranks.end(), [](auto &left, auto &right) {
629 if (std::isnan(left.second) && !std::isnan(right.second))
630 return false;
631 if (!std::isnan(left.second) && std::isnan(right.second))
632 return true;
633 return fabs(left.second) > fabs(right.second);
634 });
635
636 RooArgList out;
637 out.setName("rankings");
638 for (auto &[n, v] : ranks) {
639 out.addClone(*get()->floatParsFinal().find(n.c_str()));
640 auto vv = static_cast<RooRealVar *>(out.at(out.size() - 1));
641 vv->setVal(v);
642 vv->removeError();
643 vv->removeRange();
644 }
645 return out;
646}
647
648xRooNLLVar::xRooFitResult xRooNLLVar::minimize(const std::shared_ptr<ROOT::Fit::FitConfig> &_config)
649{
650 auto &nll = *get();
651 auto out = xRooFit::minimize(nll, (_config) ? _config : fitConfig(), fOpts);
652 // add any pars that are const here that aren't in constPars list because they may have been
653 // const-optimized and their values cached with the dataset, so if subsequently floated the
654 // nll wont evaluate correctly
655 // fConstVars.reset( fFuncVars->selectByAttrib("Constant",true) );
656
657 // before returning, flag which of the constPars were actually global observables
658 if (out) {
659 const_cast<RooArgList &>(out->constPars()).setAttribAll("global", false);
660 if (fGlobs)
661 std::unique_ptr<RooAbsCollection>(out->constPars().selectCommon(*fGlobs))->setAttribAll("global", true);
662 }
663 return xRooFitResult(std::make_shared<xRooNode>(out, fPdf), std::make_shared<xRooNLLVar>(*this));
664}
665
666std::shared_ptr<ROOT::Fit::FitConfig> xRooNLLVar::fitConfig()
667{
668 if (!fFitConfig)
670 return fFitConfig;
671}
672
674{
675 if (auto conf = fitConfig(); conf)
676 return const_cast<ROOT::Math::IOptions *>(conf->MinimizerOptions().ExtraOptions());
677 return nullptr;
678}
679
680double xRooNLLVar::getEntryVal(size_t entry) const
681{
682 auto _data = data();
683 if (!_data)
684 return 0;
685 if (size_t(_data->numEntries()) <= entry)
686 return 0;
687 auto _pdf = pdf();
688 *std::unique_ptr<RooAbsCollection>(_pdf->getObservables(_data)) = *_data->get(entry);
689 // if (auto s = dynamic_cast<RooSimultaneous*>(_pdf.get());s) return
690 // -_data->weight()*s->getPdf(s->indexCat().getLabel())->getLogVal(_data->get());
691 return -_data->weight() * _pdf->getLogVal(_data->get());
692}
693
694std::set<std::string> xRooNLLVar::binnedChannels() const
695{
696 std::set<std::string> out;
697
698 auto binnedOpt = dynamic_cast<RooCmdArg *>(fOpts->find("Binned")); // the binned option, if explicitly specified
699
700 if (auto s = dynamic_cast<RooSimultaneous *>(pdf().get())) {
701 xRooNode simPdf(*s);
702 bool allChannels = true;
703 for (auto c : simPdf.bins()) {
704 // see if there's a RooRealSumPdf in the channel - if there is, if it has BinnedLikelihood set
705 // then assume is a BinnedLikelihood channel
706 RooArgSet nodes;
707 c->get<RooAbsArg>()->treeNodeServerList(&nodes, nullptr, true, false);
708 bool isBinned = false;
709 for (auto a : nodes) {
710 if (a->InheritsFrom("RooRealSumPdf") &&
711 ((binnedOpt && binnedOpt->getInt(0)) || (!binnedOpt && a->getAttribute("BinnedLikelihood")))) {
712 TString chanName(c->GetName());
713 out.insert(chanName(chanName.Index("=") + 1, chanName.Length()).Data());
714 isBinned = true;
715 break;
716 }
717 }
718 if (!isBinned) {
719 allChannels = false;
720 }
721 }
722 if (allChannels) {
723 out.clear();
724 out.insert("*");
725 }
726 } else {
727 RooArgSet nodes;
728 pdf()->treeNodeServerList(&nodes, nullptr, true, false);
729 for (auto a : nodes) {
730 if (a->InheritsFrom("RooRealSumPdf") &&
731 ((binnedOpt && binnedOpt->getInt(0)) || (!binnedOpt && a->getAttribute("BinnedLikelihood")))) {
732 out.insert("*");
733 break;
734 }
735 }
736 }
737 return out;
738}
739
740double xRooNLLVar::getEntryBinWidth(size_t entry) const
741{
742
743 auto _data = data();
744 if (!_data)
745 return 0;
746 if (size_t(_data->numEntries()) <= entry)
747 return 0;
748 auto _pdf = pdf().get();
749 std::unique_ptr<RooAbsCollection> _robs(_pdf->getObservables(_data->get()));
750 *_robs = *_data->get(entry); // only set robs
751 if (auto s = dynamic_cast<RooSimultaneous *>(_pdf); s) {
752 _pdf = s->getPdf(s->indexCat().getCurrentLabel());
753 }
754 double volume = 1.;
755 for (auto o : *_robs) {
756
757 if (auto a = dynamic_cast<RooAbsRealLValue *>(o);
758 a && _pdf->dependsOn(*a)) { // dependsOn check needed until ParamHistFunc binBoundaries method fixed
759 std::unique_ptr<std::list<double>> bins(
760 _pdf->binBoundaries(*a, -std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity()));
761 if (bins) {
762 double lowEdge = -std::numeric_limits<double>::infinity();
763 for (auto b : *bins) {
764 if (b > a->getVal()) {
765 volume *= (b - lowEdge);
766 break;
767 }
768 lowEdge = b;
769 }
770 }
771 }
772 }
773
774 return volume;
775}
776
778{
779 // for each global observable in the dataset, determine which constraint term is associated to it
780 // and given its type, add the necessary saturated term...
781
782 double out = 0;
783
784 if (!fGlobs)
785 return 0;
786
787 auto cTerm = constraintTerm();
788 if (!cTerm)
789 return 0;
790
791 for (auto c : cTerm->list()) {
792 if (std::string(c->ClassName()) == "RooAbsPdf") {
793 // in ROOT 6.32 the constraintTerm is full of RooNormalizedPdfs which aren't public
794 // in that case use the first server
795 c = c->servers()[0];
796 }
797 if (auto gaus = dynamic_cast<RooGaussian *>(c)) {
798 auto v = dynamic_cast<RooAbsReal *>(fGlobs->find(gaus->getX().GetName()));
799 if (!v) {
800 v = dynamic_cast<RooAbsReal *>(fGlobs->find(
801 gaus->getMean().GetName())); // shouldn't really happen but does for at least ws made by pyhf
802 }
803 if (!v)
804 continue;
805 out -= std::log(ROOT::Math::gaussian_pdf(v->getVal(), gaus->getSigma().getVal(), v->getVal()));
806 } else if (auto pois = dynamic_cast<RooPoisson *>(c)) {
807 auto v = dynamic_cast<RooAbsReal *>(fGlobs->find(pois->getX().GetName()));
808 if (!v)
809 continue;
810 out -= std::log(TMath::Poisson(v->getVal(), v->getVal()));
811 }
812 }
813
814 return out;
815}
816
817double xRooNLLVar::ndof() const
818{
819 return data()->numEntries() + (fFuncGlobs ? fFuncGlobs->size() : 0) -
820 std::unique_ptr<RooAbsCollection>(pars()->selectByAttrib("Constant", false))->size();
821}
822
823double xRooNLLVar::pgof() const
824{
825 // note that if evaluating this for a single channel, until 6.30 is available if you are using Binned mode the pdf
826 // will need to be part of a Simultaneous
827 return TMath::Prob(2. * (get()->getVal() - saturatedVal()), ndof());
828}
829
831{
832 // need to count number of floating unconstrained parameters
833 // which are floating parameters not featured in the constraintTerm
834 std::unique_ptr<RooAbsCollection> _floats(pars()->selectByAttrib("Constant", false));
835 if (auto _constraintTerm = constraintTerm()) {
836 _floats->remove(*std::unique_ptr<RooAbsCollection>(_constraintTerm->getVariables()));
837 }
838 return data()->numEntries() - _floats->size();
839}
840
842{
843 // using totVal - constraintTerm while new evalbackend causes mainTerm() to return nullptr
844 double val = get()->getVal();
845 if (auto _constraintTerm = constraintTerm()) {
846 val -= _constraintTerm->getVal();
847 }
848
849 return TMath::Prob(2. * (val - saturatedMainTerm()), mainTermNdof());
850}
851
853{
855}
856
858{
859
860 // Use this term to create a goodness-of-fit metric, which is approx chi2 distributed with numEntries (data) d.o.f:
861 // prob = TMath::Prob( 2.*(nll.mainTerm()->getVal() - nll.saturatedNllTerm()), nll.data()->numEntries() )
862
863 // note that need to construct nll with explicit Binned(1 or 0) option otherwise will pick up nll eval
864 // from attributes in model already, so many get binned mainTerm eval when thinking not binned because didnt specify
865 // Binned(1)
866
867 auto _data = data();
868 if (!_data)
869 return std::numeric_limits<double>::quiet_NaN();
870
871 std::set<std::string> _binnedChannels = binnedChannels();
872
873 // for binned case each entry is: -(-N + Nlog(N) - std::lgamma(N+1))
874 // for unbinned case each entry is: -(N*log(N/(sumN*binW))) = -N*logN + N*log(sumN) + N*log(binW)
875 // but unbinned gets extendedTerm = sumN - sumN*log(sumN)
876 // so resulting sum is just sumN - sum[ N*logN - N*log(binW) ]
877 // which is the same as the binned case without the LnGamma part and with the extra sum[N*log(binW)] part
878
879 const RooAbsCategoryLValue *cat = (dynamic_cast<RooSimultaneous *>(pdf().get()))
880 ? &dynamic_cast<RooSimultaneous *>(pdf().get())->indexCat()
881 : nullptr;
882
883 double out = _data->sumEntries();
884 for (int i = 0; i < _data->numEntries(); i++) {
885 _data->get(i);
886 double w = _data->weight();
887 out -= w * std::log(w);
888 if (_binnedChannels.count("*")) {
889 out += std::lgamma(w + 1);
890 } else if (_binnedChannels.empty()) {
891 out += w * std::log(getEntryBinWidth(i));
892 } else if (cat) {
893 // need to determine which channel we are in for this entry to decide if binned or unbinned active
894 if (_binnedChannels.count(_data->get()->getCatLabel(cat->GetName()))) {
895 out += std::lgamma(w + 1);
896 } else {
897 out += w * std::log(getEntryBinWidth(i));
898 }
899 } else {
900 throw std::runtime_error("Cannot determine category of RooSimultaneous pdf");
901 }
902 }
903
904 out += simTerm();
905
906 return out;
907}
908
909std::shared_ptr<RooArgSet> xRooNLLVar::pars(bool stripGlobalObs) const
910{
911 auto out = std::shared_ptr<RooArgSet>(get()->getVariables());
912 if (stripGlobalObs && fGlobs) {
913 out->remove(*fGlobs, true, true);
914 }
915 return out;
916}
917
918TObject *
919xRooNLLVar::Scan(const char *scanPars, const std::vector<std::vector<double>> &coords, const RooArgList &profilePars)
920{
921 return Scan(*std::unique_ptr<RooAbsCollection>(get()->getVariables()->selectByName(scanPars)), coords, profilePars);
922}
923
924TObject *xRooNLLVar::Scan(const RooArgList &scanPars, const std::vector<std::vector<double>> &coords,
925 const RooArgList &profilePars)
926{
927
928 if (scanPars.size() > 2 || scanPars.empty())
929 return nullptr;
930
931 TGraph2D *out2d = (scanPars.size() == 2) ? new TGraph2D() : nullptr;
932 TGraph *out1d = (out2d) ? nullptr : new TGraph();
933 TNamed *out = (out2d) ? static_cast<TNamed *>(out2d) : static_cast<TNamed *>(out1d);
934 out->SetName(get()->GetName());
935 out->SetTitle(TString::Format("%s;%s%s%s", get()->GetTitle(), scanPars.first()->GetTitle(), out2d ? ";" : "",
936 out2d ? scanPars.at(1)->GetTitle() : ""));
937
938 std::unique_ptr<RooAbsCollection> funcVars(get()->getVariables());
939 AutoRestorer snap(*funcVars);
940
941 for (auto &coord : coords) {
942 if (coord.size() != scanPars.size()) {
943 throw std::runtime_error("Invalid coordinate");
944 }
945 for (size_t i = 0; i < coord.size(); i++) {
946 static_cast<RooAbsRealLValue &>(scanPars[i]).setVal(coord[i]);
947 }
948
949 if (profilePars.empty()) {
950 // just evaluate
951 if (out2d) {
952 out2d->SetPoint(out2d->GetN(), coord[0], coord[1], get()->getVal());
953 } else {
954 out1d->SetPoint(out1d->GetN(), coord[0], get()->getVal());
955 }
956 }
957 }
958
959 return out;
960}
961
963{
964 TString sOpt(opt);
965
966 auto _pars = pars();
967
968 if (sOpt == "sensitivity") {
969
970 // will make a plot of DeltaNLL
971 }
972
973 if (sOpt == "floating") {
974 // start scanning floating pars
975 auto floats = std::unique_ptr<RooAbsCollection>(_pars->selectByAttrib("Constant", false));
976 TVirtualPad *pad = gPad;
977 if (!pad) {
979 pad = gPad;
980 }
982 gr->SetName("multigraph");
983 gr->SetTitle(TString::Format("%s;Normalized Parameter Value;#Delta NLL", get()->GetTitle()));
984 /*((TPad*)pad)->DivideSquare(floats->size());
985 int i=0;
986 for(auto a : *floats) {
987 i++;
988 pad->cd(i);
989 Draw(a->GetName());
990 }*/
991 return;
992 }
993
994 RooArgList vars;
995 TStringToken pattern(sOpt, ":");
996 while (pattern.NextToken()) {
997 TString s(pattern);
998 if (auto a = _pars->find(s); a)
999 vars.add(*a);
1000 }
1001
1002 if (vars.size() == 1) {
1003 TGraph *out = new TGraph;
1004 out->SetBit(kCanDelete);
1005 TGraph *bad = new TGraph;
1006 bad->SetBit(kCanDelete);
1007 bad->SetMarkerColor(kRed);
1008 bad->SetMarkerStyle(5);
1009 TMultiGraph *gr = (gPad) ? dynamic_cast<TMultiGraph *>(gPad->GetPrimitive("multigraph")) : nullptr;
1010 bool normRange = false;
1011 if (!gr) {
1012 gr = new TMultiGraph;
1013 gr->Add(out, "LP");
1015 } else {
1016 normRange = true;
1017 }
1018 out->SetName(get()->GetName());
1019 gr->SetTitle(TString::Format("%s;%s;#Delta NLL", get()->GetTitle(), vars.at(0)->GetTitle()));
1020 // scan outwards from current value towards limits
1021 auto v = dynamic_cast<RooRealVar *>(vars.at(0));
1022 double low = v->getVal();
1023 double high = low;
1024 double step = (v->getMax() - v->getMin()) / 100;
1025 double init = v->getVal();
1026 double initVal = func()->getVal();
1027 // double xscale = (normRange) ? (2.*(v->getMax() - v->getMin())) : 1.;
1028 auto currTime = std::chrono::steady_clock::now();
1029 while (out->GetN() < 100 && (low > v->getMin() || high < v->getMax())) {
1030 if (out->GetN() == 0) {
1031 out->SetPoint(out->GetN(), low, 0);
1032 low -= step;
1033 high += step;
1034 if (!normRange) {
1035 gr->Draw("A");
1036 gPad->SetGrid();
1037 }
1038 continue;
1039 }
1040 if (low > v->getMin()) {
1041 v->setVal(low);
1042 auto _v = func()->getVal();
1043 if (std::isnan(_v) || std::isinf(_v)) {
1044 if (bad->GetN() == 0)
1045 gr->Add(bad, "P");
1046 bad->SetPoint(bad->GetN(), low, out->GetPointY(0));
1047 } else {
1048 out->SetPoint(out->GetN(), low, _v - initVal);
1049 }
1050 low -= step;
1051 }
1052 if (high < v->getMax()) {
1053 v->setVal(high);
1054 auto _v = func()->getVal();
1055 if (std::isnan(_v) || std::isinf(_v)) {
1056 if (bad->GetN() == 0)
1057 gr->Add(bad, "P");
1058 bad->SetPoint(bad->GetN(), high, out->GetPointY(0));
1059 } else {
1060 out->SetPoint(out->GetN(), high, _v - initVal);
1061 }
1062 high += step;
1063 }
1064 out->Sort();
1065 // should only do processEvents once every second in case using x11 (which is slow)
1066 gPad->Modified();
1067 if (std::chrono::steady_clock::now() - currTime > std::chrono::seconds(1)) {
1068 currTime = std::chrono::steady_clock::now();
1069 gPad->Update();
1071 }
1072 }
1073 // add arrow to show current par value
1074 TArrow a;
1075 a.DrawArrow(init, 0, init, -0.1);
1076 gPad->Update();
1077#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
1078 gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
1079#endif
1081 v->setVal(init);
1082 } else {
1083 Error("Draw", "Name a parameter to scan over: Draw(<name>) , choose from: %s",
1084 _pars->empty() ? "" : _pars->contentsString().c_str());
1085 }
1086}
1087
1088std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> xRooNLLVar::getData() const
1089{
1090 return std::make_pair(fData, fGlobs);
1091}
1092
1094{
1095 if (data.fComp && !data.get<RooAbsData>()) {
1096 return false;
1097 }
1098 return setData(std::dynamic_pointer_cast<RooAbsData>(data.fComp),
1099 std::shared_ptr<const RooAbsCollection>(data.globs().argList().snapshot()));
1100}
1101
1102bool xRooNLLVar::setData(const std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> &_data)
1103{
1104
1105 if (fData == _data.first && fGlobs == _data.second)
1106 return true;
1107
1108 auto _globs = fGlobs; // done to keep globs alive while NLL might still be alive.
1109
1110 auto _dglobs = (_data.second) ? _data.second
1111 : std::shared_ptr<const RooAbsCollection>(_data.first->getGlobalObservables(),
1112 [](const RooAbsCollection *) {});
1113
1114 if (fGlobs && !(fGlobs->empty() && !_dglobs) && _data.first &&
1115 fGlobs != _dglobs) { // second condition allows for no globs being a nullptr, third allow globs to remain if
1116 // nullifying data
1117 if (!_dglobs)
1118 throw std::runtime_error("Missing globs");
1119 // ignore 'extra' globs
1120 RooArgSet s;
1121 s.add(*fGlobs);
1122 std::unique_ptr<RooAbsCollection> _actualGlobs(fPdf->getObservables(s));
1123 RooArgSet s2;
1124 s2.add(*_dglobs);
1125 std::unique_ptr<RooAbsCollection> _actualGlobs2(fPdf->getObservables(s2));
1126 if (!_actualGlobs->equals(*_actualGlobs2)) {
1127 RooArgSet rC;
1128 rC.add(*_actualGlobs2);
1129 rC.remove(*std::unique_ptr<RooAbsCollection>(rC.selectCommon(*_actualGlobs)));
1130 TString r = (!rC.empty()) ? rC.contentsString() : "";
1131 RooArgSet lC;
1132 lC.add(*_actualGlobs);
1133 lC.remove(*std::unique_ptr<RooAbsCollection>(lC.selectCommon(*_actualGlobs2)));
1134 TString l = (!lC.empty()) ? lC.contentsString() : "";
1135 throw std::runtime_error(TString::Format("globs mismatch: adding %s removing %s", r.Data(), l.Data()));
1136 }
1137 fGlobs = _dglobs;
1138 }
1139
1140 if (!std::shared_ptr<RooAbsReal>::get()) {
1141 fData = _data.first;
1142 return true; // not loaded yet so nothing to do
1143 }
1144
1145 try {
1147 // happens when using MP need to rebuild the nll instead
1148 // also happens if there's no mainTerm(), which is the case in 6.32 where RooNLLVar is partially deprecated
1149 AutoRestorer snap(*fFuncVars);
1150 // ensure the const state is back where it was at nll construction time;
1151 fFuncVars->setAttribAll("Constant", false);
1152 fConstVars->setAttribAll("Constant", true);
1153 std::shared_ptr<RooAbsData> __data = fData; // do this just to keep fData alive while killing previous NLLVar
1154 // (can't kill data while NLL constructed with it)
1155 fData = _data.first;
1156 reinitialize();
1157 return true;
1158 }
1159 bool out = false;
1160 if (_data.first) {
1161 if (_data.first->getGlobalObservables()) {
1162 // replace in all terms
1163 get()->setData(*_data.first, false);
1164 } else {
1165 // replace just in mainTerm ... note to self: why not just replace in all like above? should test!
1166 out = mainTerm()->setData(*_data.first, false /* clone data? */);
1167 }
1168 } else {
1169 reset();
1170 }
1171 fData = _data.first;
1172 return out;
1173 } catch (std::runtime_error &) {
1174 // happens when using MP need to rebuild the nll instead
1175 // also happens if there's no mainTerm(), which is the case in 6.32 where RooNLLVar is partially deprecated
1176 AutoRestorer snap(*fFuncVars);
1177 // ensure the const state is back where it was at nll construction time;
1178 fFuncVars->setAttribAll("Constant", false);
1179 fConstVars->setAttribAll("Constant", true);
1180 std::shared_ptr<RooAbsData> __data = fData; // do this just to keep fData alive while killing previous NLLVar
1181 // (can't kill data while NLL constructed with it)
1182 fData = _data.first;
1183 reinitialize();
1184 return true;
1185 }
1186 throw std::runtime_error("Unable to setData");
1187}
1188
1189std::shared_ptr<RooAbsReal> xRooNLLVar::func() const
1190{
1191 if (!(*this)) {
1192 const_cast<xRooNLLVar *>(this)->reinitialize();
1193 } else if (auto f = std::unique_ptr<RooAbsCollection>(fConstVars->selectByAttrib("Constant", false)); !f->empty()) {
1194 // have to reinitialize if const par values have changed - const optimization forces this
1195 // TODO: currently changes to globs also triggers this since the vars includes globs (vars are the non-obs pars)
1196 // std::cout << "Reinitializing because of change of const parameters:" << f->contentsString() << std::endl;
1197 const_cast<xRooNLLVar *>(this)->reinitialize();
1198
1199 // note ... it may be sufficient here to do:
1200 // nll.constOptimizeTestStatistic(RooAbsArg::ConfigChange, constOptimize>1 /* do tracking too if >1 */); //
1201 // trigger a re-evaluate of which nodes to cache-and-track nll.constOptimizeTestStatistic(RooAbsArg::ValueChange,
1202 // constOptimize>1); // update the cache values -- is this needed??
1203 // this forces the optimization to be redone
1204 // for now leave as a reinitialize though, until had a chance to test this properly
1205 }
1206 if (fGlobs && fFuncGlobs) {
1207 *fFuncGlobs = *fGlobs;
1208 fFuncGlobs->setAttribAll("Constant", true);
1209 }
1210 return *this;
1211}
1212
1214{
1215 fOpts->Add(opt.Clone(nullptr));
1216 if (std::shared_ptr<RooAbsReal>::get()) {
1217 reinitialize(); // do this way to keep name of nll if user set
1218 } else {
1219 reset(); // will trigger reinitialize
1220 }
1221}
1222
1224{
1225 auto _nll = mainTerm();
1226 if (!_nll)
1227 return fData.get();
1228 RooAbsData *out = &_nll->data();
1229 if (!out)
1230 return fData.get();
1231 return out;
1232}
1233
1235{
1236 auto _func = func();
1237 if (auto a = dynamic_cast<RooNLLVar *>(_func.get()); a)
1238 return a;
1239 for (auto s : _func->servers()) {
1240 if (auto a = dynamic_cast<RooNLLVar *>(s); a)
1241 return a;
1242 }
1243 return nullptr;
1244}
1245
1247{
1248 // returns Nexp - Nobs*log(Nexp)
1249 return fPdf->extendedTerm(fData->sumEntries(), fData->get());
1250}
1251
1253{
1254 if (auto s = dynamic_cast<RooSimultaneous *>(fPdf.get()); s) {
1255 return fData->sumEntries() * log(1.0 * (s->servers().size() - 1)); // one of the servers is the cat
1256 }
1257 return 0;
1258}
1259
1261{
1262 // this is only relevant if BinnedLikelihood active
1263 // = sum[ N_i! ] since LnGamma(N_i+1) ~= N_i!
1264 // need to also subtract off sum[ N_i*log(width_i) ] in order to have formula: binnedLL = unbinnedLL + binnedDataTerm
1265 // note this is 0 if all the bin widths are 1
1266 double out = 0;
1267 for (int i = 0; i < fData->numEntries(); i++) {
1268 fData->get(i);
1269 out += std::lgamma(fData->weight() + 1) - fData->weight() * std::log(getEntryBinWidth(i));
1270 }
1271
1272 return out;
1273}
1274
1276{
1277 auto _func = func();
1278 if (auto a = dynamic_cast<RooConstraintSum *>(_func.get()); a)
1279 return a;
1280 for (auto s : _func->servers()) {
1281 if (auto a = dynamic_cast<RooConstraintSum *>(s); a)
1282 return a;
1283 // allow one more depth to support 6.32 (where sum is hidden inside the first server)
1284 for (auto s2 : s->servers()) {
1285 if (auto a2 = dynamic_cast<RooConstraintSum *>(s2); a2)
1286 return a2;
1287 }
1288 }
1289 return nullptr;
1290}
1291
1292/*xRooNLLVar::operator RooAbsReal &() const {
1293 // this works in c++ but not in python
1294 std::cout << "implicit conversion" << std::endl;
1295 return *fFunc;
1296}*/
1297
1299{
1300 TString sWhat(what);
1301 sWhat.ToLower();
1302 bool doTS = sWhat.Contains("ts");
1303 bool doCLs = sWhat.Contains("pcls");
1304 bool doNull = sWhat.Contains("pnull");
1305 bool doAlt = sWhat.Contains("palt");
1306 double nSigma = (sWhat.Contains("exp"))
1307 ? (TString(sWhat(sWhat.Index("exp") + 3, sWhat.Index(" ", sWhat.Index("exp")) == -1
1308 ? sWhat.Length()
1309 : sWhat.Index(" ", sWhat.Index("exp"))))
1310 .Atof())
1311 : std::numeric_limits<double>::quiet_NaN();
1312
1313 bool toys = sWhat.Contains("toys");
1314
1315 // bool asymp = sWhat.Contains("asymp");
1316
1317 bool readOnly = sWhat.Contains("readonly");
1318
1319 if (!readOnly) {
1320 if (toys) {
1321 sigma_mu(); // means we will be able to evaluate the asymptotic values too
1322 }
1323 // only add toys if actually required
1324 if (getVal(sWhat + " readonly").second != 0) {
1325 if (sWhat.Contains("toys=")) {
1326 // extract number of toys required ... format is "nullToys.altToysFraction" if altToysFraction=0 then use
1327 // same for both
1328 size_t nToys = TString(sWhat(sWhat.Index("toys=") + 5, sWhat.Length())).Atoi();
1329 size_t nToysAlt = (TString(sWhat(sWhat.Index("toys=") + 5, sWhat.Length())).Atof() - nToys) * nToys;
1330 if (nToysAlt == 0)
1331 nToysAlt = nToys;
1332 if (nullToys.size() < nToys) {
1333 addNullToys(nToys - nullToys.size());
1334 }
1335 if (altToys.size() < nToysAlt) {
1336 addAltToys(nToysAlt - altToys.size());
1337 }
1338 } else if (doCLs && toys) {
1339 // auto toy-generating for limits .. do in blocks of 100
1340 addCLsToys(100, 0, 0.05, nSigma);
1341 }
1342 }
1343 }
1344
1345 struct RestoreNll {
1346 RestoreNll(std::shared_ptr<xRooNLLVar> &v, bool r) : rr(r), var(v)
1347 {
1348 if (rr && var && var->get()) {
1349 _readOnly = var->get()->getAttribute("readOnly");
1350 var->get()->setAttribute("readOnly", rr);
1351 } else {
1352 rr = false;
1353 }
1354 };
1355 ~RestoreNll()
1356 {
1357 if (rr)
1358 var->get()->setAttribute("readOnly", _readOnly);
1359 };
1360
1361 bool rr = false;
1362 bool _readOnly = false;
1363
1364 std::shared_ptr<xRooNLLVar> &var;
1365 };
1366
1367 RestoreNll rest(nllVar, readOnly);
1368
1369 if (doTS)
1370 return (toys) ? ts_toys(nSigma) : ts_asymp(nSigma);
1371 if (doNull)
1372 return (toys) ? pNull_toys(nSigma) : pNull_asymp(nSigma);
1373 if (doAlt)
1374 return (toys) ? pAlt_toys(nSigma) : pAlt_asymp(nSigma);
1375 if (doCLs)
1376 return (toys) ? pCLs_toys(nSigma) : pCLs_asymp(nSigma);
1377
1378 throw std::runtime_error(std::string("Unknown: ") + what);
1379}
1380
1382{
1383 RooArgList out;
1384 out.setName("poi");
1385 out.add(*std::unique_ptr<RooAbsCollection>(coords->selectByAttrib("poi", true)));
1386 return out;
1387}
1388
1390{
1391 RooArgList out;
1392 out.setName("alt_poi");
1393 out.addClone(*std::unique_ptr<RooAbsCollection>(coords->selectByAttrib("poi", true)));
1394 for (auto a : out) {
1395 auto v = dynamic_cast<RooAbsRealLValue *>(a);
1396 if (!v)
1397 continue;
1398 if (auto s = a->getStringAttribute("altVal"); s && strlen(s)) {
1399 v->setVal(TString(s).Atof());
1400 } else {
1401 v->setVal(std::numeric_limits<double>::quiet_NaN());
1402 }
1403 }
1404 return out;
1405}
1406
1408{
1409 auto &me = const_cast<xRooHypoPoint &>(*this);
1410 int out = 0;
1411 if (me.ufit(true) && !allowedStatusCodes.count(me.ufit(true)->status()))
1412 out += 1;
1413 if (me.cfit_null(true) && !allowedStatusCodes.count(me.cfit_null(true)->status()))
1414 out += 1 << 1;
1415 if (me.cfit_alt(true) && !allowedStatusCodes.count(me.cfit_alt(true)->status()))
1416 out += 1 << 2;
1417 if (me.asimov(true))
1418 out += me.asimov(true)->status() << 3;
1419 return out;
1420}
1421
1423{
1424 std::cout << "POI: " << const_cast<xRooHypoPoint *>(this)->poi().contentsString()
1425 << " , null: " << dynamic_cast<RooAbsReal *>(const_cast<xRooHypoPoint *>(this)->poi().first())->getVal()
1426 << " , alt: "
1427 << dynamic_cast<RooAbsReal *>(const_cast<xRooHypoPoint *>(this)->alt_poi().first())->getVal();
1428 std::cout << " , pllType: " << fPllType << std::endl;
1429
1430 std::cout << " - ufit: ";
1431 if (fUfit) {
1432 std::cout << fUfit->GetName() << " " << fUfit->minNll() << " (status=" << fUfit->status() << ") ("
1433 << const_cast<xRooHypoPoint *>(this)->mu_hat().GetName()
1434 << "_hat: " << const_cast<xRooHypoPoint *>(this)->mu_hat().getVal() << " +/- "
1435 << const_cast<xRooHypoPoint *>(this)->mu_hat().getError() << ")" << std::endl;
1436 } else {
1437 std::cout << "Not calculated" << std::endl;
1438 }
1439 std::cout << " - null cfit: ";
1440 if (fNull_cfit) {
1441 std::cout << fNull_cfit->GetName() << " " << fNull_cfit->minNll() << " (status=" << fNull_cfit->status() << ")";
1442 } else {
1443 std::cout << "Not calculated";
1444 }
1445 if (!std::isnan(dynamic_cast<RooAbsReal *>(const_cast<xRooHypoPoint *>(this)->alt_poi().first())->getVal())) {
1446 std::cout << std::endl << " - alt cfit: ";
1447 if (fAlt_cfit) {
1448 std::cout << fAlt_cfit->GetName() << " " << fAlt_cfit->minNll() << " (status=" << fAlt_cfit->status() << ")"
1449 << std::endl;
1450 } else {
1451 std::cout << "Not calculated" << std::endl;
1452 }
1453 std::cout << " sigma_mu: ";
1454 const_cast<xRooHypoPoint *>(this)->asimov(true); // will trigger construction of fAsimov hypoPoint if possible
1455 if (!fAsimov || !fAsimov->fUfit || !fAsimov->fNull_cfit) {
1456 std::cout << "Not calculated";
1457 } else {
1458 std::cout << const_cast<xRooHypoPoint *>(this)->sigma_mu().first << " +/- "
1459 << const_cast<xRooHypoPoint *>(this)->sigma_mu().second;
1460 }
1461 if (fAsimov) {
1462 std::cout << std::endl;
1463 std::cout << " - asimov ufit: ";
1464 if (fAsimov->fUfit) {
1465 std::cout << fAsimov->fUfit->GetName() << " " << fAsimov->fUfit->minNll()
1466 << " (status=" << fAsimov->fUfit->status() << ")";
1467 } else {
1468 std::cout << "Not calculated";
1469 }
1470 std::cout << std::endl << " - asimov null cfit: ";
1471 if (fAsimov->fNull_cfit) {
1472 std::cout << fAsimov->fNull_cfit->GetName() << " " << fAsimov->fNull_cfit->minNll()
1473 << " (status=" << fAsimov->fNull_cfit->status() << ")";
1474 } else {
1475 std::cout << "Not calculated";
1476 }
1477 }
1478 std::cout << std::endl;
1479 } else {
1480 std::cout << std::endl;
1481 }
1482 if (fGenFit)
1483 std::cout << " - genFit: " << fGenFit->GetName() << std::endl;
1484 if (!nullToys.empty() || !altToys.empty()) {
1485 std::cout << " * null toys: " << nullToys.size();
1486 size_t firstToy = 0;
1487 while (firstToy < nullToys.size() && std::isnan(std::get<1>(nullToys[firstToy])))
1488 firstToy++;
1489 if (firstToy > 0)
1490 std::cout << " [ of which " << firstToy << " are bad]";
1491 std::cout << " , alt toys: " << altToys.size();
1492 firstToy = 0;
1493 while (firstToy < altToys.size() && std::isnan(std::get<1>(altToys[firstToy])))
1494 firstToy++;
1495 if (firstToy > 0)
1496 std::cout << " [ of which " << firstToy << " are bad]";
1497 std::cout << std::endl;
1498 }
1499 // std::cout << " nllVar: " << nllVar << std::endl;
1500}
1501
1503{
1504 if (ufit()) {
1505 auto var = dynamic_cast<RooRealVar *>(ufit()->floatParsFinal().find(fPOIName()));
1506 if (var) {
1507 return *var;
1508 } else {
1509 throw std::runtime_error(TString::Format("Cannot find POI: %s", fPOIName()));
1510 }
1511 }
1512 throw std::runtime_error("Unconditional fit unavailable");
1513}
1514
1515std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> xRooNLLVar::xRooHypoPoint::data()
1516{
1517 if (fData.first)
1518 return fData;
1519 if (fGenFit && isExpected) {
1520 // std::cout << "Generating asimov" << std::endl;poi().Print("v");
1521 fData = xRooFit::generateFrom(*nllVar->fPdf, *fGenFit, true);
1522 }
1523 return fData;
1524}
1525
1526xRooNLLVar::xRooHypoPoint::xRooHypoPoint(std::shared_ptr<RooStats::HypoTestResult> htr, const RooAbsCollection *_coords)
1527 : hypoTestResult(htr)
1528{
1529 if (hypoTestResult) {
1530 // load the pllType
1531 fPllType =
1532 xRooFit::Asymptotics::PLLType(hypoTestResult->GetFitInfo()->getGlobalObservables()->getCatIndex("pllType"));
1533 isExpected = hypoTestResult->GetFitInfo()->getGlobalObservables()->getRealValue("isExpected");
1534
1535 // load the toys
1536 auto toys = hypoTestResult->GetNullDetailedOutput();
1537 if (toys) {
1538 // load coords from the nullDist globs list
1539 if (toys->getGlobalObservables()) {
1540 coords = std::shared_ptr<RooAbsCollection>(toys->getGlobalObservables()->snapshot());
1541 }
1542 for (int i = 0; i < toys->numEntries(); i++) {
1543 auto toy = toys->get(i);
1544 nullToys.emplace_back(
1545 std::make_tuple(int(toy->getRealValue("seed")), toy->getRealValue("ts"), toys->weight()));
1546 }
1547 }
1548 toys = hypoTestResult->GetAltDetailedOutput();
1549 if (toys) {
1550 for (int i = 0; i < toys->numEntries(); i++) {
1551 auto toy = toys->get(i);
1552 altToys.emplace_back(
1553 std::make_tuple(int(toy->getRealValue("seed")), toy->getRealValue("ts"), toys->weight()));
1554 }
1555 }
1556 }
1557 if (!coords && _coords)
1558 coords.reset(_coords->snapshot());
1559}
1560
1561std::shared_ptr<xRooNLLVar::xRooHypoPoint> xRooNLLVar::xRooHypoPoint::asimov(bool readOnly)
1562{
1563
1564 if (!fAsimov && (nllVar || hypoTestResult)) {
1565 auto theFit = (!fData.first && fGenFit && !isExpected)
1566 ? fGenFit
1567 : cfit_alt(readOnly); // first condition allows genFit to be used as the altFit *if* the data is
1568 // entirely absent, provided not expected data because we postpone data
1569 // creation til later in that case (see below)
1570 if (!theFit || allowedStatusCodes.find(theFit->status()) == allowedStatusCodes.end())
1571 return fAsimov;
1572 fAsimov = std::make_shared<xRooHypoPoint>(*this);
1573 fAsimov->coords.reset(fAsimov->coords->snapshot()); // create a copy so can remove the physical range below
1574 fAsimov->hypoTestResult.reset();
1575 fAsimov->fPllType = xRooFit::Asymptotics::TwoSided;
1576 for (auto p : fAsimov->poi()) {
1577 // dynamic_cast<RooRealVar *>(p)->removeRange("physical"); -- can't use this as will modify shared property
1578 if (auto v = dynamic_cast<RooRealVar *>(p)) {
1579 v->deleteSharedProperties(); // effectively removes all custom ranges
1580 }
1581 }
1582
1583 fAsimov->nullToys.clear();
1584 fAsimov->altToys.clear();
1585 fAsimov->fUfit = retrieveFit(3);
1586 fAsimov->fNull_cfit = retrieveFit(4);
1587 fAsimov->fAlt_cfit.reset();
1588 fAsimov->fData =
1589 std::make_pair(nullptr, nullptr); // postpone generating expected data until we definitely need it
1590 fAsimov->fGenFit = theFit;
1591 fAsimov->isExpected = true;
1592 }
1593
1594 return fAsimov;
1595}
1596
1598{
1599 if (fPllType != xRooFit::Asymptotics::Uncapped && ts_asymp(nSigma).first == 0)
1600 return std::pair<double, double>(1, 0);
1601 auto first_poi = dynamic_cast<RooRealVar *>(poi().first());
1602 if (!first_poi)
1603 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1604 auto _sigma_mu = sigma_mu();
1605 double nom = xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first, fNullVal(), fNullVal(), _sigma_mu.first,
1606 first_poi->getMin("physical"), first_poi->getMax("physical"));
1607 double up =
1608 xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first + ts_asymp(nSigma).second, fNullVal(), fNullVal(),
1609 _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical"));
1610 double down =
1611 xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first - ts_asymp(nSigma).second, fNullVal(), fNullVal(),
1612 _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical"));
1613 return std::pair(nom, std::max(std::abs(up - nom), std::abs(down - nom)));
1614}
1615
1617{
1618 if (fPllType != xRooFit::Asymptotics::Uncapped && ts_asymp(nSigma).first == 0)
1619 return std::pair<double, double>(1, 0);
1620 auto first_poi = dynamic_cast<RooRealVar *>(poi().first());
1621 if (!first_poi)
1622 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1623 auto _sigma_mu = sigma_mu();
1624 double nom = xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first, fNullVal(), fAltVal(), _sigma_mu.first,
1625 first_poi->getMin("physical"), first_poi->getMax("physical"));
1626 double up =
1627 xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first + ts_asymp(nSigma).second, fNullVal(), fAltVal(),
1628 _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical"));
1629 double down =
1630 xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first - ts_asymp(nSigma).second, fNullVal(), fAltVal(),
1631 _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical"));
1632
1633 return std::pair(nom, std::max(std::abs(up - nom), std::abs(down - nom)));
1634}
1635
1637{
1638 if (fNullVal() == fAltVal())
1639 return std::pair<double, double>(1, 0); // by construction
1640
1641 if (fPllType != xRooFit::Asymptotics::Uncapped && ts_asymp(nSigma).first == 0)
1642 return std::pair<double, double>(1, 0);
1643 auto first_poi = dynamic_cast<RooRealVar *>(poi().first());
1644 if (!first_poi)
1645 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1646
1647 auto _ts_asymp = ts_asymp(nSigma);
1648 auto _sigma_mu = sigma_mu();
1649 double nom1 = xRooFit::Asymptotics::PValue(fPllType, _ts_asymp.first, fNullVal(), fNullVal(), _sigma_mu.first,
1650 first_poi->getMin("physical"), first_poi->getMax("physical"));
1651 double up1 =
1652 xRooFit::Asymptotics::PValue(fPllType, _ts_asymp.first + _ts_asymp.second, fNullVal(), fNullVal(),
1653 _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical"));
1654 double down1 =
1655 xRooFit::Asymptotics::PValue(fPllType, _ts_asymp.first - _ts_asymp.second, fNullVal(), fNullVal(),
1656 _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical"));
1657 double nom2 = xRooFit::Asymptotics::PValue(fPllType, _ts_asymp.first, fNullVal(), fAltVal(), _sigma_mu.first,
1658 first_poi->getMin("physical"), first_poi->getMax("physical"));
1659 double up2 =
1660 xRooFit::Asymptotics::PValue(fPllType, _ts_asymp.first + _ts_asymp.second, fNullVal(), fAltVal(), _sigma_mu.first,
1661 first_poi->getMin("physical"), first_poi->getMax("physical"));
1662 double down2 =
1663 xRooFit::Asymptotics::PValue(fPllType, _ts_asymp.first - _ts_asymp.second, fNullVal(), fAltVal(), _sigma_mu.first,
1664 first_poi->getMin("physical"), first_poi->getMax("physical"));
1665
1666 auto nom = (nom1 == 0) ? 0 : nom1 / nom2;
1667 auto up = (up1 == 0) ? 0 : up1 / up2;
1668 auto down = (down1 == 0) ? 0 : down1 / down2;
1669
1670 return std::pair(nom, std::max(std::abs(up - nom), std::abs(down - nom)));
1671}
1672
1674{
1675 if (std::isnan(nSigma))
1676 return pll();
1677 auto first_poi = dynamic_cast<RooRealVar *>(poi().first());
1678 auto _sigma_mu = sigma_mu();
1679 if (!first_poi || (!std::isnan(nSigma) && std::isnan(_sigma_mu.first)))
1680 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1681 double nom = xRooFit::Asymptotics::k(fPllType, ROOT::Math::gaussian_cdf(nSigma), fNullVal(), fAltVal(),
1682 _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical"));
1683 double up = xRooFit::Asymptotics::k(fPllType, ROOT::Math::gaussian_cdf(nSigma), fNullVal(), fAltVal(),
1684 _sigma_mu.first + _sigma_mu.second, first_poi->getMin("physical"),
1685 first_poi->getMax("physical"));
1686 double down = xRooFit::Asymptotics::k(fPllType, ROOT::Math::gaussian_cdf(nSigma), fNullVal(), fAltVal(),
1687 _sigma_mu.first - _sigma_mu.second, first_poi->getMin("physical"),
1688 first_poi->getMax("physical"));
1689 return std::pair<double, double>(nom, std::max(std::abs(nom - up), std::abs(nom - down)));
1690}
1691
1693{
1694 if (std::isnan(nSigma))
1695 return pll();
1696 // nans should appear in the alt toys first ... so loop until past nans
1697 size_t firstToy = 0;
1698 while (firstToy < altToys.size() && std::isnan(std::get<1>(altToys[firstToy])))
1699 firstToy++;
1700 if (firstToy >= altToys.size())
1701 return std::pair(std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN());
1702 int targetIdx =
1703 (altToys.size() - firstToy) * ROOT::Math::gaussian_cdf(nSigma) + firstToy; // TODO: Account for weights
1704 return std::pair(std::get<1>(altToys[targetIdx]), (std::get<1>(altToys[std::min(int(altToys.size()), targetIdx)]) -
1705 std::get<1>(altToys[std::max(0, targetIdx)])) /
1706 2.);
1707}
1708
1710{
1711 auto _ufit = ufit(readOnly);
1712 if (!_ufit) {
1713 if (hypoTestResult)
1714 return std::pair<double, double>(hypoTestResult->GetTestStatisticData(), 0);
1715 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1716 }
1717 if (allowedStatusCodes.find(_ufit->status()) == allowedStatusCodes.end()) {
1718 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1719 }
1720 if (auto _first_poi = dynamic_cast<RooRealVar *>(poi().first());
1721 _first_poi && _first_poi->getMin("physical") > _first_poi->getMin() &&
1722 mu_hat().getVal() < _first_poi->getMin("physical")) {
1723 // replace _ufit with fit "boundary" conditional fit
1724 _ufit = cfit_lbound(readOnly);
1725 if (!_ufit) {
1726 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1727 }
1728 }
1729 auto cFactor = (fPllType == xRooFit::Asymptotics::TwoSided)
1730 ? 1.
1731 : xRooFit::Asymptotics::CompatFactor(fPllType, fNullVal(), mu_hat().getVal());
1732 if (cFactor == 0)
1733 return std::pair<double, double>(0, 0);
1734 if (!cfit_null(readOnly) || allowedStatusCodes.find(cfit_null(readOnly)->status()) == allowedStatusCodes.end())
1735 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
1736 // std::cout << cfit->minNll() << ":" << cfit->edm() << " " << ufit->minNll() << ":" << ufit->edm() << std::endl;
1737 return std::pair<double, double>(2. * cFactor * (cfit_null(readOnly)->minNll() - _ufit->minNll()),
1738 2. * cFactor * sqrt(pow(cfit_null(readOnly)->edm(), 2) + pow(_ufit->edm(), 2)));
1739 // return 2.*cFactor*(cfit->minNll()+cfit->edm() - ufit->minNll()+ufit->edm());
1740}
1741
1742std::shared_ptr<const RooFitResult> xRooNLLVar::xRooHypoPoint::retrieveFit(int type)
1743{
1744 if (!hypoTestResult)
1745 return nullptr;
1746 // see if can retrieve from that ....
1747 if (auto fits = hypoTestResult->GetFitInfo()) {
1748 for (int i = 0; i < fits->numEntries(); i++) {
1749 auto fit = fits->get(i);
1750 if (fit->getCatIndex("type") != type)
1751 continue;
1752 // found ufit ... construct
1753 std::string _name =
1754 fits->getGlobalObservables()->getStringValue(TString::Format("%s.name", fit->getCatLabel("type")));
1755 // see if can retrieve from any open file ....
1756 TDirectory *tmp = gDirectory;
1757 for (auto file : *gROOT->GetListOfFiles()) {
1758 if (auto k = static_cast<TDirectory *>(file)->FindKeyAny(_name.c_str())) {
1759 // use pre-retrieved fits if available
1760 xRooFit::StoredFitResult *storedFr =
1761 k->GetMotherDir()->GetList()
1762 ? dynamic_cast<xRooFit::StoredFitResult *>(k->GetMotherDir()->GetList()->FindObject(k->GetName()))
1763 : nullptr;
1764 if (auto cachedFit = (storedFr) ? storedFr->fr.get() : k->ReadObject<RooFitResult>(); cachedFit) {
1765 if (!storedFr) {
1766 storedFr = new xRooFit::StoredFitResult(cachedFit);
1767 k->GetMotherDir()->Add(storedFr);
1768 }
1769 gDirectory = tmp; // one of the above calls moves to key's directory ... i didn't check which
1770 return storedFr->fr;
1771 }
1772 }
1773 }
1774 auto rfit = std::make_shared<RooFitResult>(_name.c_str(), TUUID(_name.c_str()).GetTime().AsString());
1775 rfit->setStatus(fit->getRealValue("status"));
1776 rfit->setMinNLL(fit->getRealValue("minNll"));
1777 rfit->setEDM(fit->getRealValue("edm"));
1778 rfit->setCovQual(fit->getRealValue("covQual"));
1779 if (type == 0) {
1780 std::unique_ptr<RooAbsCollection> par_hats(
1781 hypoTestResult->GetFitInfo()->getGlobalObservables()->selectByName(coords->contentsString().c_str()));
1782 par_hats->setName("floatParsFinal");
1783 rfit->setFinalParList(*par_hats);
1784 } else {
1785 rfit->setFinalParList(RooArgList());
1786 }
1787 rfit->setConstParList(RooArgList());
1788 rfit->setInitParList(RooArgList());
1789 TMatrixDSym cov(0);
1790 rfit->setCovarianceMatrix(cov);
1791 return rfit;
1792 }
1793 }
1794 return nullptr;
1795}
1796
1797std::shared_ptr<const RooFitResult> xRooNLLVar::xRooHypoPoint::ufit(bool readOnly)
1798{
1799 if (fUfit)
1800 return fUfit;
1801 if (auto rfit = retrieveFit(0)) {
1802 return fUfit = rfit;
1803 }
1804 if (!nllVar || (readOnly && nllVar->get() && !nllVar->get()->getAttribute("readOnly")))
1805 return nullptr;
1806 if (!nllVar->fFuncVars)
1807 nllVar->reinitialize();
1808 AutoRestorer snap(*nllVar->fFuncVars, nllVar.get());
1809 if (!fData.first) {
1810 if (!readOnly && isExpected && fGenFit) {
1811 // can try to do a readOnly in case can load from cache
1812 bool tmp = nllVar->get()->getAttribute("readOnly");
1813 nllVar->get()->setAttribute("readOnly");
1814 auto out = ufit(true);
1815 nllVar->get()->setAttribute("readOnly", tmp);
1816 if (out) {
1817 // retrieve from cache worked, no need to generate dataset
1818 return out;
1819 } else if (!tmp) { // don't need to setData if doing a readOnly fit
1820 nllVar->setData(data());
1821 }
1822 }
1823 } else if (!nllVar->get()->getAttribute("readOnly")) { // don't need to setData if doing a readOnly fit
1824 nllVar->setData(fData);
1825 }
1826 nllVar->fFuncVars->setAttribAll("Constant", false);
1827 *nllVar->fFuncVars = *coords; // will reconst the coords
1828 if (nllVar->fFuncGlobs)
1829 nllVar->fFuncGlobs->setAttribAll("Constant", true);
1830 std::unique_ptr<RooAbsCollection>(nllVar->fFuncVars->selectCommon(poi()))
1831 ->setAttribAll("Constant", false); // float the poi
1832 if (fGenFit) {
1833 // make initial guess same as pars we generated with
1834 nllVar->fFuncVars->assignValueOnly(fGenFit->constPars());
1835 nllVar->fFuncVars->assignValueOnly(fGenFit->floatParsFinal());
1836 // rename nll so if caching fit results will cache into subdir
1837 nllVar->get()->SetName(
1838 TString::Format("%s/%s_%s", nllVar->get()->GetName(), fGenFit->GetName(), (isExpected) ? "asimov" : "toys"));
1839 if (!isExpected)
1840 nllVar->get()->SetName(TString::Format("%s/%s", nllVar->get()->GetName(), fData.first->GetName()));
1841
1842 } else if (!std::isnan(fAltVal())) {
1843 // guess data given is expected to align with alt value, unless initVal attribute specified
1844 for (auto _poiCoord : poi()) {
1845 auto _poi = dynamic_cast<RooRealVar *>(nllVar->fFuncVars->find(_poiCoord->GetName()));
1846 if (_poi) {
1847 _poi->setVal(_poi->getStringAttribute("initVal") ? TString(_poi->getStringAttribute("initVal")).Atof()
1848 : fAltVal());
1849 }
1850 }
1851 }
1852 return (fUfit = nllVar->minimize());
1853}
1854
1855std::string collectionContents(const RooAbsCollection &coll)
1856{
1857 std::string out;
1858 for (auto &c : coll) {
1859 if (!out.empty())
1860 out += ",";
1861 out += c->GetName();
1862 if (auto v = dynamic_cast<RooAbsReal *>(c); v) {
1863 out += TString::Format("=%g", v->getVal());
1864 } else if (auto cc = dynamic_cast<RooAbsCategory *>(c); cc) {
1865 out += TString::Format("=%s", cc->getLabel());
1866 } else if (auto s = dynamic_cast<RooStringVar *>(c); v) {
1867 out += TString::Format("=%s", s->getVal());
1868 }
1869 }
1870 return out;
1871}
1872
1873std::shared_ptr<const RooFitResult> xRooNLLVar::xRooHypoPoint::cfit_null(bool readOnly)
1874{
1875 if (fNull_cfit)
1876 return fNull_cfit;
1877 if (auto rfit = retrieveFit(1)) {
1878 return fNull_cfit = rfit;
1879 }
1880 if (!nllVar || (readOnly && nllVar->get() && !nllVar->get()->getAttribute("readOnly")))
1881 return nullptr;
1882 if (!nllVar->fFuncVars)
1883 nllVar->reinitialize();
1884 AutoRestorer snap(*nllVar->fFuncVars, nllVar.get());
1885 if (!fData.first) {
1886 if (!readOnly && isExpected && fGenFit) {
1887 // can try to do a readOnly in case can load from cache
1888 bool tmp = nllVar->get()->getAttribute("readOnly");
1889 nllVar->get()->setAttribute("readOnly");
1890 auto out = cfit_null(true);
1891 nllVar->get()->setAttribute("readOnly", tmp);
1892 if (out) {
1893 // retrieve from cache worked, no need to generate dataset
1894 return out;
1895 } else if (!tmp) { // don't need to setData if doing a readOnly fit
1896 nllVar->setData(data());
1897 }
1898 }
1899 } else if (!nllVar->get()->getAttribute("readOnly")) { // don't need to setData if doing a readOnly fit
1900 nllVar->setData(fData);
1901 }
1902 if (fUfit) {
1903 // move to ufit coords before evaluating
1904 *nllVar->fFuncVars = fUfit->floatParsFinal();
1905 }
1906 nllVar->fFuncVars->setAttribAll("Constant", false);
1907 *nllVar->fFuncVars = *coords; // will reconst the coords
1908 if (nllVar->fFuncGlobs)
1909 nllVar->fFuncGlobs->setAttribAll("Constant", true);
1910 if (fPOIName()) {
1911 nllVar->fFuncVars->find(fPOIName())
1912 ->setStringAttribute("altVal", (!std::isnan(fAltVal())) ? TString::Format("%g", fAltVal()) : nullptr);
1913 }
1914 if (fGenFit) {
1915 nllVar->get()->SetName(
1916 TString::Format("%s/%s_%s", nllVar->get()->GetName(), fGenFit->GetName(), (isExpected) ? "asimov" : "toys"));
1917 if (!isExpected)
1918 nllVar->get()->SetName(TString::Format("%s/%s", nllVar->get()->GetName(), fData.first->GetName()));
1919 }
1920 nllVar->get()->setStringAttribute("fitresultTitle", collectionContents(poi()).c_str());
1921 return (fNull_cfit = nllVar->minimize());
1922}
1923
1924std::shared_ptr<const RooFitResult> xRooNLLVar::xRooHypoPoint::cfit_lbound(bool readOnly)
1925{
1926 auto _first_poi = dynamic_cast<RooRealVar *>(poi().first());
1927 if (!_first_poi)
1928 return nullptr;
1929 if (_first_poi->getMin("physical") <= _first_poi->getMin())
1930 return nullptr;
1931 if (fLbound_cfit)
1932 return fLbound_cfit;
1933 if (auto rfit = retrieveFit(6)) {
1934 return fLbound_cfit = rfit;
1935 }
1936 if (!nllVar || (readOnly && nllVar->get() && !nllVar->get()->getAttribute("readOnly")))
1937 return nullptr;
1938 if (!nllVar->fFuncVars)
1939 nllVar->reinitialize();
1940 AutoRestorer snap(*nllVar->fFuncVars, nllVar.get());
1941 if (!fData.first) {
1942 if (!readOnly && isExpected && fGenFit) {
1943 // can try to do a readOnly in case can load from cache
1944 bool tmp = nllVar->get()->getAttribute("readOnly");
1945 nllVar->get()->setAttribute("readOnly");
1946 auto out = cfit_lbound(true);
1947 nllVar->get()->setAttribute("readOnly", tmp);
1948 if (out) {
1949 // retrieve from cache worked, no need to generate dataset
1950 return out;
1951 } else if (!tmp) { // don't need to setData if doing a readOnly fit
1952 nllVar->setData(data());
1953 }
1954 }
1955 } else if (!nllVar->get()->getAttribute("readOnly")) { // don't need to setData if doing a readOnly fit
1956 nllVar->setData(fData);
1957 }
1958 if (fUfit) {
1959 // move to ufit coords before evaluating
1960 *nllVar->fFuncVars = fUfit->floatParsFinal();
1961 }
1962 nllVar->fFuncVars->setAttribAll("Constant", false);
1963 *nllVar->fFuncVars = *coords; // will reconst the coords
1964 nllVar->fFuncVars->setRealValue(_first_poi->GetName(), _first_poi->getMin("physical"));
1965 if (nllVar->fFuncGlobs)
1966 nllVar->fFuncGlobs->setAttribAll("Constant", true);
1967 if (fPOIName()) {
1968 nllVar->fFuncVars->find(fPOIName())
1969 ->setStringAttribute("altVal", (!std::isnan(fAltVal())) ? TString::Format("%g", fAltVal()) : nullptr);
1970 }
1971 if (fGenFit) {
1972 nllVar->get()->SetName(
1973 TString::Format("%s/%s_%s", nllVar->get()->GetName(), fGenFit->GetName(), (isExpected) ? "asimov" : "toys"));
1974 if (!isExpected)
1975 nllVar->get()->SetName(TString::Format("%s/%s", nllVar->get()->GetName(), fData.first->GetName()));
1976 }
1977 nllVar->get()->setStringAttribute(
1978 "fitresultTitle",
1979 collectionContents(*std::unique_ptr<RooAbsCollection>(nllVar->fFuncVars->selectCommon(poi()))).c_str());
1980 return (fLbound_cfit = nllVar->minimize());
1981}
1982
1983std::shared_ptr<const RooFitResult> xRooNLLVar::xRooHypoPoint::cfit_alt(bool readOnly)
1984{
1985 if (std::isnan(fAltVal()))
1986 return nullptr;
1987 if (fAlt_cfit)
1988 return fAlt_cfit;
1989 if (auto rfit = retrieveFit(2)) {
1990 return fAlt_cfit = rfit;
1991 }
1992 if (!nllVar || (readOnly && nllVar->get() && !nllVar->get()->getAttribute("readOnly")))
1993 return nullptr;
1994 if (!nllVar->fFuncVars)
1995 nllVar->reinitialize();
1996 AutoRestorer snap(*nllVar->fFuncVars, nllVar.get());
1997 if (!fData.first) {
1998 if (!readOnly && isExpected && fGenFit) {
1999 // can try to do a readOnly in case can load from cache
2000 bool tmp = nllVar->get()->getAttribute("readOnly");
2001 nllVar->get()->setAttribute("readOnly");
2002 auto out = cfit_alt(true);
2003 nllVar->get()->setAttribute("readOnly", tmp);
2004 if (out) {
2005 // retrieve from cache worked, no need to generate dataset
2006 return out;
2007 } else if (!tmp) { // don't need to setData if doing a readOnly fit
2008 nllVar->setData(data());
2009 }
2010 }
2011 } else if (!nllVar->get()->getAttribute("readOnly")) { // don't need to setData if doing a readOnly fit
2012 nllVar->setData(fData);
2013 }
2014 if (fUfit) {
2015 // move to ufit coords before evaluating
2016 *nllVar->fFuncVars = fUfit->floatParsFinal();
2017 }
2018 nllVar->fFuncVars->setAttribAll("Constant", false);
2019 *nllVar->fFuncVars = *coords; // will reconst the coords
2020 if (nllVar->fFuncGlobs)
2021 nllVar->fFuncGlobs->setAttribAll("Constant", true);
2022 *nllVar->fFuncVars = alt_poi();
2023 if (fGenFit) {
2024 nllVar->get()->SetName(
2025 TString::Format("%s/%s_%s", nllVar->get()->GetName(), fGenFit->GetName(), (isExpected) ? "asimov" : "toys"));
2026 if (!isExpected)
2027 nllVar->get()->SetName(TString::Format("%s/%s", nllVar->get()->GetName(), fData.first->GetName()));
2028 }
2029 nllVar->get()->setStringAttribute("fitresultTitle", collectionContents(alt_poi()).c_str());
2030 return (fAlt_cfit = nllVar->minimize());
2031}
2032
2034{
2035
2036 auto asi = asimov(readOnly);
2037
2038 if (!asi) {
2039 return std::pair<double, double>(std::numeric_limits<double>::quiet_NaN(), 0);
2040 }
2041
2042 auto out = asi->pll(readOnly);
2043 return std::pair<double, double>(std::abs(fNullVal() - fAltVal()) / sqrt(out.first),
2044 out.second * 0.5 * std::abs(fNullVal() - fAltVal()) /
2045 (out.first * sqrt(out.first)));
2046}
2047
2049{
2050 auto _ts = ts_toys(nSigma);
2051 if (std::isnan(_ts.first))
2052 return _ts;
2053 if (fPllType != xRooFit::Asymptotics::Uncapped && _ts.first == 0)
2054 return std::pair<double, double>(1, 0); // don't need toys to compute this point!
2055
2056 TEfficiency eff("", "", 1, 0, 1);
2057
2058 auto &_theToys = (alt) ? altToys : nullToys;
2059
2060 if (_theToys.empty()) {
2061 return std::pair(0.5, std::numeric_limits<double>::infinity());
2062 }
2063
2064 // loop over toys, count how many are > ts value
2065 // nans (mean bad ts evaluations) will count towards uncertainty
2066 int nans = 0;
2067 double result = 0;
2068 double result_err_up = 0;
2069 double result_err_down = 0;
2070 for (auto &toy : _theToys) {
2071 if (std::isnan(std::get<1>(toy))) {
2072 nans++;
2073 } else {
2074 bool res = std::get<1>(toy) >= _ts.first;
2075 if (std::get<2>(toy) != 1) {
2076 eff.FillWeighted(res, 0.5, std::get<2>(toy));
2077 } else {
2078 eff.Fill(res, 0.5);
2079 }
2080 if (res)
2081 result += std::get<2>(toy);
2082 if (std::get<1>(toy) >= _ts.first - _ts.second)
2083 result_err_up += std::get<2>(toy);
2084 if (std::get<1>(toy) >= _ts.first - _ts.second)
2085 result_err_down += std::get<2>(toy);
2086 }
2087 }
2088 // symmetrize the error
2089 result_err_up -= result;
2090 result_err_down -= result;
2091 double result_err = std::max(std::abs(result_err_up), std::abs(result_err_down));
2092 // assume the nans would "add" to the p-value, conservative scenario
2093 result_err += nans;
2094 result_err /= _theToys.size();
2095
2096 // don't include the nans for the central value though
2097 result /= (_theToys.size() - nans);
2098
2099 // add to the result_err (in quadrature) the uncert due to limited stats
2100 result_err = sqrt(result_err * result_err + eff.GetEfficiencyErrorUp(1) * eff.GetEfficiencyErrorUp(1));
2101 return std::pair<double, double>(result, result_err);
2102}
2103
2105{
2106 return pX_toys(false, nSigma);
2107}
2108
2110{
2111 if (!std::isnan(nSigma)) {
2112 return std::pair<double, double>(ROOT::Math::gaussian_cdf(nSigma), 0); // by construction
2113 }
2114 return pX_toys(true, nSigma);
2115}
2116
2118{
2119 xRooHypoPoint out;
2120 out.coords = coords;
2121 out.fPllType = fPllType; // out.fPOIName = fPOIName; out.fNullVal=fNullVal; out.fAltVal = fAltVal;
2122 out.nllVar = nllVar;
2123 if (!nllVar)
2124 return out;
2125 auto _cfit = cfit_null();
2126 if (!_cfit)
2127 return out;
2128 if (!nllVar->fFuncVars)
2129 nllVar->reinitialize();
2130 //*nllVar->fFuncVars = cfit_null()->floatParsFinal();
2131 //*nllVar->fFuncVars = cfit_null()->constPars();
2132 out.fData = xRooFit::generateFrom(*nllVar->fPdf, *_cfit, false, seed); // nllVar->generate(false,seed);
2133 out.fGenFit = _cfit;
2134 return out;
2135}
2136
2138{
2139 xRooHypoPoint out;
2140 out.coords = coords;
2141 out.fPllType = fPllType; // out.fPOIName = fPOIName; out.fNullVal=fNullVal; out.fAltVal = fAltVal;
2142 out.nllVar = nllVar;
2143 if (!nllVar)
2144 return out;
2145 if (!cfit_alt())
2146 return out;
2147 if (!nllVar->fFuncVars)
2148 nllVar->reinitialize();
2149 //*nllVar->fFuncVars = cfit_alt()->floatParsFinal();
2150 //*nllVar->fFuncVars = cfit_alt()->constPars();
2151 out.fData =
2152 xRooFit::generateFrom(*nllVar->fPdf, *cfit_alt(), false, seed); // out.data = nllVar->generate(false,seed);
2153 out.fGenFit = cfit_alt();
2154 return out;
2155}
2156
2157size_t xRooNLLVar::xRooHypoPoint::addToys(bool alt, int nToys, int initialSeed, double target, double target_nSigma,
2158 bool targetCLs, double relErrThreshold, size_t maxToys)
2159{
2160 if ((alt && !cfit_alt()) || (!alt && !cfit_null())) {
2161 throw std::runtime_error("Cannot add toys, invalid conditional fit");
2162 }
2163
2164 auto condition = [&]() { // returns true if need more toys
2165 if (std::isnan(target))
2166 return false;
2167 auto obs = targetCLs ? pCLs_toys(target_nSigma) : (alt ? pAlt_toys(target_nSigma) : pNull_toys(target_nSigma));
2168 if (!std::isnan(obs.first)) {
2169 double diff = (target < 0) ? obs.first : std::abs(obs.first - target);
2170 double err = obs.second;
2171 if (err > 1e-4 && diff <= relErrThreshold * obs.second) {
2172 // return true; // more toys needed
2173 if (targetCLs) {
2174 // decide which type we'd want to generate and update alt flag
2175 auto pNull = pNull_toys(target_nSigma);
2176 auto pAlt = pAlt_toys(target_nSigma);
2177 // std::cout << obs.first << " +/- " << obs.second << ": " << pNull.first << " +/- " << pNull.second << "
2178 // , " << pAlt.first << " +/- " << pAlt.second << std::endl;
2179 alt = (pAlt.second * pNull.first > pNull.second * pAlt.first);
2180 if ((alt ? pAlt.second : pNull.second) < 1e-4)
2181 return false; // stop if error gets too small
2182 }
2183 return true;
2184 }
2185 }
2186 return false;
2187 };
2188
2189 if (!std::isnan(target) && std::isnan(ts_toys(target_nSigma).first)) {
2190 if (std::isnan(target_nSigma)) {
2191 throw std::runtime_error("Cannot target obs p-value because ts value unavailable");
2192 }
2193 if (targetCLs && pCLs_toys(target_nSigma).second == 0) {
2194 // this happens if the mu_test=mu_alt ... no toys needed
2195 return 0;
2196 }
2197
2198 // try generating 100 alt toys
2199 Info("addToys", "First generating 100 alt toys in order to determine expected ts value");
2200 addToys(true, 100, initialSeed);
2201 // if still null then exit
2202 if (std::isnan(ts_toys(target_nSigma).first)) {
2203 throw std::runtime_error("Unable to determine expected ts value");
2204 }
2205 }
2206
2207 size_t nans = 0;
2208 float lastTime = 0;
2209 int lasti = 0;
2210 auto g = gROOT->Get<TGraph>("toyTime");
2211 if (!g) {
2212 g = new TGraph;
2213 g->SetNameTitle("toyTime", "Time per toy;Toy;time [s]");
2214 gROOT->Add(g);
2215 }
2216 g->Set(0);
2217 TStopwatch s2;
2218 s2.Start();
2219 TStopwatch s;
2220 s.Start();
2221
2222 size_t toysAdded(0);
2223 size_t altToysAdded(0);
2224 if (initialSeed) {
2225 RooRandom::randomGenerator()->SetSeed(initialSeed);
2226 }
2227 do {
2228 auto &toys = (alt) ? altToys : nullToys;
2229 if (toys.size() >= maxToys) {
2230 // cannot generate more toys, reached limit already
2231 break;
2232 }
2233 // don't generate toys if reached target
2234 if (!std::isnan(target) && !condition()) {
2235 break;
2236 }
2237 auto currVal = std::isnan(target) ? std::pair(0., 0.)
2238 : (targetCLs ? pCLs_toys(target_nSigma)
2239 : (alt ? pAlt_toys(target_nSigma) : pNull_toys(target_nSigma)));
2240 size_t nnToys = std::min(size_t(nToys), (maxToys - toys.size()));
2241
2242 for (size_t i = 0; i < nnToys; i++) {
2243 int seed = RooRandom::randomGenerator()->Integer(std::numeric_limits<uint32_t>::max());
2244 auto toy = ((alt) ? generateAlt(seed) : generateNull(seed));
2245 TDirectory *tmp = gDirectory;
2246 gDirectory = nullptr; // disables any saving of fit results for toys
2247 toys.push_back(std::make_tuple(seed, toy.pll().first, 1.));
2248 gDirectory = tmp;
2249 (alt ? altToysAdded : toysAdded)++;
2250 if (std::isnan(std::get<1>(toys.back())))
2251 nans++;
2252 g->SetPoint(g->GetN(), g->GetN(), s.RealTime() - lastTime); // stops the clock
2253 lastTime = s.RealTime();
2254 if (s.RealTime() > 10) {
2255 std::cout << "\r"
2256 << TString::Format("Generated %d/%d %s hypothesis toys [%.2f toys/s]",
2257 int(alt ? altToysAdded : toysAdded), int(nnToys), alt ? "alt" : "null",
2258 double(altToysAdded + toysAdded - lasti) / s.RealTime());
2259 if (!std::isnan(target)) {
2260 std::cout << " [current=" << currVal.first << "+/-" << currVal.second << " target=" << target
2261 << " nSigma=" << target_nSigma << "]";
2262 }
2263 std::cout << "..." << std::flush;
2264 lasti = altToysAdded + toysAdded;
2265 s.Reset();
2266 Draw();
2267 if (gPad) {
2268 gPad->Update();
2270 }
2271 s.Start();
2272 // std::cout << "Generated " << i << "/" << nToys << (alt ? " alt " : " null ") << " hypothesis toys " ..."
2273 // << std::endl;
2274 }
2275 s.Continue();
2276 }
2277 // sort the toys ... put nans first - do by setting all as negative inf (is that still necessary with the custom
2278 // sort below??)
2279 for (auto &t : toys) {
2280 if (std::isnan(std::get<1>(t)))
2281 std::get<1>(t) = -std::numeric_limits<double>::infinity();
2282 }
2283 std::sort(toys.begin(), toys.end(),
2284 [](const decltype(nullToys)::value_type &a, const decltype(nullToys)::value_type &b) -> bool {
2285 if (std::isnan(std::get<1>(a)))
2286 return true;
2287 if (std::isnan(std::get<1>(b)))
2288 return false;
2289 return std::get<1>(a) < std::get<1>(b);
2290 });
2291 for (auto &t : toys) {
2292 if (std::isinf(std::get<1>(t)))
2293 std::get<1>(t) = std::numeric_limits<double>::quiet_NaN();
2294 }
2295 if (std::isnan(target)) {
2296 break; // no more toys if not doing a target
2297 }
2298 // if(condition()) {
2299 // Info("addToys","Generating more toys to determine p-value ... currently: %f +/-
2300 // %f",pNull_toys(target_nSigma).first,pNull_toys(target_nSigma).second);
2301 // }
2302 } while (condition());
2303 if (lasti) {
2304 std::cout << "\r"
2305 << "Finished Generating ";
2306 if (toysAdded) {
2307 std::cout << toysAdded << " null ";
2308 }
2309 if (altToysAdded) {
2310 std::cout << altToysAdded << " alt ";
2311 }
2312 std::cout << "toys " << TString::Format("[%.2f toys/s overall]", double(toysAdded + altToysAdded) / s2.RealTime())
2313 << std::endl;
2314 Draw();
2315 if (gPad) {
2316 gPad->Update();
2317#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00)
2318 gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook
2319#endif
2321 }
2322 }
2323
2324 if (nans > 0) {
2325 std::cout << "Warning: " << nans << " toys were bad" << std::endl;
2326 }
2327 return toysAdded;
2328}
2329
2330void xRooNLLVar::xRooHypoPoint::addNullToys(int nToys, int seed, double target, double target_nSigma)
2331{
2332 addToys(false, nToys, seed, target, target_nSigma);
2333}
2334void xRooNLLVar::xRooHypoPoint::addAltToys(int nToys, int seed, double target, double target_nSigma)
2335{
2336 addToys(true, nToys, seed, target, target_nSigma);
2337}
2338
2339void xRooNLLVar::xRooHypoPoint::addCLsToys(int nToys, int seed, double target, double target_nSigma)
2340{
2341 addToys(false, nToys, seed, target, target_nSigma, true);
2342 return;
2343 //
2344 // auto condition = [&](bool doingAlt=false) { // returns true if need more toys
2345 // if(std::isnan(target)) return false;
2346 // auto pval = pCLs_toys(target_nSigma);
2347 // if (!std::isnan(pval.first)) {
2348 // double diff = std::abs(pval.first - target);
2349 // double err = pval.second;
2350 // if (err > 1e-4 && diff <= 2 * pval.second) {
2351 // return true; // more toys needed
2352 // // decide which type we'd want to generate
2353 // // if it matches the type we are generating, then return true
2354 // auto pNull = pNull_toys(target_nSigma);
2355 // auto pAlt = pAlt_toys(target_nSigma);
2356 // if ((doingAlt ? pAlt.second : pNull.second) < 1e-4) return false; // stop if error gets too small
2357 // bool doAlt = (pAlt.second * pNull.first > pNull.second * pAlt.first);
2358 // return doAlt == doingAlt;
2359 // }
2360 // }
2361 // return false;
2362 // };
2363 // while(condition()) {
2364 // bool doAlt = false;
2365 // double relErrThreshold = 2;
2366 // if(nullToys.size()<size_t(nToys)) {
2367 // addToys(false,nToys);continue;
2368 // } else if(altToys.size()<size_t(nToys)) {
2369 // addToys(true,nToys);continue;
2370 // } else {
2371 // // see which have bigger errors ... generate more of that ...
2372 // auto pNull = pNull_toys(target_nSigma);
2373 // auto pAlt = pAlt_toys(target_nSigma);
2374 // doAlt = (pAlt.second*pNull.first > pNull.second*pAlt.first);
2375 // if( (doAlt ? pAlt.second : pNull.second) < 1e-4 ) break; // stop if error gets too small
2376 // auto pCLs = pCLs_toys(target_nSigma);
2377 // relErrThreshold = (doAlt) ? (pNull.second/pNull.first) : (pAlt.second/pAlt.first);
2378 // relErrThreshold = std::min(2.,std::abs(relErrThreshold));
2379 // std::cout << "Current pCLs = " << pCLs.first << " +/- " << pCLs.second
2380 // << " (pNull = " << pNull.first << " +/- " << pNull.second
2381 // << " , pAlt = " << pAlt.first << " +/- " << pAlt.second << ") ... generating more " << (doAlt ?
2382 // "alt" : "null") << " toys " << relErrThreshold << std::endl;
2383 //
2384 // }
2385 // if( addToys(doAlt, nToys/*, seed, -1, target_nSigma,relErrThreshold*/) == 0) {
2386 // break; // no toys got added, so stop looping
2387 // }
2388 // }
2389}
2390
2392xRooNLLVar::hypoPoint(const char *poiValues, double alt_value, const xRooFit::Asymptotics::PLLType &pllType)
2393{
2394 xRooHypoPoint out;
2395 // out.fPOIName = parName; out.fNullVal = value; out.fAltVal = alt_value;
2396
2397 if (!fFuncVars) {
2398 reinitialize();
2399 }
2400 AutoRestorer snap(*fFuncVars);
2401
2402 out.nllVar = std::make_shared<xRooNLLVar>(*this);
2403 out.fData = getData();
2404
2405 TStringToken pattern(poiValues, ",");
2406 TString poiNames;
2407 while (pattern.NextToken()) {
2408 TString s = pattern.Data();
2409 TString cName = s;
2410 double val = std::numeric_limits<double>::quiet_NaN();
2411 auto i = s.Index("=");
2412 if (i != -1) {
2413 cName = s(0, i);
2414 TString cVal = s(i + 1, s.Length());
2415 if (!cVal.IsFloat())
2416 throw std::runtime_error("poiValues must contain value");
2417 val = cVal.Atof();
2418 }
2419 auto v = dynamic_cast<RooRealVar *>(fFuncVars->find(cName));
2420 if (!v)
2421 throw std::runtime_error("Cannot find poi");
2422 if (!std::isnan(val))
2423 v->setVal(val);
2424 v->setConstant(); // because will select constants as coords
2425 if (poiNames != "") {
2426 poiNames += ",";
2427 }
2428 poiNames += cName;
2429 }
2430 if (poiNames == "") {
2431 throw std::runtime_error("No poi");
2432 }
2433 if (!std::isnan(alt_value)) {
2434 std::unique_ptr<RooAbsCollection> thePoi(fFuncVars->selectByName(poiNames));
2435 for (auto b : *thePoi) {
2436 if (!static_cast<RooRealVar *>(b)->hasRange("physical")) {
2437 static_cast<RooRealVar *>(b)->setRange("physical", 0, std::numeric_limits<double>::infinity());
2438 }
2439 }
2440 }
2441 auto _snap = std::unique_ptr<RooAbsCollection>(fFuncVars->selectByAttrib("Constant", true))->snapshot();
2442 _snap->setAttribAll("poi", false);
2443 std::unique_ptr<RooAbsCollection> _poi(_snap->selectByName(poiNames));
2444 _poi->setAttribAll("poi", true);
2445 if (std::isnan(alt_value)) {
2446 for (auto a : *_poi)
2447 a->setStringAttribute("altVal", nullptr);
2448 } else {
2449 for (auto a : *_poi)
2450 a->setStringAttribute("altVal", TString::Format("%g", alt_value));
2451 }
2452 if (fGlobs)
2453 _snap->remove(*fGlobs, true, true);
2454 out.coords.reset(_snap);
2455
2456 auto _type = pllType;
2457 if (_type == xRooFit::Asymptotics::Unknown) {
2458 // decide based on values
2459 if (std::isnan(alt_value)) {
2461 } else if (dynamic_cast<RooRealVar *>(_poi->first())->getVal() >= alt_value) {
2463 } else {
2465 }
2466 }
2467
2468 out.fPllType = _type;
2469
2470 return out;
2471}
2472
2473xRooNLLVar::xRooHypoPoint
2474xRooNLLVar::hypoPoint(double value, double alt_value, const xRooFit::Asymptotics::PLLType &pllType)
2475{
2476 if (!fFuncVars) {
2477 reinitialize();
2478 }
2479 std::unique_ptr<RooAbsCollection> _poi(fFuncVars->selectByAttrib("poi", true));
2480 if (_poi->empty()) {
2481 throw std::runtime_error("No POI specified in model");
2482 } else if (_poi->size() != 1) {
2483 throw std::runtime_error("Multiple POI specified in model");
2484 }
2485 return hypoPoint(_poi->first()->GetName(), value, alt_value, pllType);
2486}
2487
2489xRooNLLVar::hypoPoint(const char *parName, double value, double alt_value, const xRooFit::Asymptotics::PLLType &pllType)
2490{
2491 return hypoPoint(TString::Format("%s=%f", parName, value), alt_value, pllType);
2492}
2493
2495{
2496
2497 if (!nllVar && !hypoTestResult)
2498 return;
2499
2500 TString sOpt(opt);
2501 sOpt.ToLower();
2502 bool hasSame = sOpt.Contains("same");
2503 sOpt.ReplaceAll("same", "");
2504
2505 TVirtualPad *pad = gPad;
2506
2507 TH1 *hAxis = nullptr;
2508
2509 auto clearPad = []() {
2510 gPad->Clear();
2511 if (gPad->GetNumber() == 0) {
2512 gPad->SetBottomMargin(gStyle->GetPadBottomMargin());
2513 gPad->SetTopMargin(gStyle->GetPadTopMargin());
2514 gPad->SetLeftMargin(gStyle->GetPadLeftMargin());
2515 gPad->SetRightMargin(gStyle->GetPadRightMargin());
2516 }
2517 };
2518
2519 if (!hasSame || !pad) {
2520 if (!pad) {
2522 pad = gPad;
2523 }
2524 clearPad();
2525 } else {
2526 // get the histogram representing the axes
2527 hAxis = dynamic_cast<TH1 *>(pad->GetPrimitive(".axis"));
2528 if (!hAxis) {
2529 for (auto o : *pad->GetListOfPrimitives()) {
2530 if (hAxis = dynamic_cast<TH1 *>(o); hAxis)
2531 break;
2532 }
2533 }
2534 }
2535
2536 // get min and max values
2537 double _min = std::numeric_limits<double>::quiet_NaN();
2538 double _max = -std::numeric_limits<double>::quiet_NaN();
2539
2540 for (auto &p : nullToys) {
2541 if (std::get<2>(p) == 0)
2542 continue;
2543 if (std::isnan(std::get<1>(p)))
2544 continue;
2545 _min = std::min(std::get<1>(p), _min);
2546 _max = std::max(std::get<1>(p), _max);
2547 }
2548 for (auto &p : altToys) {
2549 if (std::get<2>(p) == 0)
2550 continue;
2551 if (std::isnan(std::get<1>(p)))
2552 continue;
2553 _min = std::min(std::get<1>(p), _min);
2554 _max = std::max(std::get<1>(p), _max);
2555 }
2556
2557 auto obs = pll();
2558 if (!std::isnan(obs.first)) {
2559 _min = std::min(obs.first - std::abs(obs.first) * 0.1, _min);
2560 _max = std::max(obs.first + std::abs(obs.first) * 0.1, _max);
2561 }
2562 // these are used down below to add obs p-values to legend, but up here because can trigger fits that create asimov
2563 auto pNull = pNull_toys();
2564 auto pAlt = pAlt_toys();
2565 auto pNullA = pNull_asymp();
2566 auto pAltA = pAlt_asymp();
2567 sigma_mu(true);
2568 auto asi = (fAsimov && fAsimov->fUfit && fAsimov->fNull_cfit) ? fAsimov->pll().first
2569 : std::numeric_limits<double>::quiet_NaN();
2570 if (!std::isnan(asi) && asi > 0 && fPllType != xRooFit::Asymptotics::Unknown) {
2571 // can calculate asymptotic distributions,
2572 _min = std::min(asi - std::abs(asi), _min);
2573 _max = std::max(asi + std::abs(asi), _max);
2574 }
2575 if (_min > 0)
2576 _min = 0;
2577
2578 auto _poi = dynamic_cast<RooRealVar *>(poi().first());
2579
2580 auto makeHist = [&](bool isAlt) {
2581 TString title;
2582 auto h = new TH1D((isAlt) ? "alt_toys" : "null_toys", "", 100, _min, _max + (_max - _min) * 0.01);
2583 h->SetDirectory(nullptr);
2584 size_t nBadOrZero = 0;
2585 for (auto &p : (isAlt) ? altToys : nullToys) {
2586 double w = std::isnan(std::get<1>(p)) ? 0 : std::get<2>(p);
2587 if (w == 0)
2588 nBadOrZero++;
2589 if (!std::isnan(std::get<1>(p)))
2590 h->Fill(std::get<1>(p), w);
2591 }
2592 if (h->GetEntries() > 0)
2593 h->Scale(1. / h->Integral(0, h->GetNbinsX() + 1));
2594
2595 // add POI values to identify hypos
2596 // for(auto p : *fPOI) {
2597 // if (auto v = dynamic_cast<RooRealVar*>(p)) {
2598 // if (auto v2 = dynamic_cast<RooRealVar*>(fAltPoint->fCoords->find(*v)); v2 &&
2599 // v2->getVal()!=v->getVal()) {
2600 // // found point that differs in poi and altpoint value, so print my coords value for this
2601 // title += TString::Format("%s' = %g,
2602 // ",v->GetTitle(),dynamic_cast<RooRealVar*>(fCoords->find(*v))->getVal());
2603 // }
2604 // }
2605 // }
2606 if (fPOIName())
2607 title += TString::Format("%s' = %g", fPOIName(), (isAlt) ? fAltVal() : fNullVal());
2608 title += TString::Format(" , N_{toys}=%d", int((isAlt) ? altToys.size() : nullToys.size()));
2609 if (nBadOrZero > 0)
2610 title += TString::Format(" (N_{bad/0}=%d)", int(nBadOrZero));
2611 title += ";";
2612 title += tsTitle();
2613 title += TString::Format(";Probability Mass");
2614 h->SetTitle(title);
2615 h->SetLineColor(isAlt ? kRed : kBlue);
2616 h->SetLineWidth(2);
2617 h->SetMarkerSize(0);
2618 h->SetBit(kCanDelete);
2619 return h;
2620 };
2621
2622 auto nullHist = makeHist(false);
2623 auto altHist = makeHist(true);
2624
2625 TLegend *l = nullptr;
2626 auto h = (nullHist->GetEntries()) ? nullHist : altHist;
2627 if (!hasSame) {
2628 gPad->SetLogy();
2629 auto axis = static_cast<TH1 *>(h->Clone(".axis"));
2630 axis->SetBit(kCanDelete);
2631 axis->SetStats(false);
2632 axis->Reset("ICES");
2633 axis->SetTitle(TString::Format("%s HypoPoint", collectionContents(poi()).c_str()));
2634 axis->SetLineWidth(0);
2635 axis->Draw(""); // h->Draw("axis"); cant use axis option if want title drawn
2636 axis->SetMinimum(1e-7);
2637 axis->GetYaxis()->SetRangeUser(1e-7, 10);
2638 axis->SetMaximum(h->GetMaximum());
2639 hAxis = axis;
2640 l = new TLegend(0.4, 0.7, 1. - gPad->GetRightMargin(), 1. - gPad->GetTopMargin());
2641 l->SetName("legend");
2642 l->SetFillStyle(0);
2643 l->SetBorderSize(0);
2645 l->Draw();
2646 } else {
2647 for (auto o : *gPad->GetListOfPrimitives()) {
2648 l = dynamic_cast<TLegend *>(o);
2649 if (l)
2650 break;
2651 }
2652 }
2653
2654 if (h->GetEntries() > 0) {
2655 h->Draw("esame");
2656 } else {
2657 h->Draw("axissame"); // for unknown reason if second histogram empty it still draws with two weird bars???
2658 }
2659 h = altHist;
2660 if (h->GetEntries() > 0) {
2661 h->Draw("esame");
2662 } else {
2663 h->Draw("axissame"); // for unknown reason if second histogram empty it still draws with two weird bars???
2664 }
2665
2666 if (l) {
2667 l->AddEntry(nullHist);
2668 l->AddEntry(altHist);
2669 }
2670
2671 if (fAsimov && fAsimov->fUfit && fAsimov->fNull_cfit && !std::isnan(sigma_mu().first) && !std::isnan(fAltVal())) {
2672 auto hh = static_cast<TH1 *>(nullHist->Clone("null_asymp"));
2673 hh->SetBit(kCanDelete);
2674 hh->SetStats(false);
2675 hh->SetLineStyle(2);
2676 hh->Reset();
2677 for (int i = 1; i <= hh->GetNbinsX(); i++) {
2678 hh->SetBinContent(
2679 i, xRooFit::Asymptotics::PValue(fPllType, hh->GetBinLowEdge(i), fNullVal(), fNullVal(), sigma_mu().first,
2680 _poi->getMin("physical"), _poi->getMax("physical")) -
2681 xRooFit::Asymptotics::PValue(fPllType, hh->GetBinLowEdge(i + 1), fNullVal(), fNullVal(),
2682 sigma_mu().first, _poi->getMin("physical"), _poi->getMax("physical")));
2683 }
2684 hh->Draw("lsame");
2685 hh = static_cast<TH1 *>(altHist->Clone("alt_asymp"));
2686 hh->SetBit(kCanDelete);
2687 hh->SetStats(false);
2688 hh->SetLineStyle(2);
2689 hh->Reset();
2690 for (int i = 1; i <= hh->GetNbinsX(); i++) {
2691 hh->SetBinContent(
2692 i, xRooFit::Asymptotics::PValue(fPllType, hh->GetBinLowEdge(i), fNullVal(), fAltVal(), sigma_mu().first,
2693 _poi->getMin("physical"), _poi->getMax("physical")) -
2694 xRooFit::Asymptotics::PValue(fPllType, hh->GetBinLowEdge(i + 1), fNullVal(), fAltVal(),
2695 sigma_mu().first, _poi->getMin("physical"), _poi->getMax("physical")));
2696 }
2697 hh->Draw("lsame");
2698 }
2699
2700 // draw observed points
2701 TLine ll;
2702 ll.SetLineStyle(1);
2703 ll.SetLineWidth(3);
2704 // for(auto p : fObs) {
2705 auto tl = ll.DrawLine(obs.first, hAxis->GetMinimum(), obs.first, 0.1);
2706 auto label = TString::Format("obs ts = %.4f", obs.first);
2707 if (obs.second)
2708 label += TString::Format(" #pm %.4f", obs.second);
2709
2710 l->AddEntry(tl, label, "l");
2711 label = "";
2712 if (!std::isnan(pNull.first) || !std::isnan(pAlt.first)) {
2713 auto pCLs = pCLs_toys();
2714 label += " p_{toy}=(";
2715 label += (std::isnan(pNull.first)) ? "-" : TString::Format("%.4f #pm %.4f", pNull.first, pNull.second);
2716 label += (std::isnan(pAlt.first)) ? ",-" : TString::Format(",%.4f #pm %.4f", pAlt.first, pAlt.second);
2717 label += (std::isnan(pCLs.first)) ? ",-)" : TString::Format(",%.4f #pm %.4f)", pCLs.first, pCLs.second);
2718 }
2719 if (label.Length() > 0)
2720 l->AddEntry("", label, "");
2721 label = "";
2722 if (!std::isnan(pNullA.first) || !std::isnan(pAltA.first)) {
2723 auto pCLs = pCLs_asymp();
2724 label += " p_{asymp}=(";
2725 label += (std::isnan(pNullA.first)) ? "-" : TString::Format("%.4f #pm %.4f", pNullA.first, pNullA.second);
2726 label += (std::isnan(pAltA.first)) ? ",-" : TString::Format(",%.4f #pm %.4f", pAltA.first, pAltA.second);
2727 label += (std::isnan(pCLs.first)) ? ",-)" : TString::Format(",%.4f #pm %.4f)", pCLs.first, pCLs.second);
2728 }
2729 if (label.Length() > 0)
2730 l->AddEntry("", label, "");
2731
2732 if (auto ax = dynamic_cast<TH1 *>(gPad->GetPrimitive(".axis")))
2733 ax->GetYaxis()->SetRangeUser(1e-7, 1);
2734}
2735
2737{
2738 auto v = dynamic_cast<RooRealVar *>(poi().empty() ? nullptr : poi().first());
2740 if (v && v->hasRange("physical") && v->getMin("physical") != -std::numeric_limits<double>::infinity()) {
2741 return (inWords) ? TString::Format("Lower-Bound One-Sided Limit PLR")
2742 : TString::Format("#tilde{q}_{%s=%g}", v->GetTitle(), v->getVal());
2743 } else if (v) {
2744 return (inWords) ? TString::Format("One-Sided Limit PLR")
2745 : TString::Format("q_{%s=%g}", v->GetTitle(), v->getVal());
2746 } else {
2747 return "q";
2748 }
2749 } else if (fPllType == xRooFit::Asymptotics::TwoSided) {
2750 if (v && v->hasRange("physical") && v->getMin("physical") != -std::numeric_limits<double>::infinity()) {
2751 return (inWords) ? TString::Format("Lower-Bound PLR")
2752 : TString::Format("#tilde{t}_{%s=%g}", v->GetTitle(), v->getVal());
2753 } else if (v) {
2754 return (inWords) ? TString::Format("-2log[L(%s,#hat{#hat{#theta}})/L(#hat{%s},#hat{#theta})]", v->GetTitle(),
2755 v->GetTitle())
2756 : TString::Format("t_{%s=%g}", v->GetTitle(), v->getVal());
2757 } else
2758 return "t";
2759 } else if (fPllType == xRooFit::Asymptotics::OneSidedNegative) {
2760 if (v && v->hasRange("physical") && v->getMin("physical") != -std::numeric_limits<double>::infinity()) {
2761 return (inWords) ? TString::Format("Lower-Bound One-Sided Discovery PLR")
2762 : TString::Format("#tilde{r}_{%s=%g}", v->GetTitle(), v->getVal());
2763 } else if (v) {
2764 return (inWords) ? TString::Format("One-Sided Discovery PLR")
2765 : TString::Format("r_{%s=%g}", v->GetTitle(), v->getVal());
2766 } else {
2767 return "r";
2768 }
2769 } else if (fPllType == xRooFit::Asymptotics::Uncapped) {
2770 if (v && v->hasRange("physical") && v->getMin("physical") != -std::numeric_limits<double>::infinity()) {
2771 return (inWords) ? TString::Format("Lower-Bound Uncapped PLR")
2772 : TString::Format("#tilde{s}_{%s=%g}", v->GetTitle(), v->getVal());
2773 } else if (v) {
2774 return (inWords) ? TString::Format("Uncapped PLR") : TString::Format("s_{%s=%g}", v->GetTitle(), v->getVal());
2775 } else {
2776 return "s";
2777 }
2778 } else {
2779 return "Test Statistic";
2780 }
2781}
2782
2784{
2785 return (poi().empty()) ? nullptr : (poi().first())->GetName();
2786}
2788{
2789 auto first_poi = dynamic_cast<RooAbsReal *>(poi().first());
2790 return (first_poi == nullptr) ? std::numeric_limits<double>::quiet_NaN() : first_poi->getVal();
2791}
2793{
2794 auto _alt_poi = alt_poi(); // need to keep alive as alt_poi owns its contents
2795 auto first_poi = dynamic_cast<RooAbsReal *>(_alt_poi.first());
2796 return (first_poi == nullptr) ? std::numeric_limits<double>::quiet_NaN() : first_poi->getVal();
2797}
2798
2799xRooNLLVar::xRooHypoSpace xRooNLLVar::hypoSpace(const char *parName, int nPoints, double low, double high,
2800 double alt_value, const xRooFit::Asymptotics::PLLType &pllType)
2801{
2802 if (nPoints < 0) {
2803 // catches case where pyROOT has converted TestStatistic enum to int
2804 int tsType = nPoints;
2805 double alt_val = std::numeric_limits<double>::quiet_NaN();
2807 alt_val = 0;
2808 } else if (tsType == xRooFit::TestStatistic::q0 || tsType == xRooFit::TestStatistic::uncappedq0) {
2809 alt_val = 1;
2810 }
2811
2812 auto out = hypoSpace(parName, pllType, alt_val);
2813
2814 // TODO: things like the physical range and alt value can't be stored on the poi
2815 // because if they change they will change for all hypoSpaces at once, so cannot have
2816 // two hypoSpace with e.g. different physical ranges.
2817 // the hypoSpace should make a copy of them at some point
2818 for (auto p : out.poi()) {
2819 if (tsType == xRooFit::TestStatistic::qmutilde) {
2820 dynamic_cast<RooRealVar *>(p)->setRange("physical", 0, std::numeric_limits<double>::infinity());
2821 Info("xRooNLLVar::hypoSpace", "Setting physical range of %s to [0,inf]", p->GetName());
2822 } else if (dynamic_cast<RooRealVar *>(p)->hasRange("physical")) {
2823 dynamic_cast<RooRealVar *>(p)->removeRange("physical");
2824 Info("xRooNLLVar::hypoSpace", "Setting physical range of %s to [-inf,inf] (i.e. removed range)",
2825 p->GetName());
2826 }
2827 }
2828
2829 // ensure pll type is set explicitly if known at this point
2831 out.fTestStatType = xRooFit::Asymptotics::OneSidedPositive;
2832 } else if (tsType == xRooFit::TestStatistic::uncappedq0) {
2833 out.fTestStatType = xRooFit::Asymptotics::Uncapped;
2834 } else if (tsType == xRooFit::TestStatistic::q0) {
2835 out.fTestStatType = xRooFit::Asymptotics::OneSidedNegative;
2836 }
2837
2838 // in this case the arguments are shifted over by one
2839 if (int(low + 0.5) > 0) {
2840 out.AddPoints(parName, int(low + 0.5), high, alt_value);
2841 } else {
2842 if (!std::isnan(high) && !std::isnan(alt_value) && !(std::isinf(high) && std::isinf(alt_value))) {
2843 for (auto p : out.poi()) {
2844 dynamic_cast<RooRealVar *>(p)->setRange("scan", high, alt_value);
2845 }
2846 }
2847 }
2848 return out;
2849 }
2850
2851 xRooNLLVar::xRooHypoSpace hs = hypoSpace(parName, pllType, alt_value);
2852 if (nPoints > 0)
2853 hs.AddPoints(parName, nPoints, low, high);
2854 else {
2855 if (!std::isnan(low) && !std::isnan(high) && !(std::isinf(low) && std::isinf(high))) {
2856 for (auto p : hs.poi()) {
2857 dynamic_cast<RooRealVar *>(p)->setRange("scan", low, high);
2858 }
2859 }
2860 }
2861 return hs;
2862}
2863
2864xRooNLLVar::xRooHypoSpace xRooNLLVar::hypoSpace(int nPoints, double low, double high, double alt_value,
2865 const xRooFit::Asymptotics::PLLType &pllType)
2866{
2867 auto _poi = std::unique_ptr<RooAbsCollection>(
2868 std::unique_ptr<RooAbsCollection>(pdf()->getVariables())->selectByAttrib("poi", true));
2869 if (_poi->empty())
2870 throw std::runtime_error("You must specify a POI for the hypoSpace");
2871 return hypoSpace(_poi->first()->GetName(), nPoints, low, high, alt_value, pllType);
2872}
2873
2875xRooNLLVar::hypoSpace(const char *parName, const xRooFit::Asymptotics::PLLType &pllType, double alt_value)
2876{
2877 xRooNLLVar::xRooHypoSpace s(parName, parName);
2878
2879 s.AddModel(pdf());
2880 if (strlen(parName)) {
2881 std::unique_ptr<RooAbsCollection> axes(s.pars()->selectByName(parName));
2882 if (axes->empty())
2883 throw std::runtime_error("parameter not found");
2884 axes->setAttribAll("axis", true);
2885 }
2886 /*if (std::unique_ptr<RooAbsCollection>(s.pars()->selectByAttrib("poi", true))->empty()) {
2887 throw std::runtime_error("You must specify at least one POI for the hypoSpace");
2888 }*/
2889 s.fNlls[s.fPdfs.begin()->second] = std::make_shared<xRooNLLVar>(*this);
2890 s.fTestStatType = pllType;
2891
2892 for (auto poi : s.poi()) {
2893 poi->setStringAttribute("altVal", std::isnan(alt_value) ? nullptr : TString::Format("%f", alt_value));
2894 }
2895
2896 return s;
2897}
2898
2900{
2901 if (hypoTestResult) {
2902 return *hypoTestResult;
2903 }
2905 out.SetBackgroundAsAlt(true);
2906 out.SetName(TUUID().AsString());
2907 out.SetTitle(TString::Format("%s HypoPoint", collectionContents(poi()).c_str()));
2908
2909 bool setReadonly = false;
2910 if (nllVar && !nllVar->get()->getAttribute("readOnly")) {
2911 setReadonly = true;
2912 nllVar->get()->setAttribute("readOnly");
2913 }
2914
2915 auto ts_obs = ts_asymp();
2916
2917 out.SetTestStatisticData(ts_obs.first);
2918
2919 // build a ds to hold all fits ... store coords in the globs list of the nullDist
2920 // also need to store at least mu_hat value(s)
2921 RooArgList fitDetails;
2922 RooArgList fitMeta;
2923 fitMeta.addClone(RooCategory(
2924 "pllType", "test statistic type",
2925 {{"TwoSided", 0}, {"OneSidedPositive", 1}, {"OneSidedNegative", 2}, {"Uncapped", 3}, {"Unknown", 4}}));
2926 if (ufit()) {
2927 fitMeta.addClone(ufit()->floatParsFinal());
2928 }
2929 fitMeta.setCatIndex("pllType", int(fPllType));
2930 fitMeta.addClone(RooRealVar("isExpected", "isExpected", int(isExpected)));
2931 fitDetails.addClone(RooCategory("type", "fit type",
2932 {{"ufit", 0},
2933 {"cfit_null", 1},
2934 {"cfit_alt", 2},
2935 {"asimov_ufit", 3},
2936 {"asimov_cfit_null", 4},
2937 {"gen", 5},
2938 {"cfit_lbound", 6}}));
2939 // fitDetails.addClone(RooStringVar("name", "Fit Name", "")); -- not supported properly in ROOT yet
2940 fitDetails.addClone(RooRealVar("status", "status", 0));
2941 fitDetails.addClone(RooRealVar("covQual", "covQual", 0));
2942 fitDetails.addClone(RooRealVar("minNll", "minNll", 0));
2943 fitDetails.addClone(RooRealVar("edm", "edm", 0));
2944 auto fitDS = new RooDataSet("fits", "fit summary data", fitDetails);
2945 fitDS->convertToTreeStore(); // strings not stored properly in vector store, so do convert!
2946
2947 for (int i = 0; i < 7; i++) {
2948 std::shared_ptr<const RooFitResult> fit;
2949 switch (i) {
2950 case 0: fit = ufit(); break;
2951 case 1: fit = cfit_null(); break;
2952 case 2: fit = cfit_alt(); break;
2953 case 3: fit = asimov() ? asimov()->ufit(true) : nullptr; break;
2954 case 4: fit = asimov() ? asimov()->cfit_null(true) : nullptr; break;
2955 case 5: fit = fGenFit; break;
2956 case 6: fit = cfit_lbound(); break;
2957 }
2958 if (fit) {
2959 fitDetails.setCatIndex("type", i);
2960 fitMeta.addClone(RooStringVar(TString::Format("%s.name", fitDetails.getCatLabel("type")),
2961 fitDetails.getCatLabel("type"), fit->GetName()));
2962 // fitDetails.setStringValue("name",fit->GetName());
2963 fitDetails.setRealValue("status", fit->status());
2964 fitDetails.setRealValue("minNll", fit->minNll());
2965 fitDetails.setRealValue("edm", fit->edm());
2966 fitDetails.setRealValue("covQual", fit->covQual());
2967 fitDS->add(fitDetails);
2968 }
2969 }
2970 fitDS->setGlobalObservables(fitMeta);
2971
2972 out.SetFitInfo(fitDS);
2973
2974 RooArgList nullDetails;
2975 RooArgList nullMeta;
2976 nullMeta.addClone(*coords);
2977 nullDetails.addClone(RooRealVar("seed", "Toy Seed", 0));
2978 nullDetails.addClone(RooRealVar("ts", "test statistic value", 0));
2979 nullDetails.addClone(RooRealVar("weight", "weight", 1));
2980 auto nullToyDS = new RooDataSet("nullToys", "nullToys", nullDetails, RooFit::WeightVar("weight"));
2981 nullToyDS->setGlobalObservables(nullMeta);
2982 if (!nullToys.empty()) {
2983
2984 std::vector<double> values;
2985 std::vector<double> weights;
2986 values.reserve(nullToys.size());
2987 weights.reserve(nullToys.size());
2988
2989 for (auto &t : nullToys) {
2990 values.push_back(std::get<1>(t));
2991 weights.push_back(std::get<2>(t));
2992 nullDetails.setRealValue("seed", std::get<0>(t));
2993 nullDetails.setRealValue("ts", std::get<1>(t));
2994 nullToyDS->add(nullDetails, std::get<2>(t));
2995 }
2996 out.SetNullDistribution(new RooStats::SamplingDistribution("null", "Null dist", values, weights, tsTitle()));
2997#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
2998 out.fNullPValue = pNull_toys().first; // technically set above
2999 out.fNullPValueError =
3000 pNull_toys().second; // overrides binomial error used in SamplingDistribution::IntegralAndError
3001#else
3002 out.SetNullPValue(pNull_toys().first); // technically set above
3003 out.SetNullPValueError(
3004 pNull_toys().second); // overrides binomial error used in SamplingDistribution::IntegralAndError
3005#endif
3006 } else {
3007#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
3008 out.fNullPValue = pNull_asymp().first;
3009 out.fNullPValueError = pNull_asymp().second;
3010#else
3011 out.SetNullPValue(pNull_asymp().first);
3012 out.SetNullPValueError(pNull_asymp().second);
3013#endif
3014 }
3015 out.SetNullDetailedOutput(nullToyDS);
3016
3017 if (!altToys.empty()) {
3018 std::vector<double> values;
3019 std::vector<double> weights;
3020 values.reserve(altToys.size());
3021 weights.reserve(altToys.size());
3022 RooArgList altDetails;
3023 RooArgList altMeta;
3024 altDetails.addClone(RooRealVar("seed", "Toy Seed", 0));
3025 altDetails.addClone(RooRealVar("ts", "test statistic value", 0));
3026 altDetails.addClone(RooRealVar("weight", "weight", 1));
3027 auto altToyDS = new RooDataSet("altToys", "altToys", altDetails, RooFit::WeightVar("weight"));
3028 altToyDS->setGlobalObservables(altMeta);
3029 for (auto &t : altToys) {
3030 values.push_back(std::get<1>(t));
3031 weights.push_back(std::get<2>(t));
3032 altDetails.setRealValue("seed", std::get<0>(t));
3033 altDetails.setRealValue("ts", std::get<1>(t));
3034 altToyDS->add(altDetails, std::get<2>(t));
3035 }
3036 out.SetAltDistribution(new RooStats::SamplingDistribution("alt", "Alt dist", values, weights, tsTitle()));
3037 out.SetAltDetailedOutput(altToyDS);
3038#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
3039 out.fAlternatePValue = pAlt_toys().first; // technically set above
3040 out.fAlternatePValueError =
3041 pAlt_toys().second; // overrides binomial error used in SamplingDistribution::IntegralAndError
3042#else
3043 out.SetAltPValue(pAlt_toys().first); // technically set above
3044 out.SetAltPValueError(
3045 pAlt_toys().second); // overrides binomial error used in SamplingDistribution::IntegralAndError
3046#endif
3047
3048 } else {
3049#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00)
3050 out.fAlternatePValue = pAlt_asymp().first;
3051 out.fAlternatePValueError = pAlt_asymp().second;
3052#else
3053 out.SetAltPValue(pAlt_asymp().first);
3054 out.SetAltPValueError(pAlt_asymp().second);
3055#endif
3056 }
3057
3058 if (setReadonly) {
3059 nllVar->get()->setAttribute("readOnly", false);
3060 }
3061
3062 return out;
3063}
3064
3065std::string cling::printValue(const xRooNLLVar::xValueWithError *v)
3066{
3067 if (!v)
3068 return "xValueWithError: nullptr\n";
3069 return Form("%f +/- %f", v->first, v->second);
3070}
3071std::string cling::printValue(const std::map<std::string, xRooNLLVar::xValueWithError> *m)
3072{
3073 if (!m)
3074 return "nullptr\n";
3075 std::string out = "{\n";
3076 for (auto [k, v] : *m) {
3077 out += "\"" + k + "\" => " + printValue(&v) + "\n";
3078 }
3079 out += "}\n";
3080 return out;
3081}
3082
#define SafeDelete(p)
Definition RConfig.hxx:517
#define b(i)
Definition RSha256.hxx:100
#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 h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
GOFOpMode operMode() const
RooAbsData * _data
Pointer to original input dataset.
RooAbsReal * _func
Pointer to original input function.
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
@ kBlue
Definition Rtypes.h:66
#define gDirectory
Definition TDirectory.h:384
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition TError.cxx:218
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
winID h TVirtualViewer3D TVirtualGLPainter p
winID h TVirtualViewer3D vv
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 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 np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
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
@ kCanDelete
Definition TObject.h:367
#define gROOT
Definition TROOT.h:407
static char * Format(const char *format, va_list ap)
Format a string in a circular formatting buffer (using a printf style format descriptor).
Definition TString.cxx:2442
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
R__EXTERN TStyle * gStyle
Definition TStyle.h:436
R__EXTERN TSystem * gSystem
Definition TSystem.h:555
#define gPad
AutoRestorer(const RooAbsCollection &s, xRooNLLVar *nll=nullptr)
RooArgSet fPars
TString fOldTitle
TString fOldName
std::pair< std::shared_ptr< RooAbsData >, std::shared_ptr< const RooAbsCollection > > fOldData
xRooNLLVar * fNll
std::unique_ptr< RooAbsCollection > fSnap
static double k(const IncompatFunc &compatRegions, double pValue, double poiVal, double poiPrimeVal, double sigma_mu=0, double mu_low=-std::numeric_limits< double >::infinity(), double mu_high=std::numeric_limits< double >::infinity())
static int CompatFactor(const IncompatFunc &func, double mu_hat)
static double PValue(const IncompatFunc &compatRegions, double k, double mu, double mu_prime, double sigma_mu=0, double mu_low=-std::numeric_limits< double >::infinity(), double mu_high=std::numeric_limits< double >::infinity())
static std::shared_ptr< const RooFitResult > minimize(RooAbsReal &nll, const std::shared_ptr< ROOT::Fit::FitConfig > &fitConfig=nullptr, const std::shared_ptr< RooLinkedList > &nllOpts=nullptr)
Definition xRooFit.cxx:718
static std::pair< std::shared_ptr< RooAbsData >, std::shared_ptr< const RooAbsCollection > > generateFrom(RooAbsPdf &pdf, const RooFitResult &fr, bool expected=false, int seed=0)
Definition xRooFit.cxx:146
static std::shared_ptr< ROOT::Fit::FitConfig > createFitConfig()
Definition xRooFit.cxx:483
double impact(const char *poi, const char *np, bool up=true, bool prefit=false, bool approx=false)
xRooFitResult ifit(const char *np, bool up, bool prefit=false)
double conditionalError(const char *poi, const char *nps, bool up=true, bool approx=false)
RooArgList ranknp(const char *poi, bool up=true, bool prefit=false, double approxThreshold=std::numeric_limits< double >::infinity())
xRooFitResult(const std::shared_ptr< xRooNode > &in, const std::shared_ptr< xRooNLLVar > &nll=nullptr)
xRooFitResult cfit(const char *poiValues, const char *alias=nullptr)
std::shared_ptr< RooStats::HypoTestResult > hypoTestResult
Definition xRooNLLVar.h:257
std::shared_ptr< const RooFitResult > retrieveFit(int type)
std::vector< std::tuple< int, double, double > > altToys
Definition xRooNLLVar.h:254
std::shared_ptr< const RooAbsCollection > coords
Definition xRooNLLVar.h:242
std::shared_ptr< const RooFitResult > cfit_lbound(bool readOnly=false)
void Draw(Option_t *opt="") override
Default Draw method for all objects.
xValueWithError ts_toys(double nSigma=std::numeric_limits< double >::quiet_NaN())
std::shared_ptr< const RooFitResult > fUfit
Definition xRooNLLVar.h:244
xRooHypoPoint(std::shared_ptr< RooStats::HypoTestResult > htr=nullptr, const RooAbsCollection *_coords=nullptr)
xValueWithError sigma_mu(bool readOnly=false)
xValueWithError pAlt_asymp(double nSigma=std::numeric_limits< double >::quiet_NaN())
std::vector< std::tuple< int, double, double > > nullToys
Definition xRooNLLVar.h:252
std::shared_ptr< xRooHypoPoint > asimov(bool readOnly=false)
xValueWithError pAlt_toys(double nSigma=std::numeric_limits< double >::quiet_NaN())
std::shared_ptr< const RooFitResult > ufit(bool readOnly=false)
void Print(Option_t *opt="") const override
Print TNamed name and title.
std::shared_ptr< const RooFitResult > cfit_null(bool readOnly=false)
std::pair< std::shared_ptr< RooAbsData >, std::shared_ptr< const RooAbsCollection > > data()
xValueWithError pCLs_asymp(double nSigma=std::numeric_limits< double >::quiet_NaN())
std::shared_ptr< const RooFitResult > cfit_alt(bool readOnly=false)
size_t addToys(bool alt, int nToys, int initialSeed=0, double target=std::numeric_limits< double >::quiet_NaN(), double target_nSigma=std::numeric_limits< double >::quiet_NaN(), bool targetCLs=false, double relErrThreshold=2., size_t maxToys=10000)
void addAltToys(int nToys=1, int seed=0, double target=std::numeric_limits< double >::quiet_NaN(), double target_nSigma=std::numeric_limits< double >::quiet_NaN())
void addCLsToys(int nToys=1, int seed=0, double target=std::numeric_limits< double >::quiet_NaN(), double target_nSigma=std::numeric_limits< double >::quiet_NaN())
xValueWithError pX_toys(bool alt, double nSigma=std::numeric_limits< double >::quiet_NaN())
void addNullToys(int nToys=1, int seed=0, double target=std::numeric_limits< double >::quiet_NaN(), double target_nSigma=std::numeric_limits< double >::quiet_NaN())
xValueWithError ts_asymp(double nSigma=std::numeric_limits< double >::quiet_NaN())
xValueWithError pNull_asymp(double nSigma=std::numeric_limits< double >::quiet_NaN())
xValueWithError pNull_toys(double nSigma=std::numeric_limits< double >::quiet_NaN())
std::shared_ptr< RooArgSet > pars() const
Definition xRooNLLVar.h:301
bool AddModel(const xRooNode &pdf, const char *validity="")
std::map< std::shared_ptr< xRooNode >, std::shared_ptr< xRooNLLVar > > fNlls
Definition xRooNLLVar.h:366
int AddPoints(const char *parName, size_t nPoints, double low, double high)
std::set< std::pair< std::shared_ptr< RooArgList >, std::shared_ptr< xRooNode > > > fPdfs
Definition xRooNLLVar.h:368
This xRooNLLVar object has several special methods, e.g.
Definition xRooNLLVar.h:59
std::shared_ptr< RooAbsCollection > fFuncGlobs
Definition xRooNLLVar.h:471
void AddOption(const RooCmdArg &opt)
std::shared_ptr< const RooAbsCollection > fGlobs
Definition xRooNLLVar.h:464
std::shared_ptr< RooLinkedList > fOpts
Definition xRooNLLVar.h:466
std::shared_ptr< RooAbsReal > func() const
std::set< std::string > binnedChannels() const
ROOT::Math::IOptions * fitConfigOptions()
RooConstraintSum * constraintTerm() const
std::shared_ptr< ROOT::Fit::FitConfig > fFitConfig
Definition xRooNLLVar.h:467
xRooHypoSpace hypoSpace(const char *parName, int nPoints, double low, double high, double alt_value=std::numeric_limits< double >::quiet_NaN(), const xRooFit::Asymptotics::PLLType &pllType=xRooFit::Asymptotics::Unknown)
TObject * Scan(const RooArgList &scanPars, const std::vector< std::vector< double > > &coords, const RooArgList &profilePars=RooArgList())
std::shared_ptr< RooAbsCollection > fConstVars
Definition xRooNLLVar.h:470
xRooNLLVar(RooAbsPdf &pdf, const std::pair< RooAbsData *, const RooAbsCollection * > &data, const RooLinkedList &nllOpts=RooLinkedList())
std::shared_ptr< RooAbsPdf > pdf() const
Definition xRooNLLVar.h:423
std::pair< std::shared_ptr< RooAbsData >, std::shared_ptr< const RooAbsCollection > > generate(bool expected=false, int seed=0)
std::pair< std::shared_ptr< RooAbsData >, std::shared_ptr< const RooAbsCollection > > getData() const
double getEntryVal(size_t entry) const
std::shared_ptr< RooAbsCollection > fFuncVars
Definition xRooNLLVar.h:469
double getEntryBinWidth(size_t entry) const
std::shared_ptr< ROOT::Fit::FitConfig > fitConfig()
std::shared_ptr< RooArgSet > pars(bool stripGlobalObs=true) const
std::shared_ptr< RooAbsData > fData
Definition xRooNLLVar.h:463
std::shared_ptr< RooAbsPdf > fPdf
Definition xRooNLLVar.h:462
bool setData(const std::pair< std::shared_ptr< RooAbsData >, std::shared_ptr< const RooAbsCollection > > &_data)
xRooHypoPoint hypoPoint(const char *parName, double value, double alt_value=std::numeric_limits< double >::quiet_NaN(), const xRooFit::Asymptotics::PLLType &pllType=xRooFit::Asymptotics::Unknown)
xRooFitResult minimize(const std::shared_ptr< ROOT::Fit::FitConfig > &=nullptr)
The xRooNode class is designed to wrap over a TObject and provide functionality to aid with interacti...
Definition xRooNode.h:52
xRooNode bins() const
bins of a channel or sample, or channels of a multi-channel pdf
Generic interface for defining configuration options of a numerical algorithm.
Definition IOptions.h:28
void SetValue(const char *name, double val)
generic methods for retrieving options
Definition IOptions.h:42
virtual void SetNamedValue(const char *, const char *)
Definition IOptions.cxx:50
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
TIterator Use servers() and begin()
void Print(Option_t *options=nullptr) const override
Print the object to the defaultPrintStream().
Definition RooAbsArg.h:320
bool dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=nullptr, bool valueOnly=false) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Abstract base class for objects that represent a discrete value that can be set from the outside,...
A space to attach TBranches.
Abstract container object that can hold multiple RooAbsArg objects.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
Storage_t const & get() const
Const access to the underlying stl container.
RooAbsCollection * snapshot(bool deepCopy=true) const
Take a snap shot of current collection contents.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void setAttribAll(const Text_t *name, bool value=true)
Set given attribute in each element of the collection by calling each elements setAttribute() functio...
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
Storage_t::size_type size() const
RooAbsArg * first() const
bool setCatIndex(const char *name, Int_t newVal=0, bool verbose=false)
Set index value of a RooAbsCategoryLValue stored in set with given name to newVal.
bool setRealValue(const char *name, double newVal=0.0, bool verbose=false)
Set value of a RooAbsRealLValue stored in set with given name to newVal No error messages are printed...
virtual RooAbsArg * addClone(const RooAbsArg &var, bool silent=false)
Add a clone of the specified argument to list.
bool selectCommon(const RooAbsCollection &refColl, RooAbsCollection &outColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
const char * getCatLabel(const char *name, const char *defVal="", bool verbose=false) const
Get state name of a RooAbsCategory stored in set with given name.
std::string contentsString() const
Return comma separated list of contained object names as STL string.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual double weight() const =0
virtual double sumEntries() const =0
Return effective number of entries in dataset, i.e., sum all weights.
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
RooArgSet const * getGlobalObservables() const
Returns snapshot of global observables stored in this data.
Definition RooAbsData.h:288
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 ...
virtual void setVal(double value)=0
Set the current value of the object. Needs to be overridden by implementations.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
virtual bool setData(RooAbsData &, bool=true)
Definition RooAbsReal.h:374
bool setData(RooAbsData &data, bool cloneData=true) override
Change dataset that is used to given one.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Object to represent discrete states.
Definition RooCategory.h:28
Named container for two doubles, two integers two object points and three string pointers that can be...
Definition RooCmdArg.h:26
Int_t getInt(Int_t idx) const
Definition RooCmdArg.h:86
TObject * Clone(const char *newName=nullptr) const override
Make a clone of an object using the Streamer facility.
Definition RooCmdArg.h:57
const char * getString(Int_t idx) const
Return string stored in slot idx.
Definition RooCmdArg.h:95
Calculates the sum of the -(log) likelihoods of a set of RooAbsPfs that represent constraint function...
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.
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
Collection class for internal use, storing a collection of RooAbsArg pointers in a doubly linked list...
Int_t GetSize() const
TObject * At(int index) const
Return object stored in sequential position given by index.
static RooMsgService & instance()
Return reference to singleton instance.
StreamConfig & getStream(Int_t id)
Implements a -log(likelihood) calculation from a dataset and a PDF.
Definition RooNLLVar.h:50
Poisson pdf.
Definition RooPoisson.h:19
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:48
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
double getErrorLo() const
Definition RooRealVar.h:67
double getErrorHi() const
Definition RooRealVar.h:68
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
HypoTestResult is a base class for results from hypothesis tests.
This class simply holds a sampling distribution of some test statistic.
A RooAbsArg implementing string values.
Draw all kinds of Arrows.
Definition TArrow.h:29
virtual TArrow * DrawArrow(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Float_t arrowsize=0, Option_t *option="")
Draw this arrow with new coordinates.
Definition TArrow.cxx:135
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:42
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition TAttLine.h:43
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
static TCanvas * MakeDefCanvas()
Static function to build a default canvas.
Definition TCanvas.cxx:1516
const char * AsString() const
Return the date & time as a string (ctime() format).
Definition TDatime.cxx:102
Describe directory structure in memory.
Definition TDirectory.h:45
virtual TKey * FindKeyAny(const char *) const
Definition TDirectory.h:198
Class to handle efficiency histograms.
Definition TEfficiency.h:29
void FillWeighted(Bool_t bPassed, Double_t weight, Double_t x, Double_t y=0, Double_t z=0)
This function is used for filling the two histograms with a weight.
Double_t GetEfficiencyErrorUp(Int_t bin) const
Returns the upper error on the efficiency in the given global bin.
void Fill(Bool_t bPassed, Double_t x, Double_t y=0, Double_t z=0)
This function is used for filling the two histograms.
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition TGraph2D.h:41
Int_t GetN() const
Definition TGraph2D.h:120
virtual void SetPoint(Int_t point, Double_t x, Double_t y, Double_t z)
Sets point number n.
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2325
Int_t GetN() const
Definition TGraph.h:131
void SetName(const char *name="") override
Set graph name.
Definition TGraph.cxx:2364
void Draw(Option_t *chopt="") override
Draw this graph with its current attributes.
Definition TGraph.cxx:814
void SetTitle(const char *title="") override
Change (i.e.
Definition TGraph.cxx:2380
void SetNameTitle(const char *name="", const char *title="") override
Set graph name and title.
Definition TGraph.cxx:2400
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:670
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
virtual Double_t GetMinimum(Double_t minval=-FLT_MAX) const
Return minimum value larger than minval of bins in the range, unless the value has been overridden by...
Definition TH1.cxx:8626
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
void Draw(Option_t *option="") override
Draw this legend with its current attributes.
Definition TLegend.cxx:425
Use the TLine constructor to create a simple line.
Definition TLine.h:22
virtual TLine * DrawLine(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Draw this line with new coordinates.
Definition TLine.cxx:103
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
void Clear(Option_t *option="") override
Set name and title to empty strings ("").
Definition TNamed.cxx:64
Mother of all ROOT objects.
Definition TObject.h:41
virtual const char * GetName() const
Returns name of object.
Definition TObject.cxx:438
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition TObject.cxx:223
virtual void Delete(Option_t *option="")
Delete this object.
Definition TObject.cxx:248
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:780
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition TRandom.cxx:615
virtual UInt_t Integer(UInt_t imax)
Returns a random integer uniformly distributed on the interval [ 0, imax-1 ].
Definition TRandom.cxx:361
Regular expression class.
Definition TRegexp.h:31
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
Int_t Atoi() const
Return integer value of string.
Definition TString.cxx:1988
Double_t Atof() const
Return floating-point value contained in string.
Definition TString.cxx:2054
Bool_t IsFloat() const
Returns kTRUE if string contains a floating point or integer number.
Definition TString.cxx:1858
const char * Data() const
Definition TString.h:376
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition TString.h:704
Bool_t BeginsWith(const char *s, ECaseCompare cmp=kExact) const
Definition TString.h:623
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2378
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition TString.h:651
Float_t GetPadRightMargin() const
Definition TStyle.h:214
Float_t GetPadLeftMargin() const
Definition TStyle.h:213
Float_t GetPadBottomMargin() const
Definition TStyle.h:211
Float_t GetPadTopMargin() const
Definition TStyle.h:212
virtual Bool_t ProcessEvents()
Process pending events (GUI, timers, sockets).
Definition TSystem.cxx:416
This class defines a UUID (Universally Unique IDentifier), also known as GUIDs (Globally Unique IDent...
Definition TUUID.h:42
TDatime GetTime() const
Get time from UUID.
Definition TUUID.cxx:670
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition TVirtualPad.h:51
virtual TList * GetListOfPrimitives() const =0
virtual TObject * GetPrimitive(const char *name) const =0
RooCmdArg WeightVar(const char *name="weight", bool reinterpretAsWeight=false)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg GlobalObservablesSource(const char *sourceName)
double gaussian_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
const Int_t n
Definition legend1.C:16
TGraphErrors * gr
Definition legend1.C:25
double gaussian_cdf(double x, double sigma=1, double x0=0)
Alternative name for same function.
@ NumIntegration
Double_t Prob(Double_t chi2, Int_t ndf)
Computation of the probability for a certain Chi-squared (chi2) and number of degrees of freedom (ndf...
Definition TMath.cxx:637
Double_t Poisson(Double_t x, Double_t par)
Computes the Poisson distribution function for (x,par).
Definition TMath.cxx:587
#define BEGIN_XROOFIT_NAMESPACE
Definition Config.h:24
#define END_XROOFIT_NAMESPACE
Definition Config.h:25
static const char * what
Definition stlLoader.cc:5
void removeTopic(RooFit::MsgTopic oldTopic)
th1 Draw()
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4
std::string collectionContents(const RooAbsCollection &coll)
#define GETWS(a)
#define GETWSSETS(w)