Logo ROOT   master
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 
14 Hypothesis Test Calculator based on the asymptotic formulae for the profile
15 likelihood ratio.
16 
17 It performs hypothesis tests using the asymptotic formula for the profile likelihood, and
18 uses the Asimov data set to compute expected significances or limits.
19 
20 See G. Cowan, K. Cranmer, E. Gross and O. Vitells: Asymptotic formulae for
21 likelihood- based tests of new physics. Eur. Phys. J., C71:1–19, 2011.
22 It provides methods to perform hypothesis tests using the likelihood function,
23 and computes the \f$p\f$-values for the null and the alternate hypothesis using the asymptotic
24 formulae for the profile likelihood ratio described in the given paper.
25 
26 The calculator provides methods to produce the Asimov dataset, *i.e.* a dataset
27 generated where the observed values are equal to the expected ones.
28 The Asimov data set is then used to compute the observed asymptotic \f$p\f$-value for
29 the alternate hypothesis and the asymptotic expected \f$p\f$-values.
30 
31 The asymptotic formulae are valid only for one POI (parameter of interest). So
32 the calculator works only for one-dimensional (one POI) models.
33 If more than one POI exists, only the first one is used.
34 
35 The 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 
42 #include "RooStats/ModelConfig.h"
44 #include "RooStats/RooStatsUtils.h"
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"
56 #include "Math/MinimizerOptions.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 
67 #include "Math/BrentRootFinder.h"
68 #include "Math/WrappedFunction.h"
69 
70 #include "TStopwatch.h"
71 
72 using namespace RooStats;
73 using namespace std;
74 
75 
77 
79 
80 ////////////////////////////////////////////////////////////////////////////////
81 /// set print level (static function)
82 ///
83 /// - 0 minimal,
84 /// - 1 normal,
85 /// - 2 debug
86 
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
245  fAsimovData = MakeAsimovData( *GetNullModel(), nominalParams, fAsimovGlobObs);
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 
291 Double_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 
735  if (fOneSided || fOneSidedDiscovery) {
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 
794 struct 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 
808 double 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 
857 void 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 
1054 RooAbsData * 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 
1216 RooAbsData * 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();
1257  int minimPrintLevel = ROOT::Math::MinimizerOptions::DefaultPrintLevel();
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 
1304  if (verbose < 2) RooMsgService::instance().setGlobalKillBelow(msglevel);
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 
1337 RooAbsData * 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 }
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:916
virtual Double_t sumEntries() const =0
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:1910
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...
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooCmdArg Offset(Bool_t flag=kTRUE)
Poisson pdf.
Definition: RooPoisson.h:19
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
RooCmdArg WeightVar(const char *name, Bool_t reinterpretAsWeight=kFALSE)
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
static double EvaluateNLL(RooAbsPdf &pdf, RooAbsData &data, const RooArgSet *condObs, const RooArgSet *globObs, const RooArgSet *poiSet=0)
const ModelConfig * GetAlternateModel(void) const
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:249
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
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
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
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).
void optimizeConst(Int_t flag)
If flag is true, perform constant term optimization on function being minimized.
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:731
RooCmdArg Hesse(Bool_t flag=kTRUE)
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:294
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 ...
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:891
RooProdPdf is an efficient implementation of a product of PDFs of the form .
Definition: RooProdPdf.h:31
#define oocoutF(o, a)
Definition: RooMsgService.h:49
#define oocoutI(o, a)
Definition: RooMsgService.h:45
HypoTestResult is a base class for results from hypothesis tests.
RooFit::MsgLevel globalKillBelow() const
virtual Double_t sumEntries() const override
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
Basic string class.
Definition: TString.h:131
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
virtual Bool_t setIndex(Int_t index, bool printError=true) override
Set value by specifying the index code of the desired state.
#define f(i)
Definition: RSha256.hxx:104
Bool_t IsNaN(Double_t x)
Definition: TMath.h:882
static RooMsgService & instance()
Return reference to singleton instance.
STL namespace.
AsymptoticCalculator(RooAbsData &data, const ModelConfig &altModel, const ModelConfig &nullModel, bool nominalAsimov=false)
constructor for asymptotic calculator from Data set and ModelConfig
RooCmdArg GlobalObservables(const RooArgSet &globs)
static void setHideOffset(Bool_t flag)
Definition: RooAbsReal.cxx:116
TRObject operator()(const T1 &t1) const
void setEps(Double_t eps)
Change MINUIT epsilon.
void setStrategy(Int_t strat)
Change MINUIT strategy to istrat.
RooAbsArg * findServer(const char *name) const
Return server of this with name name. Returns nullptr if not found.
Definition: RooAbsArg.h:203
Template class to wrap any C++ callable object which takes one argument i.e.
static void SetPrintLevel(int level)
set print level (static function)
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.
TObject * Next() override
Common base class for the Hypothesis Test Calculators.
double sqrt(double)
static RooAbsData * MakeAsimovData(RooAbsData &data, const ModelConfig &model, const RooArgSet &poiValues, RooArgSet &globObs, const RooArgSet *genPoiValues=0)
Make Asimov data.
#define oocoutP(o, a)
Definition: RooMsgService.h:46
RooFIter serverMIterator() const
Definition: RooAbsArg.h:155
Double_t x[n]
Definition: legend1.C:17
static void FillBins(const RooAbsPdf &pdf, const RooArgList &obs, RooAbsData &data, int &index, double &binVolume, int &ibin)
fill bins by looping recursively on observables
Int_t numTypes(const char *=0) const
Return number of types defined (in range named rangeName if rangeName!=0)
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
void Class()
Definition: Class.C:29
RooCmdArg PrintLevel(Int_t code)
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
#define oocoutE(o, a)
Definition: RooMsgService.h:48
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
Plain Gaussian p.d.f.
Definition: RooGaussian.h:25
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:320
virtual HypoTestResult * GetHypoTest() const
re-implement HypoTest computation using the asymptotic
virtual const char * getCurrentLabel() const
Return label string of current state.
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name &#39;name&#39; for this variable.
Definition: RooRealVar.cxx:382
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:3288
virtual void add(const RooArgSet &row, Double_t weight=1, Double_t weightError=0)=0
RooCmdArg Constrain(const RooArgSet &params)
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:225
virtual void Add(TObject *arg)
Definition: RooLinkedList.h:62
double Root() const
Returns root value.
const RooArgSet * GetConditionalObservables() const
get RooArgSet for conditional observables (return NULL if not existing)
Definition: ModelConfig.h:252
bool Initialize() const
initialize the calculator by performing a global fit and make the Asimov data set ...
Int_t getSize() const
ROOT::R::TRInterface & r
Definition: Object.C:4
RooCmdArg EvalErrorWall(Bool_t flag)
static const std::string & DefaultMinimizerType()
const RooAbsCategoryLValue & indexCat() const
RooAbsArg * first() const
void setConstant(Bool_t value=kTRUE)
auto * a
Definition: textangle.C:12
RooFitResult * save(const char *name=0, const char *title=0)
Save and return a RooFitResult snaphot of current minimizer status.
const RooArgList & pdfList() const
Definition: RooProdPdf.h:67
RooCmdArg ConditionalObservables(const RooArgSet &set)
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
Double_t minNll() const
Definition: RooFitResult.h:98
void setGlobalKillBelow(RooFit::MsgLevel level)
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition: TClass.h:80
static bool SetObsToExpected(RooAbsPdf &pdf, const RooArgSet &obs)
set observed value to the expected one works for Gaussian, Poisson or LogNormal assumes mean paramete...
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
void setNoRounding(bool flag=kTRUE)
Definition: RooPoisson.h:33
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:23
Flat p.d.f.
Definition: RooUniform.h:24
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)...
static Bool_t hideOffset()
Definition: RooAbsReal.cxx:117
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
virtual Double_t CLs() const
is simply (not a method, but a quantity)
const Bool_t kFALSE
Definition: RtypesCore.h:90
Class for finding the root of a one dimensional function using the Brent algorithm.
RooCmdArg Strategy(Int_t code)
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
RooCmdArg CloneData(Bool_t flag)
const RooAbsData * GetData(void) const
static const std::string & DefaultMinimizerAlgo()
RooLinkedList is an collection class for internal use, storing a collection of RooAbsArg pointers in ...
Definition: RooLinkedList.h:35
Namespace for the RooStats classes.
Definition: Asimov.h:19
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:175
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
#define ClassImp(name)
Definition: Rtypes.h:361
void Print(std::ostream &os, const OptionType &opt)
void setEvalErrorWall(Bool_t flag)
Definition: RooMinimizer.h:50
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
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&#39;t match any of...
Definition: RooAbsArg.cxx:545
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
RooCmdArg Import(const char *state, TH1 &histo)
int fUseQTilde
flag to check if calculator is initialized
bool SetAllConstant(const RooAbsCollection &coll, bool constant=true)
Definition: RooStatsUtils.h:89
Int_t minimize(const char *type, const char *alg=0)
Minimise the function passed in the constructor.
const ModelConfig * GetNullModel(void) const
bool IsNLLOffset()
Test of RooStats should by default offset NLL calculations.
#define oocoutW(o, a)
Definition: RooMsgService.h:47
Mother of all ROOT objects.
Definition: TObject.h:37
Int_t setPrintLevel(Int_t newLevel)
Change the MINUIT internal printing level.
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:69
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return NULL if not existing)
Definition: ModelConfig.h:255
RooCmdArg Index(RooCategory &icat)
RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between th...
Definition: RooMinimizer.h:38
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
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:1259
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
virtual Int_t getBins(const char *name=0) const
Get number of bins of currently defined range.
bool SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup)
Sets the function for the rest of the algorithms.
const RooArgSet * GetSnapshot() const
get RooArgSet for parameters for a particular hypothesis (return NULL if not existing) ...
Bool_t contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
Definition: first.py:1
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
RooCmdArg Minimizer(const char *type, const char *alg=0)
Double_t getError() const
Definition: RooRealVar.h:62
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
const Bool_t kTRUE
Definition: RtypesCore.h:89
Hypothesis Test Calculator based on the asymptotic formulae for the profile likelihood ratio...
virtual value_type getCurrentIndex() const override final
Return current index.
Definition: RooCategory.h:35
A wrapper around TIterator derivatives.
Bool_t isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:360
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
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 &#39;data&#39; argset, to the data set...
Stopwatch class.
Definition: TStopwatch.h:28
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.