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