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