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