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"
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() ) {
257 globObs.add(*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;
304 if (condObs) conditionalObs.add(*condObs);
305 RooArgSet globalObs;
306 if (globObs) globalObs.add(*globObs);
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();
324 paramsSetConstant.add(*poiVar);
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();
340 // paramsSetConstant.add(*tmpParA);
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";
407 algorithm = "migradimproved";
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() ) {
606 globObs.add(*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
794struct PaltFunction {
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 cout << "WARNING::Detected a bin with negative expected events! Please check your inputs." << endl;
892 else
893 cout << "WARNING::Detected a bin with zero expected events- skip it" << endl;
894 }
895 // have a cut off for overflows ??
896 else
897 data.add(obs, fval*expectedEvents);
898
899 if (debug) {
900 cout << "bin " << ibin << "\t";
901 for (int j=0; j < obs.getSize(); ++j) { cout << " " << ((RooRealVar&) obs[j]).getVal(); }
902 cout << " w = " << fval*expectedEvents;
903 cout << endl;
904 }
905 // RooArgSet xxx(obs);
906 // h3->Fill(((RooRealVar&) obs[0]).getVal(), ((RooRealVar&) obs[1]).getVal(), ((RooRealVar&) obs[2]).getVal() ,
907 // pdf->getVal(&xxx) );
908 ibin++;
909 }
910 }
911 //reset bin values
912 if (debug)
913 cout << "ending loop on .. " << v->GetName() << endl;
914
915 v->setBin(0);
916
917}
918
919////////////////////////////////////////////////////////////////////////////////
920/// Inpspect a product pdf to find all the Poisson or Gaussian parts to set the observed
921/// values to expected ones.
922
924{
925 RooLinkedListIter iter(prod.pdfList().iterator());
926 bool ret = true;
927 for (RooAbsArg *a = (RooAbsArg *) iter.Next(); a != 0; a = (RooAbsArg *) iter.Next()) {
928 if (!a->dependsOn(obs)) continue;
929 RooPoisson *pois = 0;
930 RooGaussian * gaus = 0;
931 if ((pois = dynamic_cast<RooPoisson *>(a)) != 0) {
932 ret &= SetObsToExpected(*pois, obs);
933 pois->setNoRounding(true); //needed since expected value is not an integer
934 } else if ((gaus = dynamic_cast<RooGaussian *>(a)) != 0) {
935 ret &= SetObsToExpected(*gaus, obs);
936 } else {
937 // should try to add also lognormal case ?
938 RooProdPdf *subprod = dynamic_cast<RooProdPdf *>(a);
939 if (subprod)
940 ret &= SetObsToExpected(*subprod, obs);
941 else {
942 oocoutE((TObject*)0,InputArguments) << "Illegal term in counting model: "
943 << "the PDF " << a->GetName()
944 << " depends on the observables, but is not a Poisson, Gaussian or Product"
945 << endl;
946 return false;
947 }
948 }
949 }
950
951 return ret;
952}
953
954////////////////////////////////////////////////////////////////////////////////
955/// set observed value to the expected one
956/// works for Gaussian, Poisson or LogNormal
957/// assumes mean parameter value is the argument not constant and not depending on observables
958/// (if more than two arguments are not constant will use first one but print a warning !)
959/// need to iterate on the components of the Poisson to get n and nu (nu can be a RooAbsReal)
960/// (code from G. Petrucciani and extended by L.M.)
961
963{
964 RooRealVar *myobs = 0;
965 RooAbsReal *myexp = 0;
966 const char * pdfName = pdf.IsA()->GetName();
967 RooFIter iter(pdf.serverMIterator());
968 for (RooAbsArg *a = iter.next(); a != 0; a = iter.next()) {
969 if (obs.contains(*a)) {
970 if (myobs != 0) {
971 oocoutF((TObject*)0,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : Has two observables ?? " << endl;
972 return false;
973 }
974 myobs = dynamic_cast<RooRealVar *>(a);
975 if (myobs == 0) {
976 oocoutF((TObject*)0,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : Observable is not a RooRealVar??" << endl;
977 return false;
978 }
979 } else {
980 if (!a->isConstant() ) {
981 if (myexp != 0) {
982 oocoutE((TObject*)0,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : Has two non-const arguments " << endl;
983 return false;
984 }
985 myexp = dynamic_cast<RooAbsReal *>(a);
986 if (myexp == 0) {
987 oocoutF((TObject*)0,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : Expected is not a RooAbsReal??" << endl;
988 return false;
989 }
990 }
991 }
992 }
993 if (myobs == 0) {
994 oocoutF((TObject*)0,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : No observable?" << endl;
995 return false;
996 }
997 if (myexp == 0) {
998 oocoutF((TObject*)0,Generation) << "AsymptoticCalculator::SetObsExpected( " << pdfName << " ) : No observable?" << endl;
999 return false;
1000 }
1001
1002 myobs->setVal(myexp->getVal());
1003
1004 if (fgPrintLevel > 2) {
1005 std::cout << "SetObsToExpected : setting " << myobs->GetName() << " to expected value " << myexp->getVal() << " of " << myexp->GetName() << std::endl;
1006 }
1007
1008 return true;
1009}
1010
1011////////////////////////////////////////////////////////////////////////////////
1012/// Generate counting Asimov data for the case when the pdf cannot be extended.
1013/// This function assumes that the pdf is a RooPoisson or can be decomposed in a product of RooPoisson,
1014/// or is a RooGaussian. Otherwise, we cannot know how to make the Asimov data sets.
1015
1017 RooArgSet obs(observables);
1018 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
1019 RooPoisson *pois = 0;
1020 RooGaussian *gaus = 0;
1021
1022 if (fgPrintLevel > 1)
1023 std::cout << "generate counting Asimov data for pdf of type " << pdf.IsA()->GetName() << std::endl;
1024
1025 bool r = false;
1026 if (prod != 0) {
1027 r = SetObsToExpected(*prod, observables);
1028 } else if ((pois = dynamic_cast<RooPoisson *>(&pdf)) != 0) {
1029 r = SetObsToExpected(*pois, observables);
1030 // we need in this case to set Poisson to real values
1031 pois->setNoRounding(true);
1032 } else if ((gaus = dynamic_cast<RooGaussian *>(&pdf)) != 0) {
1033 r = SetObsToExpected(*gaus, observables);
1034 } else {
1035 oocoutE((TObject*)0,InputArguments) << "A counting model pdf must be either a RooProdPdf or a RooPoisson or a RooGaussian" << endl;
1036 }
1037 if (!r) return 0;
1038 int icat = 0;
1039 if (channelCat) {
1040 icat = channelCat->getCurrentIndex();
1041 }
1042
1043 RooDataSet *ret = new RooDataSet(TString::Format("CountingAsimovData%d",icat),TString::Format("CountingAsimovData%d",icat), obs);
1044 ret->add(obs);
1045 return ret;
1046}
1047
1048////////////////////////////////////////////////////////////////////////////////
1049/// Compute the asimov data set for an observable of a pdf.
1050/// It generates binned data following the binning of the observables.
1051// TODO: (possibility to change number of bins)
1052// TODO: implement integration over bin content
1053
1054RooAbsData * AsymptoticCalculator::GenerateAsimovDataSinglePdf(const RooAbsPdf & pdf, const RooArgSet & allobs, const RooRealVar & weightVar, RooCategory * channelCat) {
1055
1056 int printLevel = fgPrintLevel;
1057
1058 // Get observables defined by the pdf associated with this state
1059 std::unique_ptr<RooArgSet> obs(pdf.getObservables(allobs) );
1060
1061
1062 // if pdf cannot be extended assume is then a counting experiment
1063 if (!pdf.canBeExtended() ) return GenerateCountingAsimovData(const_cast<RooAbsPdf&>(pdf), *obs, weightVar, channelCat);
1064
1065 RooArgSet obsAndWeight(*obs);
1066 obsAndWeight.add(weightVar);
1067
1068 RooDataSet* asimovData = 0;
1069 if (channelCat) {
1070 int icat = channelCat->getCurrentIndex();
1071 asimovData = new RooDataSet(TString::Format("AsimovData%d",icat),TString::Format("combAsimovData%d",icat),
1072 RooArgSet(obsAndWeight,*channelCat),RooFit::WeightVar(weightVar));
1073 }
1074 else
1075 asimovData = new RooDataSet("AsimovData","AsimovData",RooArgSet(obsAndWeight),RooFit::WeightVar(weightVar));
1076
1077 // This works only for 1D observables
1078 //RooRealVar* thisObs = ((RooRealVar*)obstmp->first());
1079
1080 RooArgList obsList(*obs);
1081
1082 // loop on observables and on the bins
1083 if (printLevel >= 2) {
1084 cout << "Generating Asimov data for pdf " << pdf.GetName() << endl;
1085 cout << "list of observables " << endl;
1086 obsList.Print();
1087 }
1088
1089 int obsIndex = 0;
1090 double binVolume = 1;
1091 int nbins = 0;
1092 FillBins(pdf, obsList, *asimovData, obsIndex, binVolume, nbins);
1093 if (printLevel >= 2)
1094 cout << "filled from " << pdf.GetName() << " " << nbins << " nbins " << " volume is " << binVolume << endl;
1095
1096 // for (int iobs = 0; iobs < obsList.getSize(); ++iobs) {
1097 // RooRealVar * thisObs = dynamic_cast<RooRealVar*> &obsList[i];
1098 // if (thisObs == 0) continue;
1099 // // loop on the bin contents
1100 // for(int ibin=0; ibin<thisObs->numBins(); ++ibin){
1101 // thisObs->setBin(ibin);
1102
1103 // thisNorm=pdftmp->getVal(obstmp)*thisObs->getBinWidth(jj);
1104 // if (thisNorm*expectedEvents <= 0)
1105 // {
1106 // cout << "WARNING::Detected bin with zero expected events! Please check your inputs." << endl;
1107 // }
1108 // // have a cut off for overflows ??
1109 // obsDataUnbinned->add(*mc->GetObservables(), thisNorm*expectedEvents);
1110 // }
1111
1112 if (printLevel >= 1)
1113 {
1114 asimovData->Print();
1115 //cout <<"sum entries "<< asimovData->sumEntries()<<endl;
1116 }
1117 if( TMath::IsNaN(asimovData->sumEntries()) ){
1118 cout << "sum entries is nan"<<endl;
1119 assert(0);
1120 delete asimovData;
1121 asimovData = 0;
1122 }
1123
1124 return asimovData;
1125
1126}
1127
1128////////////////////////////////////////////////////////////////////////////////
1129/// generate the asimov data for the observables (not the global ones)
1130/// need to deal with the case of a sim pdf
1131
1133
1134 int printLevel = fgPrintLevel;
1135
1136 unique_ptr<RooRealVar> weightVar (new RooRealVar("binWeightAsimov", "binWeightAsimov", 1, 0, 1.E30 ));
1137
1138 if (printLevel > 1) cout <<" Generate Asimov data for observables"<<endl;
1139 //RooDataSet* simData=NULL;
1140 const RooSimultaneous* simPdf = dynamic_cast<const RooSimultaneous*>(&pdf);
1141 if (!simPdf) {
1142 // generate data for non sim pdf
1143 return GenerateAsimovDataSinglePdf( pdf, observables, *weightVar, 0);
1144 }
1145
1146 std::map<std::string, RooDataSet*> asimovDataMap;
1147
1148 //look at category of simpdf
1149 RooCategory& channelCat = const_cast<RooCategory&>(dynamic_cast<const RooCategory&>(simPdf->indexCat()));
1150 int nrIndices = channelCat.numTypes();
1151 if( nrIndices == 0 ) {
1152 oocoutW((TObject*)0,Generation) << "Simultaneous pdf does not contain any categories." << endl;
1153 }
1154 for (int i=0;i<nrIndices;i++){
1155 channelCat.setIndex(i);
1156 //iFrame++;
1157 // Get pdf associated with state from simpdf
1158 RooAbsPdf* pdftmp = simPdf->getPdf(channelCat.getCurrentLabel()) ;
1159 assert(pdftmp != 0);
1160
1161 if (printLevel > 1)
1162 {
1163 cout << "on type " << channelCat.getCurrentLabel() << " " << channelCat.getCurrentIndex() << endl;
1164 }
1165
1166 RooAbsData * dataSinglePdf = GenerateAsimovDataSinglePdf( *pdftmp, observables, *weightVar, &channelCat);
1167 //((RooRealVar*)obstmp->first())->Print();
1168 //cout << "expected events " << pdftmp->expectedEvents(*obstmp) << endl;
1169 if (!dataSinglePdf) {
1170 oocoutE((TObject*)0,Generation) << "Error generating an Asimov data set for pdf " << pdftmp->GetName() << endl;
1171 return 0;
1172 }
1173
1174 if (asimovDataMap.count(string(channelCat.getCurrentLabel())) != 0) {
1175 oocoutE((TObject*)0,Generation) << "AsymptoticCalculator::GenerateAsimovData(): The PDF for " << channelCat.getCurrentLabel()
1176 << " was already defined. It will be overridden. The faulty category definitions follow:" << endl;
1177 channelCat.Print("V");
1178 }
1179
1180 asimovDataMap[string(channelCat.getCurrentLabel())] = (RooDataSet*) dataSinglePdf;
1181
1182 if (printLevel > 1)
1183 {
1184 cout << "channel: " << channelCat.getCurrentLabel() << ", data: ";
1185 dataSinglePdf->Print();
1186 cout << endl;
1187 }
1188 }
1189
1190 RooArgSet obsAndWeight(observables);
1191 obsAndWeight.add(*weightVar);
1192
1193
1194 RooDataSet* asimovData = new RooDataSet("asimovDataFullModel","asimovDataFullModel",RooArgSet(obsAndWeight,channelCat),
1195 RooFit::Index(channelCat),RooFit::Import(asimovDataMap),RooFit::WeightVar(*weightVar));
1196
1197 for (auto element : asimovDataMap) {
1198 delete element.second;
1199 }
1200
1201 return asimovData;
1202
1203}
1204
1205////////////////////////////////////////////////////////////////////////////////
1206/// Make the Asimov data from the ModelConfig and list of poi
1207/// \param realData Real data
1208/// \param model Model config defining the pdf and the parameters
1209/// \param paramValues The snapshot of POI and parameters used for finding the best nuisance parameter values (conditioned at these values)
1210/// \param[out] asimovGlobObs Global observables set to values satisfying the constraints
1211/// \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
1212/// given an observed data set, a model and a snapshot of the poi.
1213/// \return The asimov data set. The user takes ownership.
1214///
1215
1216RooAbsData * AsymptoticCalculator::MakeAsimovData(RooAbsData & realData, const ModelConfig & model, const RooArgSet & paramValues, RooArgSet & asimovGlobObs, const RooArgSet * genPoiValues ) {
1217
1218 int verbose = fgPrintLevel;
1219
1220
1221 RooArgSet poi(*model.GetParametersOfInterest());
1222 poi = paramValues;
1223 //RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
1224 // set poi constant for conditional MLE
1225 // need to fit nuisance parameters at their conditional MLE value
1226 RooLinkedListIter it = poi.iterator();
1227 RooRealVar* tmpPar = NULL;
1228 RooArgSet paramsSetConstant;
1229 while((tmpPar = (RooRealVar*)it.Next())){
1230 tmpPar->setConstant();
1231 if (verbose>0)
1232 std::cout << "MakeAsimov: Setting poi " << tmpPar->GetName() << " to a constant value = " << tmpPar->getVal() << std::endl;
1233 paramsSetConstant.add(*tmpPar);
1234 }
1235
1236 // find conditional value of the nuisance parameters
1237 bool hasFloatParams = false;
1238 RooArgSet constrainParams;
1239 if (model.GetNuisanceParameters()) {
1240 constrainParams.add(*model.GetNuisanceParameters());
1241 RooStats::RemoveConstantParameters(&constrainParams);
1242 if (constrainParams.getSize() > 0) hasFloatParams = true;
1243
1244 } else {
1245 // Do we have free parameters anyway that need fitting?
1246 std::unique_ptr<RooArgSet> params(model.GetPdf()->getParameters(realData));
1247 RooLinkedListIter iter(params->iterator());
1248 for (RooAbsArg *a = (RooAbsArg *) iter.Next(); a != 0; a = (RooAbsArg *) iter.Next()) {
1249 RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);
1250 if ( rrv != 0 && rrv->isConstant() == false ) { hasFloatParams = true; break; }
1251 }
1252 }
1253 if (hasFloatParams) {
1254 // models need to be fitted to find best nuisance parameter values
1255
1256 TStopwatch tw2; tw2.Start();
1258 if (verbose>0) {
1259 std::cout << "MakeAsimov: doing a conditional fit for finding best nuisance values " << std::endl;
1260 minimPrintLevel = verbose;
1261 if (verbose>1) {
1262 std::cout << "POI values:\n"; poi.Print("v");
1263 if (verbose > 2) {
1264 std::cout << "Nuis param values:\n";
1265 constrainParams.Print("v");
1266 }
1267 }
1268 }
1271
1272 RooArgSet conditionalObs;
1273 if (model.GetConditionalObservables()) conditionalObs.add(*model.GetConditionalObservables());
1274 RooArgSet globalObs;
1275 if (model.GetGlobalObservables()) globalObs.add(*model.GetGlobalObservables());
1276
1277 std::string minimizerType = ROOT::Math::MinimizerOptions::DefaultMinimizerType();
1278 std::string minimizerAlgo = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo();
1279 std::vector<RooCmdArg> args;
1280 args.push_back(RooFit::Minimizer(minimizerType.c_str(),minimizerAlgo.c_str()));
1282 args.push_back(RooFit::PrintLevel(minimPrintLevel-1));
1283 args.push_back(RooFit::Hesse(false));
1284 args.push_back(RooFit::Constrain(constrainParams));
1285 args.push_back(RooFit::GlobalObservables(globalObs));
1286 args.push_back(RooFit::ConditionalObservables(conditionalObs));
1287 args.push_back(RooFit::Offset(GetGlobalRooStatsConfig().useLikelihoodOffset));
1288 args.push_back(RooFit::EvalErrorWall(GetGlobalRooStatsConfig().useEvalErrorWall));
1289
1290 RooLinkedList argList;
1291 for (auto& arg : args) {
1292 argList.Add(&arg);
1293 }
1294 model.GetPdf()->fitTo(realData, argList);
1295 if (verbose>0) { std::cout << "fit time "; tw2.Print();}
1296 if (verbose > 1) {
1297 // after the fit the nuisance parameters will have their best fit value
1298 if (model.GetNuisanceParameters() ) {
1299 std::cout << "Nuisance parameters after fit for asimov dataset: " << std::endl;
1300 model.GetNuisanceParameters()->Print("V");
1301 }
1302 }
1303
1305
1306 }
1307
1308 // restore the parameters which were set constant
1309 SetAllConstant(paramsSetConstant, false);
1310
1311
1312 RooArgSet * allParams = model.GetPdf()->getParameters(realData);
1314
1315 // if a RooArgSet of poi is passed , different poi will be used for generating the Asimov data set
1316 if (genPoiValues) *allParams = *genPoiValues;
1317
1318 // now do the actual generation of the AsimovData Set
1319 // no need to pass parameters values since we have set them before
1320 RooAbsData * asimovData = MakeAsimovData(model, *allParams, asimovGlobObs);
1321
1322 delete allParams;
1323
1324 return asimovData;
1325
1326}
1327
1328////////////////////////////////////////////////////////////////////////////////
1329/// \param model ModelConfig that contains the model pdf and the model parameters
1330/// \param allParamValues The parameters fo the model will be set to the values given in this set
1331/// \param[out] asimovGlobObs Global observables set to values satisfying the constraints
1332/// \return Asimov data set. The user takes ownership.
1333///
1334/// The parameter values (including the nuisance parameter) can result from a fit to data or be at the nominal values.
1335///
1336
1337RooAbsData * AsymptoticCalculator::MakeAsimovData(const ModelConfig & model, const RooArgSet & allParamValues, RooArgSet & asimovGlobObs) {
1338
1339 int verbose = fgPrintLevel;
1340
1341 TStopwatch tw;
1342 tw.Start();
1343
1344 // set the parameter values (do I need the poi to be constant ? )
1345 // the nuisance parameter values could be set at their fitted value (the MLE)
1346 if (allParamValues.getSize() > 0) {
1347 RooArgSet * allVars = model.GetPdf()->getVariables();
1348 *allVars = allParamValues;
1349 delete allVars;
1350 }
1351
1352
1353 // generate the Asimov data set for the observables
1354 RooAbsData * asimov = GenerateAsimovData(*model.GetPdf() , *model.GetObservables() );
1355
1356 if (verbose>0) {
1357 std::cout << "Generated Asimov data for observables "; (model.GetObservables() )->Print();
1358 if (verbose > 1) {
1359 if (asimov->numEntries() == 1 ) {
1360 std::cout << "--- Asimov data values \n";
1361 asimov->get()->Print("v");
1362 }
1363 else {
1364 std::cout << "--- Asimov data numEntries = " << asimov->numEntries() << " sumOfEntries = " << asimov->sumEntries() << std::endl;
1365 }
1366 std::cout << "\ttime for generating : "; tw.Print();
1367 }
1368 }
1369
1370
1371 // Now need to have in ASIMOV the data sets also the global observables
1372 // Their values must be the one satisfying the constraint.
1373 // to do it make a nuisance pdf with all product of constraints and then
1374 // assign to each constraint a glob observable value = to the current fitted nuisance parameter value
1375 // IN general one should solve in general the system of equations f( gobs| nuispar ) = 0 where f are the
1376 // derivatives of the constraint with respect the nuisance parameter and they are evaluated at the best fit nuisance
1377 // parameter points
1378 // As simple solution assume that constrain has a direct dependence on the nuisance parameter, i.e.
1379 // Constraint (gobs, func( nuispar) ) and the condition is satisfied for
1380 // gobs = func( nuispar) where nunispar is at the MLE value
1381
1382
1383 if (model.GetGlobalObservables() && model.GetGlobalObservables()->getSize() > 0) {
1384
1385 if (verbose>1) {
1386 std::cout << "Generating Asimov data for global observables " << std::endl;
1387 }
1388
1389 RooArgSet gobs(*model.GetGlobalObservables());
1390
1391 // snapshot data global observables
1392 RooArgSet snapGlobalObsData;
1393 SetAllConstant(gobs, true);
1394 gobs.snapshot(snapGlobalObsData);
1395
1396
1397 RooArgSet nuis;
1398 if (model.GetNuisanceParameters()) nuis.add(*model.GetNuisanceParameters());
1399 if (nuis.getSize() == 0) {
1400 oocoutW((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData: model does not have nuisance parameters but has global observables"
1401 << " set global observables to model values " << endl;
1402 asimovGlobObs = gobs;
1403 return asimov;
1404 }
1405
1406 // part 1: create the nuisance pdf
1407 std::unique_ptr<RooAbsPdf> nuispdf(RooStats::MakeNuisancePdf(model,"TempNuisPdf") );
1408 if (nuispdf.get() == 0) {
1409 oocoutF((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData: model has nuisance parameters and global obs but no nuisance pdf "
1410 << std::endl;
1411 }
1412 // unfold the nuisance pdf if it is a prod pdf
1413 RooArgList pdfList;
1414 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(nuispdf.get());
1415 if (prod ) {
1416 pdfList.add(prod->pdfList());
1417 }
1418 else
1419 // nothing to unfold - just use the pdf
1420 pdfList.add(*nuispdf.get());
1421
1422 RooLinkedListIter iter(pdfList.iterator());
1423 for (RooAbsArg *a = (RooAbsArg *) iter.Next(); a != 0; a = (RooAbsArg *) iter.Next()) {
1424 RooAbsPdf *cterm = dynamic_cast<RooAbsPdf *>(a);
1425 assert(cterm && "AsimovUtils: a factor of the nuisance pdf is not a Pdf!");
1426 if (!cterm->dependsOn(nuis)) continue; // dummy constraints
1427 // skip also the case of uniform components
1428 if (typeid(*cterm) == typeid(RooUniform)) continue;
1429
1430 std::unique_ptr<RooArgSet> cpars(cterm->getParameters(&gobs));
1431 std::unique_ptr<RooArgSet> cgobs(cterm->getObservables(&gobs));
1432 if (cgobs->getSize() > 1) {
1433 oocoutE((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData: constraint term " << cterm->GetName()
1434 << " has multiple global observables -cannot generate - skip it" << std::endl;
1435 continue;
1436 }
1437 else if (cgobs->getSize() == 0) {
1438 oocoutW((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData: constraint term " << cterm->GetName()
1439 << " has no global observables - skip it" << std::endl;
1440 continue;
1441 }
1442 // the variable representing the global observable
1443 RooRealVar &rrv = dynamic_cast<RooRealVar &>(*cgobs->first());
1444
1445 // remove the constant parameters in cpars
1447 if (cpars->getSize() != 1) {
1448 oocoutE((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData:constraint term "
1449 << cterm->GetName() << " has multiple floating params - cannot generate - skip it " << std::endl;
1450 continue;
1451 }
1452
1453 bool foundServer = false;
1454 // note : this will work only for this type of constraints
1455 // expressed as RooPoisson, RooGaussian, RooLognormal, RooGamma
1456 TClass * cClass = cterm->IsA();
1457 if (verbose > 2) std::cout << "Constraint " << cterm->GetName() << " of type " << cClass->GetName() << std::endl;
1458 if ( cClass != RooGaussian::Class() && cClass != RooPoisson::Class() &&
1459 cClass != RooGamma::Class() && cClass != RooLognormal::Class() &&
1460 cClass != RooBifurGauss::Class() ) {
1461 TString className = (cClass) ? cClass->GetName() : "undefined";
1462 oocoutW((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData:constraint term "
1463 << cterm->GetName() << " of type " << className
1464 << " is a non-supported type - result might be not correct " << std::endl;
1465 }
1466
1467 // in case of a Poisson constraint make sure the rounding is not set
1468 if (cClass == RooPoisson::Class() ) {
1469 RooPoisson * pois = dynamic_cast<RooPoisson*>(cterm);
1470 assert(pois);
1471 pois->setNoRounding(true);
1472 }
1473
1474 // look at server of the constraint term and check if the global observable is part of the server
1475 RooAbsArg * arg = cterm->findServer(rrv);
1476 if (!arg) {
1477 // special case is for the Gamma where one might define the global observable n and you have a Gamma(b, n+1, ...._
1478 // 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
1479 // so in case of the Gamma ignore this test
1480 if ( cClass != RooGamma::Class() ) {
1481 oocoutE((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData:constraint term "
1482 << cterm->GetName() << " has no direct dependence on global observable- cannot generate it " << std::endl;
1483 continue;
1484 }
1485 }
1486
1487 // loop on the server of the constraint term
1488 // need to treat the Gamma as a special case
1489 // the mode of the Gamma is (k-1)*theta where theta is the inverse of the rate parameter.
1490 // 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
1491 RooAbsReal * thetaGamma = 0;
1492 if ( cClass == RooGamma::Class() ) {
1493 RooFIter itc(cterm->serverMIterator() );
1494 for (RooAbsArg *a2 = itc.next(); a2 != 0; a2 = itc.next()) {
1495 if (TString(a2->GetName()).Contains("theta") ) {
1496 thetaGamma = dynamic_cast<RooAbsReal*>(a2);
1497 break;
1498 }
1499 }
1500 if (thetaGamma == 0) {
1501 oocoutI((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData:constraint term "
1502 << cterm->GetName() << " is a Gamma distribution and no server named theta is found. Assume that the Gamma scale is 1 " << std::endl;
1503 }
1504 else {
1505 if (verbose>2)
1506 std::cout << "Gamma constraint has a scale " << thetaGamma->GetName() << " = " << thetaGamma->getVal() << std::endl;
1507 }
1508 }
1509 RooFIter iter2(cterm->serverMIterator() );
1510 for (RooAbsArg *a2 = iter2.next(); a2 != 0; a2 = iter2.next()) {
1511 RooAbsReal * rrv2 = dynamic_cast<RooAbsReal *>(a2);
1512 if (verbose > 2) std::cout << "Loop on constraint server term " << a2->GetName() << std::endl;
1513 if (rrv2 && rrv2->dependsOn(nuis) ) {
1514
1515
1516 // found server depending on nuisance
1517 if (foundServer) {
1518 oocoutE((TObject*)0,Generation) << "AsymptoticCalculator::MakeAsimovData:constraint term "
1519 << cterm->GetName() << " constraint term has more server depending on nuisance- cannot generate it " <<
1520 std::endl;
1521 foundServer = false;
1522 break;
1523 }
1524 if (thetaGamma && thetaGamma->getVal() > 0)
1525 rrv.setVal( rrv2->getVal() / thetaGamma->getVal() );
1526 else
1527 rrv.setVal( rrv2->getVal() );
1528 foundServer = true;
1529
1530 if (verbose>2)
1531 std::cout << "setting global observable "<< rrv.GetName() << " to value " << rrv.getVal()
1532 << " which comes from " << rrv2->GetName() << std::endl;
1533 }
1534 }
1535
1536 if (!foundServer) {
1537 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;
1538 std::cerr << "Parameters: " << std::endl;
1539 cpars->Print("V");
1540 std::cerr << "Observables: " << std::endl;
1541 cgobs->Print("V");
1542 }
1543// rrv.setVal(match->getVal());
1544 }
1545
1546 // make a snapshot of global observables
1547 // needed this ?? (LM)
1548
1549 asimovGlobObs.removeAll();
1550 SetAllConstant(gobs, true);
1551 gobs.snapshot(asimovGlobObs);
1552
1553 // revert global observables to the data value
1554 gobs = snapGlobalObsData;
1555
1556 if (verbose>0) {
1557 std::cout << "Generated Asimov data for global observables ";
1558 if (verbose == 1) gobs.Print();
1559 }
1560
1561 if (verbose > 1) {
1562 std::cout << "\nGlobal observables for data: " << std::endl;
1563 gobs.Print("V");
1564 std::cout << "\nGlobal observables for asimov: " << std::endl;
1565 asimovGlobObs.Print("V");
1566 }
1567
1568
1569 }
1570
1571 return asimov;
1572
1573}
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:90
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
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 (of arbitrary type) an...
Definition: RooAbsArg.h:73
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:276
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:730
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:302
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) 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:544
RooFIter serverMIterator() const
Definition: RooAbsArg.h:155
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:1909
Bool_t isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:342
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.
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
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:44
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
virtual Double_t sumEntries() const =0
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:175
virtual void add(const RooArgSet &row, Double_t weight=1, Double_t weightError=0)=0
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:918
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events from this p.d.f for use in extended likelihood calculations.
Definition: RooAbsPdf.cxx:3283
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:60
static void setHideOffset(Bool_t flag)
Definition: RooAbsReal.cxx:116
static Bool_t hideOffset()
Definition: RooAbsReal.cxx:117
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:126
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:23
virtual value_type getCurrentIndex() const override final
Return current index.
Definition: RooCategory.h:35
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:33
virtual Double_t sumEntries() const override
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
Definition: RooFitResult.h:98
Plain Gaussian p.d.f.
Definition: RooGaussian.h:25
A wrapper around TIterator derivatives.
TObject * Next() override
RooLinkedList is an collection class for internal use, storing a collection of RooAbsArg pointers in ...
Definition: RooLinkedList.h:35
virtual void Add(TObject *arg)
Definition: RooLinkedList.h:62
RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between th...
Definition: RooMinimizer.h:38
void setEvalErrorWall(Bool_t flag)
Definition: RooMinimizer.h:50
RooFitResult * save(const char *name=0, const char *title=0)
Save and return a RooFitResult snaphot 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)
Definition: RooPoisson.h:33
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
const RooArgList & pdfList() const
Definition: RooProdPdf.h:67
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name 'name' for this variable.
Definition: RooRealVar.cxx:418
Double_t getError() const
Definition: RooRealVar.h:62
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:261
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 * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
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:131
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition: TString.cxx:2311
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
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(const RooArgSet &globs)
RooCmdArg Hesse(Bool_t flag=kTRUE)
RooCmdArg CloneData(Bool_t flag)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg ConditionalObservables(const RooArgSet &set)
RooCmdArg Offset(Bool_t flag=kTRUE)
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:65
@ Minimization
Definition: RooGlobalFunc.h:67
@ Generation
Definition: RooGlobalFunc.h:67
@ InputArguments
Definition: RooGlobalFunc.h:68
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:882
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:891
Definition: first.py:1
auto * a
Definition: textangle.C:12