Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
AsymptoticCalculator.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer, Sven Kreiss 23/05/10
3/*************************************************************************
4 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11/** \class RooStats::AsymptoticCalculator
12 \ingroup Roostats
13
14Hypothesis Test Calculator based on the asymptotic formulae for the profile
15likelihood ratio.
16
17It performs hypothesis tests using the asymptotic formula for the profile likelihood, and
18uses the Asimov data set to compute expected significances or limits.
19
20See G. Cowan, K. Cranmer, E. Gross and O. Vitells: Asymptotic formulae for
21likelihood- based tests of new physics. Eur. Phys. J., C71:1–19, 2011.
22It provides methods to perform hypothesis tests using the likelihood function,
23and computes the \f$p\f$-values for the null and the alternate hypothesis using the asymptotic
24formulae for the profile likelihood ratio described in the given paper.
25
26The calculator provides methods to produce the Asimov dataset, *i.e.* a dataset
27generated where the observed values are equal to the expected ones.
28The Asimov data set is then used to compute the observed asymptotic \f$p\f$-value for
29the alternate hypothesis and the asymptotic expected \f$p\f$-values.
30
31The asymptotic formulae are valid only for one POI (parameter of interest). So
32the calculator works only for one-dimensional (one POI) models.
33If more than one POI exists, only the first one is used.
34
35The calculator can generate Asimov datasets from two kinds of PDFs:
36- "Counting" distributions: RooPoisson, RooGaussian, or products of RooPoissons.
37- Extended, *i.e.* number of events can be read off from extended likelihood term.
38*/
39
40
45
46#include "RooArgSet.h"
47#include "RooArgList.h"
48#include "RooProdPdf.h"
49#include "RooSimultaneous.h"
50#include "RooDataSet.h"
51#include "RooCategory.h"
52#include "RooRealVar.h"
53#include "RooMinimizer.h"
54#include "RooFitResult.h"
56#include "RooPoisson.h"
57#include "RooUniform.h"
58#include "RooGamma.h"
59#include "RooGaussian.h"
60#include "RooBifurGauss.h"
61#include "RooLognormal.h"
62#include "RooDataHist.h"
63#include <cmath>
64#include <typeinfo>
65
68
69#include "TStopwatch.h"
70
71using namespace RooStats;
72using std::cout, std::endl, std::string, std::unique_ptr;
73
75
77
78////////////////////////////////////////////////////////////////////////////////
79/// set print level (static function)
80///
81/// - 0 minimal,
82/// - 1 normal,
83/// - 2 debug
84
86 fgPrintLevel = level;
87}
88
89////////////////////////////////////////////////////////////////////////////////
90/// constructor for asymptotic calculator from Data set and ModelConfig
91
94 const ModelConfig &altModel,
95 const ModelConfig &nullModel, bool nominalAsimov) :
96 HypoTestCalculatorGeneric(data, altModel, nullModel, nullptr),
97 fOneSided(false), fOneSidedDiscovery(false), fNominalAsimov(nominalAsimov),
98 fUseQTilde(-1),
99 fNLLObs(0), fNLLAsimov(0),
100 fAsimovData(nullptr)
101{
102 if (!Initialize()) return;
103
104 int verbose = fgPrintLevel;
105 // try to guess default configuration
106 // (this part should be only in constructor because the null snapshot might change during HypoTestInversion
107 const RooArgSet * nullSnapshot = GetNullModel()->GetSnapshot();
108 assert(nullSnapshot);
109 RooRealVar * muNull = dynamic_cast<RooRealVar*>(nullSnapshot->first() );
110 assert(muNull);
111 if (muNull->getVal() == muNull->getMin()) {
112 fOneSidedDiscovery = true;
113 if (verbose > 0)
114 oocoutI(nullptr,InputArguments) << "AsymptotiCalculator: Minimum of POI is " << muNull->getMin() << " corresponds to null snapshot - default configuration is one-sided discovery formulae " << std::endl;
115 }
116
117}
118
119////////////////////////////////////////////////////////////////////////////////
120/// Initialize the calculator
121/// The initialization will perform a global fit of the model to the data
122/// and build an Asimov data set.
123/// It will then also fit the model to the Asimov data set to find the likelihood value
124/// of the Asimov data set
125/// nominalAsimov is an option for using Asimov data set obtained using nominal nuisance parameter values
126/// By default the nuisance parameters are fitted to the data
127/// NOTE: If a fit has been done before, one for speeding up could set all the initial parameters
128/// to the fit value and in addition set the null snapshot to the best fit
129
131
132 int verbose = fgPrintLevel;
133 if (verbose >= 0)
134 oocoutP(nullptr,Eval) << "AsymptoticCalculator::Initialize...." << std::endl;
135
136
137 RooAbsPdf * nullPdf = GetNullModel()->GetPdf();
138 if (!nullPdf) {
139 oocoutE(nullptr,InputArguments) << "AsymptoticCalculator::Initialize - ModelConfig has not a pdf defined" << std::endl;
140 return false;
141 }
142 RooAbsData * obsData = const_cast<RooAbsData *>(GetData() );
143 if (!obsData ) {
144 oocoutE(nullptr,InputArguments) << "AsymptoticCalculator::Initialize - data set has not been defined" << std::endl;
145 return false;
146 }
147 RooAbsData & data = *obsData;
148
149
150
152 if (!poi || poi->empty()) {
153 oocoutE(nullptr,InputArguments) << "AsymptoticCalculator::Initialize - ModelConfig has not POI defined." << endl;
154 return false;
155 }
156 if (poi->size() > 1) {
157 oocoutW(nullptr,InputArguments) << "AsymptoticCalculator::Initialize - ModelConfig has more than one POI defined \n\t"
158 << "The asymptotic calculator works for only one POI - consider as POI only the first parameter"
159 << std::endl;
160 }
161
162
163 // This will set the poi value to the null snapshot value in the ModelConfig
164 const RooArgSet * nullSnapshot = GetNullModel()->GetSnapshot();
165 if(nullSnapshot == nullptr || nullSnapshot->empty()) {
166 oocoutE(nullptr,InputArguments) << "AsymptoticCalculator::Initialize - Null model needs a snapshot. Set using modelconfig->SetSnapshot(poi)." << endl;
167 return false;
168 }
169
170 // GetNullModel()->Print();
171 // printf("ASymptotic calc: null snapshot\n");
172 // nullSnapshot->Print("v");
173 // printf("PDF variables " );
174 // nullPdf->getVariables()->Print("v");
175
176 // keep snapshot for the initial parameter values (need for nominal Asimov)
177 RooArgSet nominalParams;
178 std::unique_ptr<RooArgSet> allParams{nullPdf->getParameters(data)};
179 RemoveConstantParameters(&*allParams);
180 if (fNominalAsimov) {
181 allParams->snapshot(nominalParams);
182 }
186
187 // evaluate the unconditional nll for the full model on the observed data
188 if (verbose >= 0)
189 oocoutP(nullptr,Eval) << "AsymptoticCalculator::Initialize - Find best unconditional NLL on observed data" << endl;
191 // fill also snapshot of best poi
192 poi->snapshot(fBestFitPoi);
193 RooRealVar * muBest = dynamic_cast<RooRealVar*>(fBestFitPoi.first());
194 assert(muBest);
195 if (verbose >= 0)
196 oocoutP(nullptr,Eval) << "Best fitted POI value = " << muBest->getVal() << " +/- " << muBest->getError() << std::endl;
197 // keep snapshot of all best fit parameters
198 allParams->snapshot(fBestFitParams);
199
200 // compute Asimov data set for the background (alt poi ) value
201 const RooArgSet * altSnapshot = GetAlternateModel()->GetSnapshot();
202 if(altSnapshot == nullptr || altSnapshot->empty()) {
203 oocoutE(nullptr,InputArguments) << "Alt (Background) model needs a snapshot. Set using modelconfig->SetSnapshot(poi)." << endl;
204 return false;
205 }
206
207 RooArgSet poiAlt(*altSnapshot); // this is the poi snapshot of B (i.e. for mu=0)
208
209 oocoutP(nullptr,Eval) << "AsymptoticCalculator: Building Asimov data Set" << endl;
210
211 // check that in case of binned models the n number of bins of the observables are consistent
212 // with the number of bins in the observed data
213 // This number will be used for making the Asimov data set so it will be more consistent with the
214 // observed data
215 int prevBins = 0;
216 RooRealVar * xobs = nullptr;
217 if (GetNullModel()->GetObservables() && GetNullModel()->GetObservables()->size() == 1 ) {
218 xobs = static_cast<RooRealVar*>((GetNullModel()->GetObservables())->first());
219 if (data.IsA() == RooDataHist::Class() ) {
220 if (data.numEntries() != xobs->getBins() ) {
221 prevBins = xobs->getBins();
222 oocoutW(nullptr,InputArguments) << "AsymptoticCalculator: number of bins in " << xobs->GetName() << " are different than data bins "
223 << " set the same data bins " << data.numEntries() << " in range "
224 << " [ " << xobs->getMin() << " , " << xobs->getMax() << " ]" << std::endl;
225 xobs->setBins(data.numEntries());
226 }
227 }
228 }
229
230 if (!fNominalAsimov) {
231 if (verbose >= 0)
232 oocoutI(nullptr,InputArguments) << "AsymptoticCalculator: Asimov data will be generated using fitted nuisance parameter values" << endl;
233 RooArgSet * tmp = (RooArgSet*) poiAlt.snapshot();
235 }
236
237 else {
238 // assume use current value of nuisance as nominal ones
239 if (verbose >= 0)
240 oocoutI(nullptr,InputArguments) << "AsymptoticCalculator: Asimovdata set will be generated using nominal (current) nuisance parameter values" << endl;
241 nominalParams.assign(poiAlt); // set poi to alt value but keep nuisance at the nominal one
243 }
244
245 if (!fAsimovData) {
246 oocoutE(nullptr,InputArguments) << "AsymptoticCalculator: Error : Asimov data set could not be generated " << endl;
247 return false;
248 }
249
250 // set global observables to their Asimov values
251 RooArgSet globObs;
252 RooArgSet globObsSnapshot;
253 if (GetNullModel()->GetGlobalObservables() ) {
254 globObs.add(*GetNullModel()->GetGlobalObservables());
255 assert(globObs.size() == fAsimovGlobObs.size() );
256 // store previous snapshot value
257 globObs.snapshot(globObsSnapshot);
258 globObs.assign(fAsimovGlobObs);
259 }
260
261
262 // evaluate the likelihood. Since we use on Asimov data , conditional and unconditional values should be the same
263 // do conditional fit since is faster
264
265 RooRealVar * muAlt = static_cast<RooRealVar*>(poiAlt.first());
266 assert(muAlt);
267 if (verbose >= 0) {
268 oocoutP(nullptr, Eval)
269 << "AsymptoticCalculator::Initialize Find best conditional NLL on ASIMOV data set for given alt POI ( "
270 << muAlt->GetName() << " ) = " << muAlt->getVal() << std::endl;
271 }
272
274 // for unconditional fit
275 //fNLLAsimov = EvaluateNLL( *nullPdf, *fAsimovData);
276 //poi->Print("v");
277
278 // restore previous value
279 globObs.assign(globObsSnapshot);
280
281 // restore number of bins
282 if (prevBins > 0 && xobs) xobs->setBins(prevBins);
283
284 fIsInitialized = true;
285 return true;
286}
287
288////////////////////////////////////////////////////////////////////////////////
289
291{
292 int verbose = fgPrintLevel;
293
294 RooAbsPdf &pdf = *modelConfig.GetPdf();
295
298
299
300 std::unique_ptr<RooArgSet> allParams{pdf.getParameters(data)};
302 // add constraint terms for all non-constant parameters
303
304 // need to call constrain for RooSimultaneous until stripDisconnected problem fixed
305 auto& config = GetGlobalRooStatsConfig();
306 std::unique_ptr<RooAbsReal> nll{modelConfig.createNLL(data, RooFit::Constrain(*allParams), RooFit::Offset(config.useLikelihoodOffset))};
307
308 std::unique_ptr<RooArgSet> attachedSet{nll->getVariables()};
309
310 // if poi are specified - do a conditional fit
311 RooArgSet paramsSetConstant;
312 // support now only one POI
313 if (poiSet && !poiSet->empty()) {
314 RooRealVar * muTest = static_cast<RooRealVar*> (poiSet->first());
315 RooRealVar * poiVar = dynamic_cast<RooRealVar*>(attachedSet->find( muTest->GetName() ) );
316 if (poiVar && !poiVar->isConstant() ) {
317 poiVar->setVal( muTest->getVal() );
318 poiVar->setConstant();
319 paramsSetConstant.add(*poiVar);
320 }
321 if (poiSet->size() > 1)
322 std::cout << "Model with more than one POI are not supported - ignore extra parameters, consider only first one" << std::endl;
323
324
325
326 // This for more than one POI (not yet supported)
327 //
328 // RooLinkedListIter it = poiSet->iterator();
329 // RooRealVar* tmpPar = nullptr, *tmpParA=nullptr;
330 // while((tmpPar = (RooRealVar*)it.Next())){
331 // tmpParA = ((RooRealVar*)attachedSet->find(tmpPar->GetName()));
332 // tmpParA->setVal( tmpPar->getVal() );
333 // if (!tmpParA->isConstant() ) {
334 // tmpParA->setConstant();
335 // paramsSetConstant.add(*tmpParA);
336 // }
337 // }
338
339 // check if there are non-const parameters so it is worth to do the minimization
340
341 }
342
343 TStopwatch tw;
344 tw.Start();
345 double val = -1;
346
347 //check if needed to skip the fit
348 RooArgSet nllParams(*attachedSet);
350 bool skipFit = (nllParams.empty());
351
352 if (skipFit) {
353 val = nll->getVal(); // just evaluate nll in conditional fits with model without nuisance params
354 } else {
355
356 int minimPrintLevel = verbose;
357
358 RooMinimizer minim(*nll);
360 minim.setStrategy( strategy);
361 minim.setEvalErrorWall(config.useEvalErrorWall);
362 // use tolerance - but never smaller than 1 (default in RooMinimizer)
364 tol = std::max(tol,1.0); // 1.0 is the minimum value used in RooMinimizer
365 minim.setEps( tol );
366 //LM: RooMinimizer.setPrintLevel has +1 offset - so subtract here -1
367 minim.setPrintLevel(minimPrintLevel-1);
368 int status = -1;
369 minim.optimizeConst(2);
370 TString minimizer = ""; // empty string to take RooMinimizer default initially
372
373 if (verbose > 0) {
374 std::cout << "AsymptoticCalculator::EvaluateNLL ........ using " << minimizer << " / " << algorithm
375 << " with strategy " << strategy << " and tolerance " << tol << std::endl;
376 }
377
378 for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
379 // status = minim.minimize(fMinimizer, ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str());
380 status = minim.minimize(minimizer, algorithm);
381 // RooMinimizer::minimize returns -1 when the fit fails
382 if (status >= 0) {
383 break;
384 } else {
385 if (tries == 1) {
386 printf(" ----> Doing a re-scan first\n");
387 minim.minimize(minimizer,"Scan");
388 }
389 if (tries == 2) {
391 printf(" ----> trying with strategy = 1\n");
392 minim.setStrategy(1);
393 }
394 else
395 tries++; // skip this trial if strategy is already 1
396 }
397 if (tries == 3) {
398 printf(" ----> trying with improve\n");
399 minimizer = "Minuit";
400 algorithm = "migradimproved";
401 }
402 }
403 }
404
405 std::unique_ptr<RooFitResult> result;
406
407 // ignore errors in Hesse or in Improve and also when matrix was made pos def (status returned = 1)
408 if (status >= 0) {
409 result = std::unique_ptr<RooFitResult>{minim.save()};
410 }
411 if (result){
412 if (!RooStats::IsNLLOffset()) {
413 val = result->minNll();
414 } else {
415 bool previous = RooAbsReal::hideOffset();
417 val = nll->getVal();
418 if (!previous) RooAbsReal::setHideOffset(false) ;
419 }
420
421 }
422 else {
423 oocoutE(nullptr,Fitting) << "FIT FAILED !- return a NaN NLL " << std::endl;
424 val = TMath::QuietNaN();
425 }
426
427 minim.optimizeConst(false);
428 }
429
430 double muTest = 0;
431 if (verbose > 0) {
432 std::cout << "AsymptoticCalculator::EvaluateNLL - value = " << val;
433 if (poiSet) {
434 muTest = ( static_cast<RooRealVar*>(poiSet->first()) )->getVal();
435 std::cout << " for poi fixed at = " << muTest;
436 }
437 if (!skipFit) {
438 std::cout << "\tfit time : ";
439 tw.Print();
440 }
441 else
442 std::cout << std::endl;
443 }
444
445 // reset the parameter free which where set as constant
446 if (poiSet && !paramsSetConstant.empty()) SetAllConstant(paramsSetConstant,false);
447
448
449 if (verbose < 2) RooMsgService::instance().setGlobalKillBelow(msglevel);
450
451 return val;
452}
453
454////////////////////////////////////////////////////////////////////////////////
455/// It performs an hypothesis tests using the likelihood function
456/// and computes the p values for the null and the alternate using the asymptotic
457/// formulae for the profile likelihood ratio.
458/// See G. Cowan, K. Cranmer, E. Gross and O. Vitells.
459/// Asymptotic formulae for likelihood- based tests of new physics. Eur. Phys. J., C71:1–19, 2011.
460/// The formulae are valid only for one POI. If more than one POI exists consider as POI only the
461/// first one
462
464 int verbose = fgPrintLevel;
465
466 // re-initialized the calculator in case it is needed (pdf or data modified)
467 if (!fIsInitialized) {
468 if (!Initialize() ) {
469 oocoutE(nullptr,InputArguments) << "AsymptoticCalculator::GetHypoTest - Error initializing Asymptotic calculator - return nullptr result " << endl;
470 return nullptr;
471 }
472 }
473
474 if (!fAsimovData) {
475 oocoutE(nullptr,InputArguments) << "AsymptoticCalculator::GetHypoTest - Asimov data set has not been generated - return nullptr result " << endl;
476 return nullptr;
477 }
478
479 assert(GetNullModel() );
480 assert(GetData() );
481
482 RooAbsPdf * nullPdf = GetNullModel()->GetPdf();
483 assert(nullPdf);
484
485 // make conditional fit on null snapshot of poi
486
487 const RooArgSet * nullSnapshot = GetNullModel()->GetSnapshot();
488 assert(nullSnapshot && !nullSnapshot->empty());
489
490 // use as POI the nullSnapshot
491 // if more than one POI exists, consider only the first one
492 RooArgSet poiTest(*nullSnapshot);
493
494 if (poiTest.size() > 1) {
495 oocoutW(nullptr,InputArguments) << "AsymptoticCalculator::GetHypoTest: snapshot has more than one POI - assume as POI first parameter " << std::endl;
496 }
497
498 std::unique_ptr<RooArgSet> allParams{nullPdf->getParameters(*GetData() )};
499 allParams->assign(fBestFitParams);
500
501 // set the one-side condition
502 // (this works when we have only one params of interest
503 RooRealVar * muHat = dynamic_cast<RooRealVar*> ( fBestFitPoi.first() );
504 assert(muHat && "no best fit parameter defined");
505 RooRealVar * muTest = dynamic_cast<RooRealVar*> ( nullSnapshot->find(muHat->GetName() ) );
506 assert(muTest && "poi snapshot is not existing");
507
508
509
510 if (verbose> 0) {
511 std::cout << std::endl;
512 oocoutI(nullptr,Eval) << "AsymptoticCalculator::GetHypoTest: - perform an hypothesis test for POI ( " << muTest->GetName() << " ) = " << muTest->getVal() << std::endl;
513 oocoutP(nullptr,Eval) << "AsymptoticCalculator::GetHypoTest - Find best conditional NLL on OBSERVED data set ..... " << std::endl;
514 }
515
516 // evaluate the conditional NLL on the observed data for the snapshot value
517 double condNLL = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData&>(*GetData()), &poiTest);
518
519 double qmu = 2.*(condNLL - fNLLObs);
520
521
522
523 if (verbose > 0)
524 oocoutP(nullptr,Eval) << "\t OBSERVED DATA : qmu = " << qmu << " condNLL = " << condNLL << " uncond " << fNLLObs << std::endl;
525
526
527 // this tolerance is used to avoid having negative qmu due to numerical errors
528 double tol = 2.E-3 * std::max(1.,ROOT::Math::MinimizerOptions::DefaultTolerance());
529 if (qmu < -tol || TMath::IsNaN(fNLLObs) ) {
530
531 if (qmu < 0) {
532 oocoutW(nullptr,Minimization) << "AsymptoticCalculator: Found a negative value of the qmu - retry to do the unconditional fit "
533 << std::endl;
534 } else {
535 oocoutW(nullptr, Minimization)
536 << "AsymptoticCalculator: unconditional fit failed before - retry to do it now " << std::endl;
537 }
538
539 double nll = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData&>(*GetData()));
540
541 if (nll < fNLLObs || (TMath::IsNaN(fNLLObs) && !TMath::IsNaN(nll) ) ) {
542 oocoutW(nullptr,Minimization) << "AsymptoticCalculator: Found a better unconditional minimum "
543 << " old NLL = " << fNLLObs << " old muHat " << muHat->getVal() << std::endl;
544
545 // update values
546 fNLLObs = nll;
548 assert(poi);
550 poi->snapshot(fBestFitPoi);
551 // restore also muHad since previous pointer has been deleted
552 muHat = dynamic_cast<RooRealVar*> ( fBestFitPoi.first() );
553 assert(muHat);
554
555 oocoutW(nullptr,Minimization) << "AsymptoticCalculator: New minimum found for "
556 << " NLL = " << fNLLObs << " muHat " << muHat->getVal() << std::endl;
557
558
559 qmu = 2.*(condNLL - fNLLObs);
560
561 if (verbose > 0)
562 oocoutP(nullptr,Eval) << "After unconditional refit, new qmu value is " << qmu << std::endl;
563
564 }
565 }
566
567 if (qmu < -tol ) {
568 oocoutE(nullptr,Minimization) << "AsymptoticCalculator: qmu is still < 0 for mu = "
569 << muTest->getVal() << " return a dummy result "
570 << std::endl;
571 return new HypoTestResult();
572 }
573 if (TMath::IsNaN(qmu) ) {
574 oocoutE(nullptr,Minimization) << "AsymptoticCalculator: failure in fitting for qmu or qmuA "
575 << muTest->getVal() << " return a dummy result "
576 << std::endl;
577 return new HypoTestResult();
578 }
579
580
581
582
583
584 // compute conditional ML on Asimov data set
585 // (need to const cast because it uses fitTo which is a non const method
586 // RooArgSet asimovGlobObs;
587 // RooAbsData * asimovData = (const_cast<AsymptoticCalculator*>(this))->MakeAsimovData( poi, asimovGlobObs);
588 // set global observables to their Asimov values
589 RooArgSet globObs;
590 RooArgSet globObsSnapshot;
591 if (GetNullModel()->GetGlobalObservables() ) {
592 globObs.add(*GetNullModel()->GetGlobalObservables());
593 // store previous snapshot value
594 globObs.snapshot(globObsSnapshot);
595 globObs.assign(fAsimovGlobObs);
596 }
597
598
599 if (verbose > 0) oocoutP(nullptr,Eval) << "AsymptoticCalculator::GetHypoTest -- Find best conditional NLL on ASIMOV data set .... " << std::endl;
600
601 double condNLL_A = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiTest);
602
603
604 double qmu_A = 2.*(condNLL_A - fNLLAsimov );
605
606 if (verbose > 0)
607 oocoutP(nullptr,Eval) << "\t ASIMOV data qmu_A = " << qmu_A << " condNLL = " << condNLL_A << " uncond " << fNLLAsimov << std::endl;
608
609 if (qmu_A < -tol || TMath::IsNaN(fNLLAsimov) ) {
610
611 if (qmu_A < 0) {
612 oocoutW(nullptr, Minimization)
613 << "AsymptoticCalculator: Found a negative value of the qmu Asimov- retry to do the unconditional fit "
614 << std::endl;
615 } else {
616 oocoutW(nullptr, Minimization)
617 << "AsymptoticCalculator: Fit failed for unconditional the qmu Asimov- retry unconditional fit "
618 << std::endl;
619 }
620
621 double nll = EvaluateNLL(*GetNullModel(), *fAsimovData);
622
623 if (nll < fNLLAsimov || (TMath::IsNaN(fNLLAsimov) && !TMath::IsNaN(nll) )) {
624 oocoutW(nullptr,Minimization) << "AsymptoticCalculator: Found a better unconditional minimum for Asimov data set"
625 << " old NLL = " << fNLLAsimov << std::endl;
626
627 // update values
628 fNLLAsimov = nll;
629
630 oocoutW(nullptr,Minimization) << "AsymptoticCalculator: New minimum found for "
631 << " NLL = " << fNLLAsimov << std::endl;
632 qmu_A = 2.*(condNLL_A - fNLLAsimov);
633
634 if (verbose > 0)
635 oocoutP(nullptr,Eval) << "After unconditional Asimov refit, new qmu_A value is " << qmu_A << std::endl;
636
637 }
638 }
639
640 if (qmu_A < - tol) {
641 oocoutE(nullptr,Minimization) << "AsymptoticCalculator: qmu_A is still < 0 for mu = "
642 << muTest->getVal() << " return a dummy result "
643 << std::endl;
644 return new HypoTestResult();
645 }
646 if (TMath::IsNaN(qmu) ) {
647 oocoutE(nullptr,Minimization) << "AsymptoticCalculator: failure in fitting for qmu or qmuA "
648 << muTest->getVal() << " return a dummy result "
649 << std::endl;
650 return new HypoTestResult();
651 }
652
653
654 // restore previous value of global observables
655 globObs.assign(globObsSnapshot);
656
657 // now we compute p-values using the asymptotic formulae
658 // described in the paper
659 // Cowan et al, Eur.Phys.J. C (2011) 71:1554
660
661 // first try to guess automatically if needed to use qtilde (or ttilde in case of two sided)
662 // if explicitly fUseQTilde this was not set
663 // qtilde is in this case used if poi is bounded at the value of the alt hypothesis
664 // for Qtilde (need to distinguish case when qmu > qmuA = mu^2/ sigma^2)
665 // (see Cowan et al, Eur.Phys.J. C(2011) 71:1554 paper equations 64 and 65
666 // (remember qmu_A = mu^2/sigma^2 )
667 bool useQTilde = false;
668 // default case (check if poi is limited or not to a zero value)
669 if (!fOneSidedDiscovery) { // qtilde is not a discovery test
670 if (fUseQTilde == -1 && !fOneSidedDiscovery) {
671 // alternate snapshot is value for which background is zero (for limits)
672 RooRealVar * muAlt = dynamic_cast<RooRealVar*>(GetAlternateModel()->GetSnapshot()->first() );
673 // null snapshot is value for which background is zero (for discovery)
674 //RooRealVar * muNull = dynamic_cast<RooRealVar*>(GetNullModel()->GetSnapshot()->first() );
675 assert(muAlt != nullptr );
676 if (muTest->getMin() == muAlt->getVal() ) {
677 fUseQTilde = 1;
678 oocoutI(nullptr,InputArguments) << "Minimum of POI is " << muTest->getMin() << " corresponds to alt snapshot - using qtilde asymptotic formulae " << std::endl;
679 } else {
680 fUseQTilde = 0;
681 oocoutI(nullptr,InputArguments) << "Minimum of POI is " << muTest->getMin() << " is different to alt snapshot " << muAlt->getVal()
682 << " - using standard q asymptotic formulae " << std::endl;
683 }
684 }
685 useQTilde = fUseQTilde;
686 }
687
688
689 //check for one side condition (remember this is valid only for one poi)
690 if (fOneSided ) {
691 if ( muHat->getVal() > muTest->getVal() ) {
692 oocoutI(nullptr,Eval) << "Using one-sided qmu - setting qmu to zero muHat = " << muHat->getVal()
693 << " muTest = " << muTest->getVal() << std::endl;
694 qmu = 0;
695 }
696 }
697 if (fOneSidedDiscovery ) {
698 if ( muHat->getVal() < muTest->getVal() ) {
699 oocoutI(nullptr,Eval) << "Using one-sided discovery qmu - setting qmu to zero muHat = " << muHat->getVal()
700 << " muTest = " << muTest->getVal() << std::endl;
701 qmu = 0;
702 }
703 }
704
705 // fix for negative qmu values due to numerical errors
706 if (qmu < 0 && qmu > -tol) qmu = 0;
707 if (qmu_A < 0 && qmu_A > -tol) qmu_A = 0;
708
709 // asymptotic formula for pnull and from paper Eur.Phys.J C 2011 71:1554
710 // we have 4 different cases:
711 // t(mu), t_tilde(mu) for the 2-sided
712 // q(mu) and q_tilde(mu) for the one -sided test statistics
713
714 double pnull = -1;
715 double palt = -1;
716
717 // asymptotic formula for pnull (for only one POI)
718 // From fact that qmu is a chi2 with ndf=1
719
720 double sqrtqmu = (qmu > 0) ? std::sqrt(qmu) : 0;
721 double sqrtqmu_A = (qmu_A > 0) ? std::sqrt(qmu_A) : 0;
722
723
725 // for one-sided PL (q_mu : equations 56,57)
726 if (verbose>2) {
727 if (fOneSided) {
728 oocoutI(nullptr,Eval) << "Using one-sided limit asymptotic formula (qmu)" << endl;
729 } else {
730 oocoutI(nullptr, Eval) << "Using one-sided discovery asymptotic formula (q0)" << endl;
731 }
732 }
733 pnull = ROOT::Math::normal_cdf_c( sqrtqmu, 1.);
734 palt = ROOT::Math::normal_cdf( sqrtqmu_A - sqrtqmu, 1.);
735 }
736 else {
737 // for 2-sided PL (t_mu : equations 35,36 in asymptotic paper)
738 if (verbose > 2) oocoutI(nullptr,Eval) << "Using two-sided asymptotic formula (tmu)" << endl;
739 pnull = 2.*ROOT::Math::normal_cdf_c( sqrtqmu, 1.);
740 palt = ROOT::Math::normal_cdf_c( sqrtqmu + sqrtqmu_A, 1.) +
741 ROOT::Math::normal_cdf_c( sqrtqmu - sqrtqmu_A, 1.);
742
743 }
744
745 if (useQTilde ) {
746 if (fOneSided) {
747 // for bounded one-sided (q_mu_tilde: equations 64,65)
748 if ( qmu > qmu_A && (qmu_A > 0 || qmu > tol) ) { // to avoid case 0/0
749 if (verbose > 2) oocoutI(nullptr,Eval) << "Using qmu_tilde (qmu is greater than qmu_A)" << endl;
750 pnull = ROOT::Math::normal_cdf_c( (qmu + qmu_A)/(2 * sqrtqmu_A), 1.);
751 palt = ROOT::Math::normal_cdf_c( (qmu - qmu_A)/(2 * sqrtqmu_A), 1.);
752 }
753 }
754 else {
755 // for 2 sided bounded test statistic (N.B there is no one sided discovery qtilde)
756 // t_mu_tilde: equations 43,44 in asymptotic paper
757 if ( qmu > qmu_A && (qmu_A > 0 || qmu > tol) ) {
758 if (verbose > 2) oocoutI(nullptr,Eval) << "Using tmu_tilde (qmu is greater than qmu_A)" << endl;
759 pnull = ROOT::Math::normal_cdf_c(sqrtqmu,1.) +
760 ROOT::Math::normal_cdf_c( (qmu + qmu_A)/(2 * sqrtqmu_A), 1.);
761 palt = ROOT::Math::normal_cdf_c( sqrtqmu_A + sqrtqmu, 1.) +
762 ROOT::Math::normal_cdf_c( (qmu - qmu_A)/(2 * sqrtqmu_A), 1.);
763 }
764 }
765 }
766
767
768
769 // create an HypoTest result but where the sampling distributions are set to zero
770 string resultname = "HypoTestAsymptotic_result";
771 HypoTestResult* res = new HypoTestResult(resultname.c_str(), pnull, palt);
772
773 if (verbose > 0) {
774 //std::cout
775 oocoutP(nullptr, Eval) << "poi = " << muTest->getVal() << " qmu = " << qmu << " qmu_A = " << qmu_A
776 << " sigma = " << muTest->getVal() / sqrtqmu_A << " CLsplusb = " << pnull
777 << " CLb = " << palt << " CLs = " << res->CLs() << std::endl;
778 }
779
780 return res;
781
782}
783
785 PaltFunction( double offset, double pval, int icase) :
786 fOffset(offset), fPval(pval), fCase(icase) {}
787 double operator() (double x) const {
788 return ROOT::Math::normal_cdf_c(x + fOffset) + ROOT::Math::normal_cdf_c(fCase*(x - fOffset)) - fPval;
789 }
790 double fOffset;
791 double fPval;
792 int fCase;
793};
794
795////////////////////////////////////////////////////////////////////////////////
796/// function given the null and the alt p value - return the expected one given the N - sigma value
797
798double AsymptoticCalculator::GetExpectedPValues(double pnull, double palt, double nsigma, bool useCls, bool oneSided ) {
799 if (oneSided) {
800 double sqrtqmu = ROOT::Math::normal_quantile_c( pnull,1.);
801 double sqrtqmu_A = ROOT::Math::normal_quantile( palt,1.) + sqrtqmu;
802 double clsplusb = ROOT::Math::normal_cdf_c( sqrtqmu_A - nsigma, 1.);
803 if (!useCls) return clsplusb;
804 double clb = ROOT::Math::normal_cdf( nsigma, 1.);
805 return (clb == 0) ? -1 : clsplusb / clb;
806 }
807
808 // case of 2 sided test statistic
809 // need to compute numerically
810 double sqrttmu = ROOT::Math::normal_quantile_c( 0.5*pnull,1.);
811 if (sqrttmu == 0) {
812 // here cannot invert the function - skip the point
813 return -1;
814 }
815 // invert formula for palt to get sqrttmu_A
816 PaltFunction f( sqrttmu, palt, -1);
819 brf.SetFunction( wf, 0, 20);
820 bool ret = brf.Solve();
821 if (!ret) {
822 oocoutE(nullptr,Eval) << "Error finding expected p-values - return -1" << std::endl;
823 return -1;
824 }
825 double sqrttmu_A = brf.Root();
826
827 // now invert for expected value
828 PaltFunction f2( sqrttmu_A, ROOT::Math::normal_cdf( nsigma, 1.), 1);
830 brf.SetFunction(wf2,0,20);
831 ret = brf.Solve();
832 if (!ret) {
833 oocoutE(nullptr,Eval) << "Error finding expected p-values - return -1" << std::endl;
834 return -1;
835 }
836 return 2*ROOT::Math::normal_cdf_c( brf.Root(),1.);
837}
838
839// void GetExpectedLimit(double nsigma, double alpha, double &clsblimit, double &clslimit) {
840// // get expected limit
841// double
842// }
843
844////////////////////////////////////////////////////////////////////////////////
845/// fill bins by looping recursively on observables
846
847void AsymptoticCalculator::FillBins(const RooAbsPdf & pdf, const RooArgList &obs, RooAbsData & data, int &index, double &binVolume, int &ibin) {
848
849 bool debug = (fgPrintLevel >= 2);
850
851 RooRealVar * v = dynamic_cast<RooRealVar*>(&(obs[index]) );
852 if (!v) return;
853
854 RooArgSet obstmp(obs);
855 double expectedEvents = pdf.expectedEvents(obstmp);
856 // if (debug) {
857 // std::cout << "expected events = " << expectedEvents << std::endl;
858 // }
859
860 if (debug) cout << "looping on observable " << v->GetName() << endl;
861 for (int i = 0; i < v->getBins(); ++i) {
862 v->setBin(i);
863 if (index < int(obs.size()) -1) {
864 index++; // increase index
865 double prevBinVolume = binVolume;
866 binVolume *= v->getBinWidth(i); // increase bin volume
867 FillBins(pdf, obs, data, index, binVolume, ibin);
868 index--; // decrease index
869 binVolume = prevBinVolume; // decrease also bin volume
870 }
871 else {
872
873 // this is now a new bin - compute the pdf in this bin
874 double totBinVolume = binVolume * v->getBinWidth(i);
875 double fval = pdf.getVal(&obstmp)*totBinVolume;
876
877 //if (debug) std::cout << "pdf value in the bin " << fval << " bin volume = " << totBinVolume << " " << fval*expectedEvents << std::endl;
878 if (fval*expectedEvents <= 0)
879 {
880 if (fval*expectedEvents < 0) {
881 oocoutW(nullptr,InputArguments)
882 << "AsymptoticCalculator::" << __func__
883 << "(): Detected a bin with negative expected events! Please check your inputs." << endl;
884 }
885 else {
886 oocoutW(nullptr,InputArguments)
887 << "AsymptoticCalculator::" << __func__
888 << "(): Detected a bin with zero expected events- skip it" << endl;
889 }
890 }
891 // have a cut off for overflows ??
892 else
893 data.add(obs, fval*expectedEvents);
894
895 if (debug) {
896 cout << "bin " << ibin << "\t";
897 for (std::size_t j=0; j < obs.size(); ++j) { cout << " " << (static_cast<RooRealVar&>( obs[j])).getVal(); }
898 cout << " w = " << fval*expectedEvents;
899 cout << endl;
900 }
901 // RooArgSet xxx(obs);
902 // h3->Fill(((RooRealVar&) obs[0]).getVal(), ((RooRealVar&) obs[1]).getVal(), ((RooRealVar&) obs[2]).getVal() ,
903 // pdf->getVal(&xxx) );
904 ibin++;
905 }
906 }
907 //reset bin values
908 if (debug)
909 cout << "ending loop on .. " << v->GetName() << endl;
910
911 v->setBin(0);
912
913}
914
915////////////////////////////////////////////////////////////////////////////////
916/// Inpspect a product pdf to find all the Poisson or Gaussian parts to set the observed
917/// values to expected ones.
918
920{
921 bool ret = true;
922 for (auto *a : prod.pdfList()) {
923 if (!a->dependsOn(obs)) continue;
924 RooPoisson *pois = nullptr;
925 RooGaussian * gauss = nullptr;
926 if ((pois = dynamic_cast<RooPoisson *>(a)) != nullptr) {
927 ret &= SetObsToExpected(*pois, obs);
928 pois->setNoRounding(true); //needed since expected value is not an integer
929 } else if ((gauss = dynamic_cast<RooGaussian *>(a)) != nullptr) {
930 ret &= SetObsToExpected(*gauss, obs);
931 } else {
932 // should try to add also lognormal case ?
933 RooProdPdf *subprod = dynamic_cast<RooProdPdf *>(a);
934 if (subprod) {
935 ret &= SetObsToExpected(*subprod, obs);
936 } else {
937 oocoutE(nullptr, InputArguments)
938 << "Illegal term in counting model: "
939 << "the PDF " << a->GetName() << " depends on the observables, but is not a Poisson, Gaussian or Product"
940 << endl;
941 return false;
942 }
943 }
944 }
945
946 return ret;
947}
948
949////////////////////////////////////////////////////////////////////////////////
950/// set observed value to the expected one
951/// works for Gaussian, Poisson or LogNormal
952/// assumes mean parameter value is the argument not constant and not depending on observables
953/// (if more than two arguments are not constant will use first one but print a warning !)
954/// need to iterate on the components of the Poisson to get n and nu (nu can be a RooAbsReal)
955/// (code from G. Petrucciani and extended by L.M.)
956
958{
959 RooRealVar *myobs = nullptr;
960 RooAbsReal *myexp = nullptr;
961 const char * pdfName = pdf.ClassName();
962 for (RooAbsArg *a : pdf.servers()) {
963 if (obs.contains(*a)) {
964 if (myobs != nullptr) {
965 oocoutF(nullptr,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : Has two observables ?? " << endl;
966 return false;
967 }
968 myobs = dynamic_cast<RooRealVar *>(a);
969 if (myobs == nullptr) {
970 oocoutF(nullptr,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : Observable is not a RooRealVar??" << endl;
971 return false;
972 }
973 } else {
974 if (!a->isConstant() ) {
975 if (myexp != nullptr) {
976 oocoutE(nullptr,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : Has two non-const arguments " << endl;
977 return false;
978 }
979 myexp = dynamic_cast<RooAbsReal *>(a);
980 if (myexp == nullptr) {
981 oocoutF(nullptr,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : Expected is not a RooAbsReal??" << endl;
982 return false;
983 }
984 }
985 }
986 }
987 if (myobs == nullptr) {
988 oocoutF(nullptr,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : No observable?" << endl;
989 return false;
990 }
991 if (myexp == nullptr) {
992 oocoutF(nullptr,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : No observable?" << endl;
993 return false;
994 }
995
996 myobs->setVal(myexp->getVal());
997
998 if (fgPrintLevel > 2) {
999 std::cout << "SetObsToExpected : setting " << myobs->GetName() << " to expected value " << myexp->getVal() << " of " << myexp->GetName() << std::endl;
1000 }
1001
1002 return true;
1003}
1004
1005////////////////////////////////////////////////////////////////////////////////
1006/// Generate counting Asimov data for the case when the pdf cannot be extended.
1007/// This function assumes that the pdf is a RooPoisson or can be decomposed in a product of RooPoisson,
1008/// or is a RooGaussian. Otherwise, we cannot know how to make the Asimov data sets.
1009
1011 RooArgSet obs(observables);
1012 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
1013 RooPoisson *pois = nullptr;
1014 RooGaussian *gauss = nullptr;
1015
1016 if (fgPrintLevel > 1)
1017 std::cout << "generate counting Asimov data for pdf of type " << pdf.ClassName() << std::endl;
1018
1019 bool r = false;
1020 if (prod != nullptr) {
1021 r = SetObsToExpected(*prod, observables);
1022 } else if ((pois = dynamic_cast<RooPoisson *>(&pdf)) != nullptr) {
1023 r = SetObsToExpected(*pois, observables);
1024 // we need in this case to set Poisson to real values
1025 pois->setNoRounding(true);
1026 } else if ((gauss = dynamic_cast<RooGaussian *>(&pdf)) != nullptr) {
1027 r = SetObsToExpected(*gauss, observables);
1028 } else {
1029 oocoutE(nullptr,InputArguments) << "A counting model pdf must be either a RooProdPdf or a RooPoisson or a RooGaussian" << endl;
1030 }
1031 if (!r) return nullptr;
1032 int icat = 0;
1033 if (channelCat) {
1034 icat = channelCat->getCurrentIndex();
1035 }
1036
1037 RooDataSet *ret = new RooDataSet(std::string("CountingAsimovData") + std::to_string(icat),
1038 std::string("CountingAsimovData") + std::to_string(icat), obs);
1039 ret->add(obs);
1040 return ret;
1041}
1042
1043////////////////////////////////////////////////////////////////////////////////
1044/// Compute the asimov data set for an observable of a pdf.
1045/// It generates binned data following the binning of the observables.
1046// TODO: (possibility to change number of bins)
1047// TODO: implement integration over bin content
1048
1049RooAbsData * AsymptoticCalculator::GenerateAsimovDataSinglePdf(const RooAbsPdf & pdf, const RooArgSet & allobs, const RooRealVar & weightVar, RooCategory * channelCat) {
1050
1051 int printLevel = fgPrintLevel;
1052
1053 // Get observables defined by the pdf associated with this state
1054 std::unique_ptr<RooArgSet> obs(pdf.getObservables(allobs) );
1055
1056
1057 // if pdf cannot be extended assume is then a counting experiment
1058 if (!pdf.canBeExtended() ) return GenerateCountingAsimovData(const_cast<RooAbsPdf&>(pdf), *obs, weightVar, channelCat);
1059
1060 RooArgSet obsAndWeight(*obs);
1061 obsAndWeight.add(weightVar);
1062
1063 RooDataSet* asimovData = nullptr;
1064 if (channelCat) {
1065 int icat = channelCat->getCurrentIndex();
1066 asimovData = new RooDataSet(std::string("AsimovData") + std::to_string(icat),
1067 std::string("combAsimovData") + std::to_string(icat),
1068 RooArgSet(obsAndWeight,*channelCat),RooFit::WeightVar(weightVar));
1069 }
1070 else
1071 asimovData = new RooDataSet("AsimovData","AsimovData",RooArgSet(obsAndWeight),RooFit::WeightVar(weightVar));
1072
1073 // This works only for 1D observables
1074 //RooRealVar* thisObs = ((RooRealVar*)obstmp->first());
1075
1076 RooArgList obsList(*obs);
1077
1078 // loop on observables and on the bins
1079 if (printLevel >= 2) {
1080 cout << "Generating Asimov data for pdf " << pdf.GetName() << endl;
1081 cout << "list of observables " << endl;
1082 obsList.Print();
1083 }
1084
1085 int obsIndex = 0;
1086 double binVolume = 1;
1087 int nbins = 0;
1088 FillBins(pdf, obsList, *asimovData, obsIndex, binVolume, nbins);
1089 if (printLevel >= 2)
1090 cout << "filled from " << pdf.GetName() << " " << nbins << " nbins " << " volume is " << binVolume << endl;
1091
1092 // for (int iobs = 0; iobs < obsList.size(); ++iobs) {
1093 // RooRealVar * thisObs = dynamic_cast<RooRealVar*> &obsList[i];
1094 // if (thisObs == 0) continue;
1095 // // loop on the bin contents
1096 // for(int ibin=0; ibin<thisObs->numBins(); ++ibin){
1097 // thisObs->setBin(ibin);
1098
1099 // thisNorm=pdftmp->getVal(obstmp)*thisObs->getBinWidth(jj);
1100 // if (thisNorm*expectedEvents <= 0)
1101 // {
1102 // cout << "WARNING::Detected bin with zero expected events! Please check your inputs." << endl;
1103 // }
1104 // // have a cut off for overflows ??
1105 // obsDataUnbinned->add(*mc->GetObservables(), thisNorm*expectedEvents);
1106 // }
1107
1108 if (printLevel >= 1)
1109 {
1110 asimovData->Print();
1111 //cout <<"sum entries "<< asimovData->sumEntries()<<endl;
1112 }
1113 if( TMath::IsNaN(asimovData->sumEntries()) ){
1114 cout << "sum entries is nan"<<endl;
1115 assert(0);
1116 delete asimovData;
1117 asimovData = nullptr;
1118 }
1119
1120 return asimovData;
1121
1122}
1123
1124////////////////////////////////////////////////////////////////////////////////
1125/// generate the asimov data for the observables (not the global ones)
1126/// need to deal with the case of a sim pdf
1127
1129
1130 int printLevel = fgPrintLevel;
1131
1132 unique_ptr<RooRealVar> weightVar (new RooRealVar("binWeightAsimov", "binWeightAsimov", 1, 0, 1.E30 ));
1133
1134 if (printLevel > 1) cout <<" Generate Asimov data for observables"<<endl;
1135 //RooDataSet* simData=nullptr;
1136 const RooSimultaneous* simPdf = dynamic_cast<const RooSimultaneous*>(&pdf);
1137 if (!simPdf) {
1138 // generate data for non sim pdf
1139 return GenerateAsimovDataSinglePdf( pdf, observables, *weightVar, nullptr);
1140 }
1141
1142 std::map<std::string, RooDataSet*> asimovDataMap;
1143
1144 //look at category of simpdf
1145 RooCategory& channelCat = const_cast<RooCategory&>(dynamic_cast<const RooCategory&>(simPdf->indexCat()));
1146 int nrIndices = channelCat.numTypes();
1147 if( nrIndices == 0 ) {
1148 oocoutW(nullptr,Generation) << "Simultaneous pdf does not contain any categories." << endl;
1149 }
1150 for (int i=0;i<nrIndices;i++){
1151 channelCat.setIndex(i);
1152 //iFrame++;
1153 // Get pdf associated with state from simpdf
1154 RooAbsPdf* pdftmp = simPdf->getPdf(channelCat.getCurrentLabel()) ;
1155 assert(pdftmp != nullptr);
1156
1157 if (printLevel > 1)
1158 {
1159 cout << "on type " << channelCat.getCurrentLabel() << " " << channelCat.getCurrentIndex() << endl;
1160 }
1161
1162 RooAbsData * dataSinglePdf = GenerateAsimovDataSinglePdf( *pdftmp, observables, *weightVar, &channelCat);
1163 //((RooRealVar*)obstmp->first())->Print();
1164 //cout << "expected events " << pdftmp->expectedEvents(*obstmp) << endl;
1165 if (!dataSinglePdf) {
1166 oocoutE(nullptr,Generation) << "Error generating an Asimov data set for pdf " << pdftmp->GetName() << endl;
1167 return nullptr;
1168 }
1169
1170 if (asimovDataMap.count(string(channelCat.getCurrentLabel())) != 0) {
1171 oocoutE(nullptr,Generation) << "AsymptoticCalculator::GenerateAsimovData(): The PDF for " << channelCat.getCurrentLabel()
1172 << " was already defined. It will be overridden. The faulty category definitions follow:" << endl;
1173 channelCat.Print("V");
1174 }
1175
1176 asimovDataMap[string(channelCat.getCurrentLabel())] = static_cast<RooDataSet*>(dataSinglePdf);
1177
1178 if (printLevel > 1)
1179 {
1180 cout << "channel: " << channelCat.getCurrentLabel() << ", data: ";
1181 dataSinglePdf->Print();
1182 cout << endl;
1183 }
1184 }
1185
1186 RooArgSet obsAndWeight(observables);
1187 obsAndWeight.add(*weightVar);
1188
1189
1190 RooDataSet* asimovData = new RooDataSet("asimovDataFullModel","asimovDataFullModel",RooArgSet(obsAndWeight,channelCat),
1191 RooFit::Index(channelCat),RooFit::Import(asimovDataMap),RooFit::WeightVar(*weightVar));
1192
1193 for (auto &element : asimovDataMap) {
1194 delete element.second;
1195 }
1196
1197 return asimovData;
1198
1199}
1200
1201////////////////////////////////////////////////////////////////////////////////
1202/// Make the Asimov data from the ModelConfig and list of poi
1203/// \param realData Real data
1204/// \param model Model config defining the pdf and the parameters
1205/// \param paramValues The snapshot of POI and parameters used for finding the best nuisance parameter values (conditioned at these values)
1206/// \param[out] asimovGlobObs Global observables set to values satisfying the constraints
1207/// \param genPoiValues Optional. A different set of POI values used for generating. By default the same POI are used for generating and for finding the nuisance parameters
1208/// given an observed data set, a model and a snapshot of the poi.
1209/// \return The asimov data set. The user takes ownership.
1210///
1211
1212RooAbsData * AsymptoticCalculator::MakeAsimovData(RooAbsData & realData, const ModelConfig & model, const RooArgSet & paramValues, RooArgSet & asimovGlobObs, const RooArgSet * genPoiValues ) {
1213
1214 int verbose = fgPrintLevel;
1215
1216
1217 RooArgSet poi(*model.GetParametersOfInterest());
1218 poi.assign(paramValues);
1219
1220 // set poi constant for conditional MLE
1221 // need to fit nuisance parameters at their conditional MLE value
1222 RooArgSet paramsSetConstant;
1223 for (auto *tmpPar : static_range_cast<RooRealVar *>(poi)) {
1224 tmpPar->setConstant();
1225 if (verbose>0)
1226 std::cout << "MakeAsimov: Setting poi " << tmpPar->GetName() << " to a constant value = " << tmpPar->getVal() << std::endl;
1227 paramsSetConstant.add(*tmpPar);
1228 }
1229
1230 // find conditional value of the nuisance parameters
1231 bool hasFloatParams = false;
1232 RooArgSet constrainParams;
1233 if (model.GetNuisanceParameters()) {
1234 constrainParams.add(*model.GetNuisanceParameters());
1235 RooStats::RemoveConstantParameters(&constrainParams);
1236 if (!constrainParams.empty()) hasFloatParams = true;
1237
1238 } else {
1239 // Do we have free parameters anyway that need fitting?
1240 std::unique_ptr<RooArgSet> params(model.GetPdf()->getParameters(realData));
1241 for (auto const *rrv : dynamic_range_cast<RooRealVar *>(*params)) {
1242 if ( rrv != nullptr && rrv->isConstant() == false ) { hasFloatParams = true; break; }
1243 }
1244 }
1245 if (hasFloatParams) {
1246 // models need to be fitted to find best nuisance parameter values
1247
1248 TStopwatch tw2; tw2.Start();
1250 if (verbose>0) {
1251 std::cout << "MakeAsimov: doing a conditional fit for finding best nuisance values " << std::endl;
1252 minimPrintLevel = verbose;
1253 if (verbose>1) {
1254 std::cout << "POI values:\n"; poi.Print("v");
1255 if (verbose > 2) {
1256 std::cout << "Nuis param values:\n";
1257 constrainParams.Print("v");
1258 }
1259 }
1260 }
1263
1264 RooArgSet conditionalObs;
1265 if (model.GetConditionalObservables()) conditionalObs.add(*model.GetConditionalObservables());
1266 RooArgSet globalObs;
1267 if (model.GetGlobalObservables()) globalObs.add(*model.GetGlobalObservables());
1268
1269 std::string minimizerAlgo = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo();
1270 std::vector<RooCmdArg> args;
1271 args.push_back(RooFit::Minimizer("",minimizerAlgo.c_str())); // empty mimimizer type to select default
1273 args.push_back(RooFit::PrintLevel(minimPrintLevel-1));
1274 args.push_back(RooFit::Hesse(false));
1275 args.push_back(RooFit::Constrain(constrainParams));
1276 args.push_back(RooFit::GlobalObservables(globalObs));
1277 args.push_back(RooFit::ConditionalObservables(conditionalObs));
1278 args.push_back(RooFit::Offset(GetGlobalRooStatsConfig().useLikelihoodOffset));
1279 args.push_back(RooFit::EvalErrorWall(GetGlobalRooStatsConfig().useEvalErrorWall));
1280
1281 RooLinkedList argList;
1282 for (auto& arg : args) {
1283 argList.Add(&arg);
1284 }
1285 model.GetPdf()->fitTo(realData, argList);
1286 if (verbose>0) { std::cout << "fit time "; tw2.Print();}
1287 if (verbose > 1) {
1288 // after the fit the nuisance parameters will have their best fit value
1289 if (model.GetNuisanceParameters() ) {
1290 std::cout << "Nuisance parameters after fit for asimov dataset: " << std::endl;
1291 model.GetNuisanceParameters()->Print("V");
1292 }
1293 }
1294
1295 if (verbose < 2) RooMsgService::instance().setGlobalKillBelow(msglevel);
1296
1297 }
1298
1299 // restore the parameters which were set constant
1300 SetAllConstant(paramsSetConstant, false);
1301
1302
1303 std::unique_ptr<RooArgSet> allParams{model.GetPdf()->getParameters(realData)};
1305
1306 // if a RooArgSet of poi is passed , different poi will be used for generating the Asimov data set
1307 if (genPoiValues) {
1308 allParams->assign(*genPoiValues);
1309 }
1310
1311 // now do the actual generation of the AsimovData Set
1312 // no need to pass parameters values since we have set them before
1313 return MakeAsimovData(model, *allParams, asimovGlobObs);
1314}
1315
1316////////////////////////////////////////////////////////////////////////////////
1317/// \param model ModelConfig that contains the model pdf and the model parameters
1318/// \param allParamValues The parameters of the model will be set to the values given in this set
1319/// \param[out] asimovGlobObs Global observables set to values satisfying the constraints
1320/// \return Asimov data set. The user takes ownership.
1321///
1322/// The parameter values (including the nuisance parameter) can result from a fit to data or be at the nominal values.
1323///
1324
1325RooAbsData * AsymptoticCalculator::MakeAsimovData(const ModelConfig & model, const RooArgSet & allParamValues, RooArgSet & asimovGlobObs) {
1326
1327 int verbose = fgPrintLevel;
1328
1329 TStopwatch tw;
1330 tw.Start();
1331
1332 // set the parameter values (do I need the poi to be constant ? )
1333 // the nuisance parameter values could be set at their fitted value (the MLE)
1334 if (!allParamValues.empty()) {
1335 std::unique_ptr<RooArgSet> allVars{model.GetPdf()->getVariables()};
1336 allVars->assign(allParamValues);
1337 }
1338
1339
1340 // generate the Asimov data set for the observables
1341 RooAbsData * asimov = GenerateAsimovData(*model.GetPdf() , *model.GetObservables() );
1342
1343 if (verbose>0) {
1344 std::cout << "Generated Asimov data for observables "; (model.GetObservables() )->Print();
1345 if (verbose > 1) {
1346 if (asimov->numEntries() == 1 ) {
1347 std::cout << "--- Asimov data values \n";
1348 asimov->get()->Print("v");
1349 }
1350 else {
1351 std::cout << "--- Asimov data numEntries = " << asimov->numEntries() << " sumOfEntries = " << asimov->sumEntries() << std::endl;
1352 }
1353 std::cout << "\ttime for generating : "; tw.Print();
1354 }
1355 }
1356
1357
1358 // Now need to have in ASIMOV the data sets also the global observables
1359 // Their values must be the one satisfying the constraint.
1360 // to do it make a nuisance pdf with all product of constraints and then
1361 // assign to each constraint a glob observable value = to the current fitted nuisance parameter value
1362 // IN general one should solve in general the system of equations f( gobs| nuispar ) = 0 where f are the
1363 // derivatives of the constraint with respect the nuisance parameter and they are evaluated at the best fit nuisance
1364 // parameter points
1365 // As simple solution assume that constrain has a direct dependence on the nuisance parameter, i.e.
1366 // Constraint (gobs, func( nuispar) ) and the condition is satisfied for
1367 // gobs = func( nuispar) where nunispar is at the MLE value
1368
1369
1370 if (model.GetGlobalObservables() && !model.GetGlobalObservables()->empty()) {
1371
1372 if (verbose>1) {
1373 std::cout << "Generating Asimov data for global observables " << std::endl;
1374 }
1375
1376 RooArgSet gobs(*model.GetGlobalObservables());
1377
1378 // snapshot data global observables
1379 RooArgSet snapGlobalObsData;
1380 SetAllConstant(gobs, true);
1381 gobs.snapshot(snapGlobalObsData);
1382
1383
1384 RooArgSet nuis;
1385 if (model.GetNuisanceParameters()) nuis.add(*model.GetNuisanceParameters());
1386 if (nuis.empty()) {
1387 oocoutW(nullptr,Generation) << "AsymptoticCalculator::MakeAsimovData: model does not have nuisance parameters but has global observables"
1388 << " set global observables to model values " << endl;
1389 asimovGlobObs.assign(gobs);
1390 return asimov;
1391 }
1392
1393 // part 1: create the nuisance pdf
1394 std::unique_ptr<RooAbsPdf> nuispdf(RooStats::MakeNuisancePdf(model,"TempNuisPdf") );
1395 if (nuispdf == nullptr) {
1396 oocoutF(nullptr, Generation) << "AsymptoticCalculator::MakeAsimovData: model has nuisance parameters and "
1397 "global obs but no nuisance pdf "
1398 << std::endl;
1399 }
1400 // unfold the nuisance pdf if it is a prod pdf
1401 RooArgList pdfList;
1402 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(nuispdf.get());
1403 if (prod ) {
1404 pdfList.add(prod->pdfList());
1405 } else {
1406 // nothing to unfold - just use the pdf
1407 pdfList.add(*nuispdf);
1408 }
1409
1410 for (auto *cterm : static_range_cast<RooAbsPdf *>(pdfList)) {
1411 assert(dynamic_cast<RooAbsPdf *>(static_cast<RooAbsArg *>(cterm)) &&
1412 "AsimovUtils: a factor of the nuisance pdf is not a Pdf!");
1413
1414 if (!cterm->dependsOn(nuis)) continue; // dummy constraints
1415 // skip also the case of uniform components
1416 if (typeid(*cterm) == typeid(RooUniform)) continue;
1417
1418 std::unique_ptr<RooArgSet> cpars(cterm->getParameters(&gobs));
1419 std::unique_ptr<RooArgSet> cgobs(cterm->getObservables(&gobs));
1420 if (cgobs->size() > 1) {
1421 oocoutE(nullptr,Generation) << "AsymptoticCalculator::MakeAsimovData: constraint term " << cterm->GetName()
1422 << " has multiple global observables -cannot generate - skip it" << std::endl;
1423 continue;
1424 }
1425 else if (cgobs->empty()) {
1426 oocoutW(nullptr, Generation)
1427 << "AsymptoticCalculator::MakeAsimovData: constraint term " << cterm->GetName()
1428 << " has no global observables - skip it" << std::endl;
1429 continue;
1430 }
1431 // the variable representing the global observable
1432 RooRealVar &rrv = dynamic_cast<RooRealVar &>(*cgobs->first());
1433
1434 // remove the constant parameters in cpars
1436 if (cpars->size() != 1) {
1437 oocoutE(nullptr, Generation)
1438 << "AsymptoticCalculator::MakeAsimovData:constraint term "
1439 << cterm->GetName() << " has multiple floating params - cannot generate - skip it " << std::endl;
1440 continue;
1441 }
1442
1443 bool foundServer = false;
1444 // note : this will work only for this type of constraints
1445 // expressed as RooPoisson, RooGaussian, RooLognormal, RooGamma
1446 TClass * cClass = cterm->IsA();
1447 if (verbose > 2) std::cout << "Constraint " << cterm->GetName() << " of type " << cClass->GetName() << std::endl;
1448 if ( cClass != RooGaussian::Class() && cClass != RooPoisson::Class() &&
1449 cClass != RooGamma::Class() && cClass != RooLognormal::Class() &&
1450 cClass != RooBifurGauss::Class() ) {
1451 TString className = (cClass) ? cClass->GetName() : "undefined";
1452 oocoutW(nullptr, Generation)
1453 << "AsymptoticCalculator::MakeAsimovData:constraint term "
1454 << cterm->GetName() << " of type " << className
1455 << " is a non-supported type - result might be not correct " << std::endl;
1456 }
1457
1458 // in case of a Poisson constraint make sure the rounding is not set
1459 if (cClass == RooPoisson::Class() ) {
1460 RooPoisson * pois = static_cast<RooPoisson*>(cterm);
1461 assert(dynamic_cast<RooPoisson *>(cterm));
1462 pois->setNoRounding(true);
1463 }
1464
1465 // look at server of the constraint term and check if the global observable is part of the server
1466 RooAbsArg * arg = cterm->findServer(rrv);
1467 if (!arg) {
1468 // special case is for the Gamma where one might define the global observable n and you have a Gamma(b, n+1, ...._
1469 // in this case n+1 is the server and we don;t have a direct dependency, but we want to set n to the b value
1470 // so in case of the Gamma ignore this test
1471 if ( cClass != RooGamma::Class() ) {
1472 oocoutE(nullptr, Generation)
1473 << "AsymptoticCalculator::MakeAsimovData:constraint term "
1474 << cterm->GetName() << " has no direct dependence on global observable- cannot generate it " << std::endl;
1475 continue;
1476 }
1477 }
1478
1479 // loop on the server of the constraint term
1480 // need to treat the Gamma as a special case
1481 // the mode of the Gamma is (k-1)*theta where theta is the inverse of the rate parameter.
1482 // we assume that the global observable is defined as ngobs = k-1 and the theta parameter has the name theta otherwise we use other procedure which might be wrong
1483 RooAbsReal * thetaGamma = nullptr;
1484 if ( cClass == RooGamma::Class() ) {
1485 for (RooAbsArg *a2 : cterm->servers()) {
1486 if (TString(a2->GetName()).Contains("theta") ) {
1487 thetaGamma = dynamic_cast<RooAbsReal*>(a2);
1488 break;
1489 }
1490 }
1491 if (thetaGamma == nullptr) {
1492 oocoutI(nullptr, Generation)
1493 << "AsymptoticCalculator::MakeAsimovData:constraint term "
1494 << cterm->GetName() << " is a Gamma distribution and no server named theta is found. Assume that the Gamma scale is 1 " << std::endl;
1495 }
1496 else {
1497 if (verbose>2)
1498 std::cout << "Gamma constraint has a scale " << thetaGamma->GetName() << " = " << thetaGamma->getVal() << std::endl;
1499 }
1500 }
1501 for (RooAbsArg *a2 : cterm->servers()) {
1502 RooAbsReal * rrv2 = dynamic_cast<RooAbsReal *>(a2);
1503 if (verbose > 2) std::cout << "Loop on constraint server term " << a2->GetName() << std::endl;
1504 if (rrv2 && rrv2->dependsOn(nuis) ) {
1505
1506
1507 // found server depending on nuisance
1508 if (foundServer) {
1509 oocoutE(nullptr,Generation) << "AsymptoticCalculator::MakeAsimovData:constraint term "
1510 << cterm->GetName() << " constraint term has more server depending on nuisance- cannot generate it " <<
1511 std::endl;
1512 foundServer = false;
1513 break;
1514 }
1515 if (thetaGamma && thetaGamma->getVal() > 0) {
1516 rrv.setVal( rrv2->getVal() / thetaGamma->getVal() );
1517 } else {
1518 rrv.setVal(rrv2->getVal());
1519 }
1520 foundServer = true;
1521
1522 if (verbose > 2) {
1523 std::cout << "setting global observable " << rrv.GetName() << " to value " << rrv.getVal()
1524 << " which comes from " << rrv2->GetName() << std::endl;
1525 }
1526 }
1527 }
1528
1529 if (!foundServer) {
1530 oocoutE(nullptr,Generation) << "AsymptoticCalculator::MakeAsimovData - can't find nuisance for constraint term - global observables will not be set to Asimov value " << cterm->GetName() << std::endl;
1531 std::cerr << "Parameters: " << std::endl;
1532 cpars->Print("V");
1533 std::cerr << "Observables: " << std::endl;
1534 cgobs->Print("V");
1535 }
1536// rrv.setVal(match->getVal());
1537 }
1538
1539 // make a snapshot of global observables
1540 // needed this ?? (LM)
1541
1542 asimovGlobObs.removeAll();
1543 SetAllConstant(gobs, true);
1544 gobs.snapshot(asimovGlobObs);
1545
1546 // revert global observables to the data value
1547 gobs.assign(snapGlobalObsData);
1548
1549 if (verbose>0) {
1550 std::cout << "Generated Asimov data for global observables ";
1551 if (verbose == 1) gobs.Print();
1552 }
1553
1554 if (verbose > 1) {
1555 std::cout << "\nGlobal observables for data: " << std::endl;
1556 gobs.Print("V");
1557 std::cout << "\nGlobal observables for asimov: " << std::endl;
1558 asimovGlobObs.Print("V");
1559 }
1560
1561
1562 }
1563
1564 return asimov;
1565
1566}
#define f(i)
Definition RSha256.hxx:104
#define a(i)
Definition RSha256.hxx:99
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define oocoutW(o, a)
#define oocoutE(o, a)
#define oocoutF(o, a)
#define oocoutI(o, a)
#define oocoutP(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h offset
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void 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 char Point_t Rectangle_t WindowAttributes_t index
void Print(GNN_Data &d, std::string txt="")
TRObject operator()(const T1 &t1) const
Class for finding the root of a one dimensional function using the Brent algorithm.
bool SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup) override
Sets the function for the rest of the algorithms.
bool Solve(int maxIter=100, double absTol=1E-8, double relTol=1E-10) override
Returns the X value corresponding to the function value fy for (xmin<x<xmax).
double Root() const override
Returns root value.
static const std::string & DefaultMinimizerAlgo()
Template class to wrap any C++ callable object which takes one argument i.e.
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.
bool isConstant() const
Check if the "Constant" attribute is set.
Definition RooAbsArg.h:361
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expression tree)
RooAbsArg * findServer(const char *name) const
Return server of this with name name. Returns nullptr if not found.
Definition RooAbsArg.h:210
virtual const char * getCurrentLabel() const
Return label string of current state.
Int_t numTypes(const char *=nullptr) const
Return number of types defined (in range named rangeName if rangeName!=nullptr)
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
bool contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
const char * GetName() const override
Returns name of object.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
Storage_t::size_type size() const
RooAbsArg * first() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
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
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
Definition RooAbsData.h:225
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual double expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
bool canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition RooAbsPdf.h:219
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:157
virtual Int_t getBins(const char *name=nullptr) const
Get number of bins of currently defined range.
void setConstant(bool value=true)
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
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
static bool hideOffset()
static void setHideOffset(bool flag)
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:178
static TClass * Class()
Object to represent discrete states.
Definition RooCategory.h:28
bool setIndex(Int_t index, bool printError=true) override
Set value by specifying the index code of the desired state.
value_type getCurrentIndex() const final
Return current index.
Definition RooCategory.h:40
static TClass * Class()
Container class to hold unbinned data.
Definition RooDataSet.h:57
double sumEntries() const override
Return effective number of entries in dataset, i.e., sum all weights.
void add(const RooArgSet &row, double weight, double weightError)
Add one ore more rows of data.
static TClass * Class()
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
static TClass * Class()
Collection class for internal use, storing a collection of RooAbsArg pointers in a doubly linked list...
virtual void Add(TObject *arg)
static TClass * Class()
Wrapper class around ROOT::Fit:Fitter that provides a seamless interface between the minimizer functi...
void optimizeConst(int flag)
If flag is true, perform constant term optimization on function being minimized.
RooFit::OwningPtr< RooFitResult > save(const char *name=nullptr, const char *title=nullptr)
Save and return a RooFitResult snapshot of current minimizer status.
void setEvalErrorWall(bool flag)
void setEps(double eps)
Change MINUIT epsilon.
void setPrintLevel(int newLevel)
Change the MINUIT internal printing level.
int minimize(const char *type, const char *alg=nullptr)
Minimise the function passed in the constructor.
void setStrategy(int istrat)
Change MINUIT strategy to istrat.
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
Poisson pdf.
Definition RooPoisson.h:19
static TClass * Class()
void setNoRounding(bool flag=true)
Switch off/on rounding of x to the nearest integer.
Definition RooPoisson.h:36
Efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:33
const RooArgList & pdfList() const
Definition RooProdPdf.h:67
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
double getError() const
Definition RooRealVar.h:58
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
RooAbsPdf * getPdf(RooStringView catName) const
Return the p.d.f associated with the given index category name.
const RooAbsCategoryLValue & indexCat() const
Hypothesis Test Calculator based on the asymptotic formulae for the profile likelihood ratio.
static double GetExpectedPValues(double pnull, double palt, double nsigma, bool usecls, bool oneSided=true)
function given the null and the alt p value - return the expected one given the N - sigma value
static double EvaluateNLL(RooStats::ModelConfig const &modelConfig, RooAbsData &data, const RooArgSet *poiSet=nullptr)
static int fgPrintLevel
control print level (0 minimal, 1 normal, 2 debug)
static bool SetObsToExpected(RooAbsPdf &pdf, const RooArgSet &obs)
set observed value to the expected one works for Gaussian, Poisson or LogNormal assumes mean paramete...
static void SetPrintLevel(int level)
set print level (static function)
RooArgSet fAsimovGlobObs
snapshot of Asimov global observables
static RooAbsData * GenerateAsimovData(const RooAbsPdf &pdf, const RooArgSet &observables)
generate the asimov data for the observables (not the global ones) need to deal with the case of a si...
int fUseQTilde
flag to indicate if using qtilde or not (-1 (default based on RooRealVar)), 0 false,...
static RooAbsData * GenerateAsimovDataSinglePdf(const RooAbsPdf &pdf, const RooArgSet &obs, const RooRealVar &weightVar, RooCategory *channelCat=nullptr)
Compute the asimov data set for an observable of a pdf.
bool fIsInitialized
! flag to check if calculator is initialized
HypoTestResult * GetHypoTest() const override
re-implement HypoTest computation using the asymptotic
bool fOneSided
for one sided PL test statistic (upper limits)
RooArgSet fBestFitParams
snapshot of all best fitted Parameter values
AsymptoticCalculator(RooAbsData &data, const ModelConfig &altModel, const ModelConfig &nullModel, bool nominalAsimov=false)
constructor for asymptotic calculator from Data set and ModelConfig
static RooAbsData * GenerateCountingAsimovData(RooAbsPdf &pdf, const RooArgSet &obs, const RooRealVar &weightVar, RooCategory *channelCat=nullptr)
Generate counting Asimov data for the case when the pdf cannot be extended.
bool fOneSidedDiscovery
for one sided PL test statistic (for discovery)
RooAbsData * fAsimovData
asimov data set
RooArgSet fBestFitPoi
snapshot of best fitted POI values
static RooAbsData * MakeAsimovData(RooAbsData &data, const ModelConfig &model, const RooArgSet &poiValues, RooArgSet &globObs, const RooArgSet *genPoiValues=nullptr)
Make Asimov data.
static void FillBins(const RooAbsPdf &pdf, const RooArgList &obs, RooAbsData &data, int &index, double &binVolume, int &ibin)
fill bins by looping recursively on observables
bool fNominalAsimov
make Asimov at nominal parameter values
bool Initialize() const
initialize the calculator by performing a global fit and make the Asimov data set
Common base class for the Hypothesis Test Calculators.
const ModelConfig * GetNullModel(void) const
const ModelConfig * GetAlternateModel(void) const
HypoTestResult is a base class for results from hypothesis tests.
virtual double CLs() const
is simply (not a method, but a quantity)
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:35
const RooArgSet * GetConditionalObservables() const
get RooArgSet for conditional observables (return nullptr if not existing)
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return nullptr if not existing)
std::unique_ptr< RooAbsReal > createNLL(RooAbsData &data, CmdArgs_t const &...cmdArgs) const
Wrapper around RooAbsPdf::createNLL(), where the pdf and some configuration options are retrieved fro...
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
const RooArgSet * GetObservables() const
get RooArgSet for observables (return nullptr if not existing)
const RooArgSet * GetSnapshot() const
get RooArgSet for parameters for a particular hypothesis (return nullptr if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
Flat p.d.f.
Definition RooUniform.h:24
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition TClass.h:81
TClass * IsA() const override
Definition TClass.h:614
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:207
Stopwatch class.
Definition TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void Print(Option_t *option="") const override
Print the real and cpu time passed between the start and stop events.
Basic string class.
Definition TString.h:139
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
RooCmdArg Index(RooCategory &icat)
RooCmdArg WeightVar(const char *name="weight", bool reinterpretAsWeight=false)
RooCmdArg Import(const char *state, TH1 &histo)
RooCmdArg Offset(std::string const &mode)
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg Minimizer(const char *type, const char *alg=nullptr)
RooCmdArg Hesse(bool flag=true)
RooCmdArg Strategy(Int_t code)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg EvalErrorWall(bool flag)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg ConditionalObservables(Args_t &&... argsOrArgSet)
Create a RooCmdArg to declare conditional observables.
double normal_cdf_c(double x, double sigma=1, double x0=0)
Complement of the cumulative distribution function of the normal (Gaussian) distribution (upper tail)...
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
Double_t x[n]
Definition legend1.C:17
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Namespace for the RooStats classes.
Definition Asimov.h:19
bool SetAllConstant(const RooAbsCollection &coll, bool constant=true)
utility function to set all variable constant in a collection (from G.
void RemoveConstantParameters(RooArgSet *set)
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
extract constraint terms from pdf
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
bool IsNLLOffset()
function returning if the flag to check if the flag to use NLLOffset is set
Bool_t IsNaN(Double_t x)
Definition TMath.h:892
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754.
Definition TMath.h:902
PaltFunction(double offset, double pval, int icase)