Logo ROOT   6.10/09
Reference Guide
HypoTestInverter.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
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::HypoTestInverter
12  \ingroup Roostats
13 
14 HypoTestInverter class for performing an hypothesis test inversion by scanning
15 the hypothesis test results of an HypoTestCalculator for various values of the
16 parameter of interest. By looking at the confidence level curve of the result an
17 upper limit, where it intersects the desired confidence level, can be derived.
18 The class implements the RooStats::IntervalCalculator interface and returns an
19 RooStats::HypoTestInverterResult class. The result is a SimpleInterval, which
20 via the method UpperLimit returns to the user the upper limit value.
21 
22 The HypoTestInverter implements various option for performing the scan.
23 HypoTestInverter::RunFixedScan will scan using a fixed grid the parameter of
24 interest. HypoTestInverter::RunAutoScan will perform an automatic scan to find
25 optimally the curve and it will stop until the desired precision is obtained.
26 The confidence level value at a given point can be done via HypoTestInverter::RunOnePoint.
27 The class can scan the CLs+b values or alternatively CLs
28 (if the method HypoTestInverter::UseCLs has been called).
29 
30 Contributions to this class have been written by Giovanni Petrucciani and Annapaola Decosa
31 */
32 
33 // include other header files
34 
35 #include "RooAbsData.h"
36 #
37 #include "TMath.h"
38 
39 #include "RooStats/HybridResult.h"
40 
42 #include <cassert>
43 #include <cmath>
44 
45 #include "TF1.h"
46 #include "TFile.h"
47 #include "TH1.h"
48 #include "TLine.h"
49 #include "TCanvas.h"
50 #include "TGraphErrors.h"
51 #include "RooRealVar.h"
52 #include "RooArgSet.h"
53 #include "RooAbsPdf.h"
54 #include "RooRandom.h"
55 #include "RooConstVar.h"
56 #include "RooMsgService.h"
57 #include "RooStats/ModelConfig.h"
64 #include "RooStats/ToyMCSampler.h"
65 #include "RooStats/HypoTestPlot.h"
67 
68 #include "RooStats/ProofConfig.h"
69 
71 
72 using namespace RooStats;
73 using namespace std;
74 
75 // static variable definitions
76 double HypoTestInverter::fgCLAccuracy = 0.005;
77 unsigned int HypoTestInverter::fgNToys = 500;
78 
81 std::string HypoTestInverter::fgAlgo = "logSecant";
82 
84 
85 // helper class to wrap the functionality of the various HypoTestCalculators
86 
87 template<class HypoTestType>
88 struct HypoTestWrapper {
89 
90  static void SetToys(HypoTestType * h, int toyNull, int toyAlt) { h->SetToys(toyNull,toyAlt); }
91 
92 };
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 /// set flag to close proof for every new run
96 
98  fgCloseProof = flag;
99 }
100 
101 ////////////////////////////////////////////////////////////////////////////////
102 /// get the variable to scan
103 /// try first with null model if not go to alternate model
104 
106 
107  RooRealVar * varToScan = 0;
108  const ModelConfig * mc = hc.GetNullModel();
109  if (mc) {
110  const RooArgSet * poi = mc->GetParametersOfInterest();
111  if (poi) varToScan = dynamic_cast<RooRealVar*> (poi->first() );
112  }
113  if (!varToScan) {
114  mc = hc.GetAlternateModel();
115  if (mc) {
116  const RooArgSet * poi = mc->GetParametersOfInterest();
117  if (poi) varToScan = dynamic_cast<RooRealVar*> (poi->first() );
118  }
119  }
120  return varToScan;
121 }
122 
123 ////////////////////////////////////////////////////////////////////////////////
124 /// check the model given the given hypotestcalculator
125 
127  const ModelConfig * modelSB = hc.GetNullModel();
128  const ModelConfig * modelB = hc.GetAlternateModel();
129  if (!modelSB || ! modelB)
130  oocoutF((TObject*)0,InputArguments) << "HypoTestInverter - model are not existing" << std::endl;
131  assert(modelSB && modelB);
132 
133  oocoutI((TObject*)0,InputArguments) << "HypoTestInverter ---- Input models: \n"
134  << "\t\t using as S+B (null) model : "
135  << modelSB->GetName() << "\n"
136  << "\t\t using as B (alternate) model : "
137  << modelB->GetName() << "\n" << std::endl;
138 
139  // check if scanVariable is included in B model pdf
140  RooAbsPdf * bPdf = modelB->GetPdf();
141  const RooArgSet * bObs = modelB->GetObservables();
142  if (!bPdf || !bObs) {
143  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - B model has no pdf or observables defined" << std::endl;
144  return;
145  }
146  RooArgSet * bParams = bPdf->getParameters(*bObs);
147  if (!bParams) {
148  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - pdf of B model has no parameters" << std::endl;
149  return;
150  }
151  if (bParams->find(scanVariable.GetName() ) ) {
152  const RooArgSet * poiB = modelB->GetSnapshot();
153  if (!poiB || !poiB->find(scanVariable.GetName()) ||
154  ( (RooRealVar*) poiB->find(scanVariable.GetName()) )->getVal() != 0 )
155  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter - using a B model with POI "
156  << scanVariable.GetName() << " not equal to zero "
157  << " user must check input model configurations " << endl;
158  if (poiB) delete poiB;
159  }
160  delete bParams;
161 }
162 
163 ////////////////////////////////////////////////////////////////////////////////
164 /// default constructor (doesn't do anything)
165 
167  fTotalToysRun(0),
168  fMaxToys(0),
169  fCalculator0(0),
170  fScannedVariable(0),
171  fResults(0),
172  fUseCLs(false),
173  fScanLog(false),
174  fSize(0),
175  fVerbose(0),
176  fCalcType(kUndefined),
177  fNBins(0), fXmin(1), fXmax(1),
178  fNumErr(0)
179 {
180 }
181 
182 ////////////////////////////////////////////////////////////////////////////////
183 /// Constructor from a HypoTestCalculatorGeneric
184 /// The HypoTest calculator must be a FrequentistCalculator or HybridCalculator type
185 /// Other type of calculators are not supported.
186 /// The calculator must be created before by using the S+B model for the null and
187 /// the B model for the alt
188 /// If no variable to scan are given they are assumed to be the first variable
189 /// from the parameter of interests of the null model
190 
192  RooRealVar* scannedVariable, double size ) :
193  fTotalToysRun(0),
194  fMaxToys(0),
195  fCalculator0(0),
196  fScannedVariable(scannedVariable),
197  fResults(0),
198  fUseCLs(false),
199  fScanLog(false),
200  fSize(size),
201  fVerbose(0),
203  fNBins(0), fXmin(1), fXmax(1),
204  fNumErr(0)
205 {
206 
207  if (!fScannedVariable) {
209  }
210  if (!fScannedVariable)
211  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
212  else
214 
215  HybridCalculator * hybCalc = dynamic_cast<HybridCalculator*>(&hc);
216  if (hybCalc) {
217  fCalcType = kHybrid;
218  fCalculator0 = hybCalc;
219  return;
220  }
221  FrequentistCalculator * freqCalc = dynamic_cast<FrequentistCalculator*>(&hc);
222  if (freqCalc) {
224  fCalculator0 = freqCalc;
225  return;
226  }
227  AsymptoticCalculator * asymCalc = dynamic_cast<AsymptoticCalculator*>(&hc);
228  if (asymCalc) {
230  fCalculator0 = asymCalc;
231  return;
232  }
233  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Type of hypotest calculator is not supported " <<std::endl;
234  fCalculator0 = &hc;
235 }
236 
237 ////////////////////////////////////////////////////////////////////////////////
238 /// Constructor from a reference to a HybridCalculator
239 /// The calculator must be created before by using the S+B model for the null and
240 /// the B model for the alt
241 /// If no variable to scan are given they are assumed to be the first variable
242 /// from the parameter of interests of the null model
243 
245  RooRealVar* scannedVariable, double size ) :
246  fTotalToysRun(0),
247  fMaxToys(0),
248  fCalculator0(&hc),
249  fScannedVariable(scannedVariable),
250  fResults(0),
251  fUseCLs(false),
252  fScanLog(false),
253  fSize(size),
254  fVerbose(0),
256  fNBins(0), fXmin(1), fXmax(1),
257  fNumErr(0)
258 {
259 
260  if (!fScannedVariable) {
262  }
263  if (!fScannedVariable)
264  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
265  else
267 
268 }
269 
270 ////////////////////////////////////////////////////////////////////////////////
271 /// Constructor from a reference to a FrequentistCalculator
272 /// The calculator must be created before by using the S+B model for the null and
273 /// the B model for the alt
274 /// If no variable to scan are given they are assumed to be the first variable
275 /// from the parameter of interests of the null model
276 
278  RooRealVar* scannedVariable, double size ) :
279  fTotalToysRun(0),
280  fMaxToys(0),
281  fCalculator0(&hc),
282  fScannedVariable(scannedVariable),
283  fResults(0),
284  fUseCLs(false),
285  fScanLog(false),
286  fSize(size),
287  fVerbose(0),
289  fNBins(0), fXmin(1), fXmax(1),
290  fNumErr(0)
291 {
292 
293  if (!fScannedVariable) {
295  }
296  if (!fScannedVariable)
297  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
298  else
300 }
301 
302 ////////////////////////////////////////////////////////////////////////////////
303 /// Constructor from a reference to a AsymptoticCalculator
304 /// The calculator must be created before by using the S+B model for the null and
305 /// the B model for the alt
306 /// If no variable to scan are given they are assumed to be the first variable
307 /// from the parameter of interests of the null model
308 
310  RooRealVar* scannedVariable, double size ) :
311  fTotalToysRun(0),
312  fMaxToys(0),
313  fCalculator0(&hc),
314  fScannedVariable(scannedVariable),
315  fResults(0),
316  fUseCLs(false),
317  fScanLog(false),
318  fSize(size),
319  fVerbose(0),
321  fNBins(0), fXmin(1), fXmax(1),
322  fNumErr(0)
323 {
324 
325  if (!fScannedVariable) {
327  }
328  if (!fScannedVariable)
329  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
330  else
332 
333 }
334 
335 ////////////////////////////////////////////////////////////////////////////////
336 /// Constructor from a model for B model and a model for S+B.
337 /// An HypoTestCalculator (Hybrid of Frequentis) will be created using the
338 /// S+B model as the null and the B model as the alternate
339 /// If no variable to scan are given they are assumed to be the first variable
340 /// from the parameter of interests of the null model
341 
343  RooRealVar * scannedVariable, ECalculatorType type, double size) :
344  fTotalToysRun(0),
345  fMaxToys(0),
346  fCalculator0(0),
347  fScannedVariable(scannedVariable),
348  fResults(0),
349  fUseCLs(false),
350  fScanLog(false),
351  fSize(size),
352  fVerbose(0),
353  fCalcType(type),
354  fNBins(0), fXmin(1), fXmax(1),
355  fNumErr(0)
356 {
357  if(fCalcType==kFrequentist) fHC.reset(new FrequentistCalculator(data, bModel, sbModel));
358  if(fCalcType==kHybrid) fHC.reset( new HybridCalculator(data, bModel, sbModel)) ;
359  if(fCalcType==kAsymptotic) fHC.reset( new AsymptoticCalculator(data, bModel, sbModel));
360  fCalculator0 = fHC.get();
361  // get scanned variable
362  if (!fScannedVariable) {
364  }
365  if (!fScannedVariable)
366  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
367  else
369 
370 }
371 
372 ////////////////////////////////////////////////////////////////////////////////
373 /// copy-constructor
374 /// NOTE: this class does not copy the contained result and
375 /// the HypoTestCalculator, but only the pointers
376 /// It requires the original HTI to be alive
377 
380  fTotalToysRun(0),
381  fCalculator0(0), fScannedVariable(0), // add these for Coverity
382  fResults(0)
383 {
384  (*this) = rhs;
385 }
386 
387 ////////////////////////////////////////////////////////////////////////////////
388 /// assignment operator
389 /// NOTE: this class does not copy the contained result and
390 /// the HypoTestCalculator, but only the pointers
391 /// It requires the original HTI to be alive
392 
394  if (this == &rhs) return *this;
395  fTotalToysRun = 0;
396  fMaxToys = rhs.fMaxToys;
399  fUseCLs = rhs.fUseCLs;
400  fScanLog = rhs.fScanLog;
401  fSize = rhs.fSize;
402  fVerbose = rhs.fVerbose;
403  fCalcType = rhs.fCalcType;
404  fNBins = rhs.fNBins;
405  fXmin = rhs.fXmin;
406  fXmax = rhs.fXmax;
407  fNumErr = rhs.fNumErr;
408 
409  return *this;
410 }
411 
412 ////////////////////////////////////////////////////////////////////////////////
413 /// destructor (delete the HypoTestInverterResult)
414 
416 {
417  if (fResults) delete fResults;
418  fCalculator0 = 0;
419 }
420 
421 ////////////////////////////////////////////////////////////////////////////////
422 /// return the test statistic which is or will be used by the class
423 
425 {
428  }
429  else
430  return 0;
431 }
432 
433 ////////////////////////////////////////////////////////////////////////////////
434 /// set the test statistic to use
435 
437 {
440  return true;
441  }
442  else return false;
443 }
444 
445 ////////////////////////////////////////////////////////////////////////////////
446 /// delete contained result and graph
447 
449  if (fResults) delete fResults;
450  fResults = 0;
451  fLimitPlot.reset(nullptr);
452 }
453 
454 ////////////////////////////////////////////////////////////////////////////////
455 /// create a new HypoTestInverterResult to hold all computed results
456 
458  if (fResults == 0) {
459  TString results_name = "result_";
460  results_name += fScannedVariable->GetName();
462  TString title = "HypoTestInverter Result For ";
463  title += fScannedVariable->GetName();
464  fResults->SetTitle(title);
465  }
468  // check if one or two sided scan
469  if (fCalculator0) {
470  // if asymptotic calculator
472  if (ac)
474  else {
475  // in case of the other calculators
477  if (sampler) {
478  ProfileLikelihoodTestStat * pl = dynamic_cast<ProfileLikelihoodTestStat*>(sampler->GetTestStatistic());
479  if (pl && pl->IsTwoSided() ) fResults->fIsTwoSided = true;
480  }
481  }
482  }
483 }
484 
485 ////////////////////////////////////////////////////////////////////////////////
486 /// Run a fixed scan or the automatic scan depending on the configuration
487 /// Return if needed a copy of the result object which will be managed by the user
488 
490 
491  // if having a result with at least one point return it
492  if (fResults && fResults->ArraySize() >= 1) {
493  oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - return an already existing interval " << std::endl;
494  return (HypoTestInverterResult*)(fResults->Clone());
495  }
496 
497  if (fNBins > 0) {
498  oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - run a fixed scan" << std::endl;
499  bool ret = RunFixedScan(fNBins, fXmin, fXmax, fScanLog);
500  if (!ret)
501  oocoutE((TObject*)0,Eval) << "HypoTestInverter::GetInterval - error running a fixed scan " << std::endl;
502  }
503  else {
504  oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - run an automatic scan" << std::endl;
505  double limit(0),err(0);
506  bool ret = RunLimit(limit,err);
507  if (!ret)
508  oocoutE((TObject*)0,Eval) << "HypoTestInverter::GetInterval - error running an auto scan " << std::endl;
509  }
510 
512 
513  return (HypoTestInverterResult*) (fResults->Clone());
514 }
515 
516 ////////////////////////////////////////////////////////////////////////////////
517 /// Run the Hypothesis test at a previous configured point
518 /// (internal function called by RunOnePoint)
519 
520 HypoTestResult * HypoTestInverter::Eval(HypoTestCalculatorGeneric &hc, bool adaptive, double clsTarget) const {
521  //for debug
522  // std::cout << ">>>>>>>>>>> " << std::endl;
523  // std::cout << "alternate model " << std::endl;
524  // hc.GetAlternateModel()->GetNuisanceParameters()->Print("V");
525  // hc.GetAlternateModel()->GetParametersOfInterest()->Print("V");
526  // std::cout << "Null model " << std::endl;
527  // hc.GetNullModel()->GetNuisanceParameters()->Print("V");
528  // hc.GetNullModel()->GetParametersOfInterest()->Print("V");
529  // std::cout << "<<<<<<<<<<<<<<< " << std::endl;
530 
531  // run the hypothesis test
532  HypoTestResult * hcResult = hc.GetHypoTest();
533  if (hcResult == 0) {
534  oocoutE((TObject*)0,Eval) << "HypoTestInverter::Eval - HypoTest failed" << std::endl;
535  return hcResult;
536  }
537 
538  // since the b model is the alt need to set the flag
539  hcResult->SetBackgroundAsAlt(true);
540 
541 
542  // bool flipPvalue = false;
543  // if (flipPValues)
544  // hcResult->SetPValueIsRightTail(!hcResult->GetPValueIsRightTail());
545 
546  // adjust for some numerical error in discrete models and == is not anymore
547  if (hcResult->GetPValueIsRightTail() )
548  hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()-fNumErr); // issue with < vs <= in discrete models
549  else
550  hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()+fNumErr); // issue with < vs <= in discrete models
551 
552  double clsMid = (fUseCLs ? hcResult->CLs() : hcResult->CLsplusb());
553  double clsMidErr = (fUseCLs ? hcResult->CLsError() : hcResult->CLsplusbError());
554 
555  //if (fVerbose) std::cout << (fUseCLs ? "\tCLs = " : "\tCLsplusb = ") << clsMid << " +/- " << clsMidErr << std::endl;
556 
557  if (adaptive) {
558 
559  if (fCalcType == kHybrid) HypoTestWrapper<HybridCalculator>::SetToys((HybridCalculator*)&hc, fUseCLs ? fgNToys : 1, 4*fgNToys);
560  if (fCalcType == kFrequentist) HypoTestWrapper<FrequentistCalculator>::SetToys((FrequentistCalculator*)&hc, fUseCLs ? fgNToys : 1, 4*fgNToys);
561 
562  while (clsMidErr >= fgCLAccuracy && (clsTarget == -1 || fabs(clsMid-clsTarget) < 3*clsMidErr) ) {
563  std::unique_ptr<HypoTestResult> more(hc.GetHypoTest());
564 
565  // if (flipPValues)
566  // more->SetPValueIsRightTail(!more->GetPValueIsRightTail());
567 
568  hcResult->Append(more.get());
569  clsMid = (fUseCLs ? hcResult->CLs() : hcResult->CLsplusb());
570  clsMidErr = (fUseCLs ? hcResult->CLsError() : hcResult->CLsplusbError());
571  if (fVerbose) std::cout << (fUseCLs ? "\tCLs = " : "\tCLsplusb = ") << clsMid << " +/- " << clsMidErr << std::endl;
572  }
573 
574  }
575  if (fVerbose ) {
576  oocoutP((TObject*)0,Eval) << "P values for " << fScannedVariable->GetName() << " = " <<
577  fScannedVariable->getVal() << "\n" <<
578  "\tCLs = " << hcResult->CLs() << " +/- " << hcResult->CLsError() << "\n" <<
579  "\tCLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << "\n" <<
580  "\tCLsplusb = " << hcResult->CLsplusb() << " +/- " << hcResult->CLsplusbError() << "\n" <<
581  std::endl;
582  }
583 
584  if (fCalcType == kFrequentist || fCalcType == kHybrid) {
585  fTotalToysRun += (hcResult->GetAltDistribution()->GetSize() + hcResult->GetNullDistribution()->GetSize());
586 
587  // set sampling distribution name
588  TString nullDistName = TString::Format("%s_%s_%4.2f",hcResult->GetNullDistribution()->GetName(),
590  TString altDistName = TString::Format("%s_%s_%4.2f",hcResult->GetAltDistribution()->GetName(),
592 
593  hcResult->GetNullDistribution()->SetName(nullDistName);
594  hcResult->GetAltDistribution()->SetName(altDistName);
595  }
596 
597  return hcResult;
598 }
599 
600 ////////////////////////////////////////////////////////////////////////////////
601 /// Run a Fixed scan in npoints between min and max
602 
603 bool HypoTestInverter::RunFixedScan( int nBins, double xMin, double xMax, bool scanLog ) const
604 {
605 
606  CreateResults();
607  // interpolate the limits
608  fResults->fFittedLowerLimit = false;
609  fResults->fFittedUpperLimit = false;
610 
611  // safety checks
612  if ( nBins<=0 ) {
613  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - Please provide nBins>0\n";
614  return false;
615  }
616  if ( nBins==1 && xMin!=xMax ) {
617  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - nBins==1 -> I will run for xMin (" << xMin << ")\n";
618  }
619  if ( xMin==xMax && nBins>1 ) {
620  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMin==xMax -> I will enforce nBins==1\n";
621  nBins = 1;
622  }
623  if ( xMin>xMax ) {
624  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - Please provide xMin ("
625  << xMin << ") smaller that xMax (" << xMax << ")\n";
626  return false;
627  }
628 
629  if (xMin < fScannedVariable->getMin()) {
630  xMin = fScannedVariable->getMin();
631  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMin < lower bound, use xmin = "
632  << xMin << std::endl;
633  }
634  if (xMax > fScannedVariable->getMax()) {
635  xMax = fScannedVariable->getMax();
636  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMax > upper bound, use xmax = "
637  << xMax << std::endl;
638  }
639 
640  double thisX = xMin;
641  for (int i=0; i<nBins; i++) {
642 
643  if (i > 0) { // avoids case of nBins = 1
644  if (scanLog)
645  thisX = exp( log(xMin) + i*(log(xMax)-log(xMin))/(nBins-1) ); // scan in log x
646  else
647  thisX = xMin + i*(xMax-xMin)/(nBins-1); // linear scan in x
648  }
649 
650  bool status = RunOnePoint(thisX);
651 
652  // check if failed status
653  if ( status==false ) {
654  std::cout << "\t\tLoop interrupted because of failed status\n";
655  return false;
656  }
657  }
658 
659  return true;
660 }
661 
662 ////////////////////////////////////////////////////////////////////////////////
663 /// run only one point at the given POI value
664 
665 bool HypoTestInverter::RunOnePoint( double rVal, bool adaptive, double clTarget) const
666 {
667 
668  CreateResults();
669 
670  // check if rVal is in the range specified for fScannedVariable
671  if ( rVal < fScannedVariable->getMin() ) {
672  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunOnePoint - Out of range: using the lower bound "
674  << " on the scanned variable rather than " << rVal<< "\n";
675  rVal = fScannedVariable->getMin();
676  }
677  if ( rVal > fScannedVariable->getMax() ) {
678  // print a message when you have a significative difference since rval is computed
679  if ( rVal > fScannedVariable->getMax()*(1.+1.E-12) )
680  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunOnePoint - Out of range: using the upper bound "
682  << " on the scanned variable rather than " << rVal<< "\n";
683  rVal = fScannedVariable->getMax();
684  }
685 
686  // save old value
687  double oldValue = fScannedVariable->getVal();
688 
689  // evaluate hybrid calculator at a single point
690  fScannedVariable->setVal(rVal);
691  // need to set value of rval in hybridcalculator
692  // assume null model is S+B and alternate is B only
693  const ModelConfig * sbModel = fCalculator0->GetNullModel();
694  RooArgSet poi; poi.add(*sbModel->GetParametersOfInterest());
695  // set poi to right values
696  poi = RooArgSet(*fScannedVariable);
697  const_cast<ModelConfig*>(sbModel)->SetSnapshot(poi);
698 
699  if (fVerbose > 0)
700  oocoutP((TObject*)0,Eval) << "Running for " << fScannedVariable->GetName() << " = " << fScannedVariable->getVal() << endl;
701 
702  // compute the results
703  HypoTestResult* result = Eval(*fCalculator0,adaptive,clTarget);
704  if (!result) {
705  oocoutE((TObject*)0,Eval) << "HypoTestInverter - Error running point " << fScannedVariable->GetName() << " = " <<
706  fScannedVariable->getVal() << endl;
707  return false;
708  }
709  // in case of a dummy result
710  if (TMath::IsNaN(result->NullPValue() ) && TMath::IsNaN(result->AlternatePValue() ) ) {
711  oocoutW((TObject*)0,Eval) << "HypoTestInverter - Skip invalid result for point " << fScannedVariable->GetName() << " = " <<
712  fScannedVariable->getVal() << endl;
713  return true; // need to return true to avoid breaking the scan loop
714  }
715 
716  double lastXtested;
717  if ( fResults->ArraySize()!=0 ) lastXtested = fResults->GetXValue(fResults->ArraySize()-1);
718  else lastXtested = -999;
719 
720  if ( (std::abs(rVal) < 1 && TMath::AreEqualAbs(rVal, lastXtested,1.E-12) ) ||
721  (std::abs(rVal) >= 1 && TMath::AreEqualRel(rVal, lastXtested,1.E-12) ) ) {
722 
723  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunOnePoint - Merge with previous result for "
724  << fScannedVariable->GetName() << " = " << rVal << std::endl;
725  HypoTestResult* prevResult = fResults->GetResult(fResults->ArraySize()-1);
726  if (prevResult && prevResult->GetNullDistribution() && prevResult->GetAltDistribution()) {
727  prevResult->Append(result);
728  delete result; // we can delete the result
729  }
730  else {
731  // if it was empty we re-use it
732  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunOnePoint - replace previous empty result\n";
733  fResults->fYObjects.Remove( prevResult);
734  fResults->fYObjects.Add(result);
735  }
736 
737  } else {
738 
739  // fill the results in the HypoTestInverterResult array
740  fResults->fXValues.push_back(rVal);
741  fResults->fYObjects.Add(result);
742 
743  }
744 
745  // std::cout << "computed value for poi " << rVal << " : " << fResults->GetYValue(fResults->ArraySize()-1)
746  // << " +/- " << fResults->GetYError(fResults->ArraySize()-1) << endl;
747 
748  fScannedVariable->setVal(oldValue);
749 
750  return true;
751 }
752 
753 ////////////////////////////////////////////////////////////////////////////////
754 /// run an automatic scan until the desired accuracy is reached
755 /// Start by default from the full interval (min,max) of the POI and then via bisection find the line crossing
756 /// the target line
757 /// Optionally an hint can be provided and the scan will be done closer to that value
758 /// If by bisection the desired accuracy will not be reached a fit to the points is performed
759 
760 bool HypoTestInverter::RunLimit(double &limit, double &limitErr, double absAccuracy, double relAccuracy, const double*hint) const {
761 
762 
763  // routine from G. Petrucciani (from HiggsCombination CMS package)
764 // bool HybridNew::runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) {
765 
767  //w->loadSnapshot("clean");
768 
769  if ((hint != 0) && (*hint > r->getMin())) {
770  r->setMax(std::min<double>(3.0 * (*hint), r->getMax()));
771  r->setMin(std::max<double>(0.3 * (*hint), r->getMin()));
772  oocoutI((TObject*)0,InputArguments) << "HypoTestInverter::RunLimit - Use hint value " << *hint
773  << " search in interval " << r->getMin() << " , " << r->getMax() << std::endl;
774  }
775 
776  // if not specified use the default values for rel and absolute accuracy
777  if (absAccuracy <= 0) absAccuracy = fgAbsAccuracy;
778  if (relAccuracy <= 0) relAccuracy = fgRelAccuracy;
779 
780  typedef std::pair<double,double> CLs_t;
781  double clsTarget = fSize;
782  CLs_t clsMin(1,0), clsMax(0,0), clsMid(0,0);
783  double rMin = r->getMin(), rMax = r->getMax();
784  limit = 0.5*(rMax + rMin);
785  limitErr = 0.5*(rMax - rMin);
786  bool done = false;
787 
788  TF1 expoFit("expoFit","[0]*exp([1]*(x-[2]))", rMin, rMax);
789 
790  // if (readHybridResults_) {
791  // if (verbose > 0) std::cout << "Search for upper limit using pre-computed grid of p-values" << std::endl;
792 
793  // readAllToysFromFile();
794  // double minDist=1e3;
795  // for (int i = 0, n = limitPlot_->GetN(); i < n; ++i) {
796  // double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
797  // if (verbose > 0) std::cout << " r " << x << (CLs_ ? ", CLs = " : ", CLsplusb = ") << y << " +/- " << ey << std::endl;
798  // if (y-3*ey >= clsTarget) { rMin = x; clsMin = CLs_t(y,ey); }
799  // if (y+3*ey <= clsTarget) { rMax = x; clsMax = CLs_t(y,ey); }
800  // if (fabs(y-clsTarget) < minDist) { limit = x; minDist = fabs(y-clsTarget); }
801  // }
802  // if (verbose > 0) std::cout << " after scan x ~ " << limit << ", bounds [ " << rMin << ", " << rMax << "]" << std::endl;
803  // limitErr = std::max(limit-rMin, rMax-limit);
804  // expoFit.SetRange(rMin,rMax);
805 
806  // if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
807  // if (verbose > 1) std::cout << " reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
808  // done = true;
809  // }
810  // } else {
811 
812  fLimitPlot.reset(new TGraphErrors());
813 
814  if (fVerbose > 0) std::cout << "Search for upper limit to the limit" << std::endl;
815  for (int tries = 0; tries < 6; ++tries) {
816  //clsMax = eval(w, mc_s, mc_b, data, rMax);
817  if (! RunOnePoint(rMax) ) {
818  oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Hypotest failed" << std::endl;
819  return false;
820  }
821  clsMax = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
822  if (clsMax.first == 0 || clsMax.first + 3 * fabs(clsMax.second) < clsTarget ) break;
823  rMax += rMax;
824  if (tries == 5) {
825  oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Cannot set higher limit: at " << r->GetName()
826  << " = " << rMax << " still get "
827  << (fUseCLs ? "CLs" : "CLsplusb") << " = " << clsMax.first << std::endl;
828  return false;
829  }
830  }
831  if (fVerbose > 0) {
832  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Search for lower limit to the limit" << std::endl;
833  }
834  //clsMin = (fUseCLs && rMin == 0 ? CLs_t(1,0) : eval(w, mc_s, mc_b, data, rMin));
835  if ( fUseCLs && rMin == 0 ) {
836  clsMin = CLs_t(1,0);
837  }
838  else {
839  if (! RunOnePoint(rMin) ) return false;
840  clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
841  }
842  if (clsMin.first != 1 && clsMin.first - 3 * fabs(clsMin.second) < clsTarget) {
843  if (fUseCLs) {
844  rMin = 0;
845  clsMin = CLs_t(1,0); // this is always true for CLs
846  } else {
847  rMin = -rMax / 4;
848  for (int tries = 0; tries < 6; ++tries) {
849  if (! RunOnePoint(rMax) ) return false;
850  clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
851  if (clsMin.first == 1 || clsMin.first - 3 * fabs(clsMin.second) > clsTarget) break;
852  rMin += rMin;
853  if (tries == 5) {
854  oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Cannot set lower limit: at " << r->GetName()
855  << " = " << rMin << " still get " << (fUseCLs ? "CLs" : "CLsplusb")
856  << " = " << clsMin.first << std::endl;
857  return false;
858  }
859  }
860  }
861  }
862 
863  if (fVerbose > 0)
864  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Now doing proper bracketing & bisection" << std::endl;
865  do {
866 
867  // break loop in case max toys is reached
868  if (fMaxToys > 0 && fTotalToysRun > fMaxToys ) {
869  oocoutW((TObject*)0,Eval) << "HypoTestInverter::RunLimit - maximum number of toys reached " << std::endl;
870  done = false; break;
871  }
872 
873 
874  // determine point by bisection or interpolation
875  limit = 0.5*(rMin+rMax); limitErr = 0.5*(rMax-rMin);
876  if (fgAlgo == "logSecant" && clsMax.first != 0) {
877  double logMin = log(clsMin.first), logMax = log(clsMax.first), logTarget = log(clsTarget);
878  limit = rMin + (rMax-rMin) * (logTarget - logMin)/(logMax - logMin);
879  if (clsMax.second != 0 && clsMin.second != 0) {
880  limitErr = hypot((logTarget-logMax) * (clsMin.second/clsMin.first), (logTarget-logMin) * (clsMax.second/clsMax.first));
881  limitErr *= (rMax-rMin)/((logMax-logMin)*(logMax-logMin));
882  }
883  }
884  r->setError(limitErr);
885 
886  // exit if reached accuracy on r
887  if (limitErr < std::max(absAccuracy, relAccuracy * limit)) {
888  if (fVerbose > 1)
889  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - reached accuracy " << limitErr
890  << " below " << std::max(absAccuracy, relAccuracy * limit) << std::endl;
891  done = true; break;
892  }
893 
894  // evaluate point
895 
896  //clsMid = eval(w, mc_s, mc_b, data, limit, true, clsTarget);
897  if (! RunOnePoint(limit, true, clsTarget) ) return false;
898  clsMid = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
899 
900  if (clsMid.second == -1) {
901  std::cerr << "Hypotest failed" << std::endl;
902  return false;
903  }
904 
905  // if sufficiently far away, drop one of the points
906  if (fabs(clsMid.first-clsTarget) >= 2*clsMid.second) {
907  if ((clsMid.first>clsTarget) == (clsMax.first>clsTarget)) {
908  rMax = limit; clsMax = clsMid;
909  } else {
910  rMin = limit; clsMin = clsMid;
911  }
912  } else {
913  if (fVerbose > 0) std::cout << "Trying to move the interval edges closer" << std::endl;
914  double rMinBound = rMin, rMaxBound = rMax;
915  // try to reduce the size of the interval
916  while (clsMin.second == 0 || fabs(rMin-limit) > std::max(absAccuracy, relAccuracy * limit)) {
917  rMin = 0.5*(rMin+limit);
918  if (!RunOnePoint(rMin,true, clsTarget) ) return false;
919  clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
920  //clsMin = eval(w, mc_s, mc_b, data, rMin, true, clsTarget);
921  if (fabs(clsMin.first-clsTarget) <= 2*clsMin.second) break;
922  rMinBound = rMin;
923  }
924  while (clsMax.second == 0 || fabs(rMax-limit) > std::max(absAccuracy, relAccuracy * limit)) {
925  rMax = 0.5*(rMax+limit);
926 // clsMax = eval(w, mc_s, mc_b, data, rMax, true, clsTarget);
927  if (!RunOnePoint(rMax,true,clsTarget) ) return false;
928  clsMax = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
929  if (fabs(clsMax.first-clsTarget) <= 2*clsMax.second) break;
930  rMaxBound = rMax;
931  }
932  expoFit.SetRange(rMinBound,rMaxBound);
933  break;
934  }
935  } while (true);
936 
937  if (!done) { // didn't reach accuracy with scan, now do fit
938  if (fVerbose) {
939  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Before fit --- \n";
940  std::cout << "Limit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " [" << rMin << ", " << rMax << "]\n";
941  }
942 
943  expoFit.FixParameter(0,clsTarget);
944  expoFit.SetParameter(1,log(clsMax.first/clsMin.first)/(rMax-rMin));
945  expoFit.SetParameter(2,limit);
946  double rMinBound, rMaxBound; expoFit.GetRange(rMinBound, rMaxBound);
947  limitErr = std::max(fabs(rMinBound-limit), fabs(rMaxBound-limit));
948  int npoints = 0;
949 
950  HypoTestInverterPlot plot("plot","plot",fResults);
951  fLimitPlot.reset(plot.MakePlot() );
952 
953 
954  for (int j = 0; j < fLimitPlot->GetN(); ++j) {
955  if (fLimitPlot->GetX()[j] >= rMinBound && fLimitPlot->GetX()[j] <= rMaxBound) npoints++;
956  }
957  for (int i = 0, imax = /*(readHybridResults_ ? 0 : */ 8; i <= imax; ++i, ++npoints) {
958  fLimitPlot->Sort();
959  fLimitPlot->Fit(&expoFit,(fVerbose <= 1 ? "QNR EX0" : "NR EXO"));
960  if (fVerbose) {
961  oocoutI((TObject*)0,Eval) << "Fit to " << npoints << " points: " << expoFit.GetParameter(2) << " +/- " << expoFit.GetParError(2) << std::endl;
962  }
963  if ((rMin < expoFit.GetParameter(2)) && (expoFit.GetParameter(2) < rMax) && (expoFit.GetParError(2) < 0.5*(rMaxBound-rMinBound))) {
964  // sanity check fit result
965  limit = expoFit.GetParameter(2);
966  limitErr = expoFit.GetParError(2);
967  if (limitErr < std::max(absAccuracy, relAccuracy * limit)) break;
968  }
969  // add one point in the interval.
970  double rTry = RooRandom::uniform()*(rMaxBound-rMinBound)+rMinBound;
971  if (i != imax) {
972  if (!RunOnePoint(rTry,true,clsTarget) ) return false;
973  //eval(w, mc_s, mc_b, data, rTry, true, clsTarget);
974  }
975 
976  }
977  }
978 
979 //if (!plot_.empty() && fLimitPlot.get()) {
980  if (fLimitPlot.get() && fLimitPlot->GetN() > 0) {
981  //new TCanvas("c1","c1");
982  fLimitPlot->Sort();
983  fLimitPlot->SetLineWidth(2);
984  double xmin = r->getMin(), xmax = r->getMax();
985  for (int j = 0; j < fLimitPlot->GetN(); ++j) {
986  if (fLimitPlot->GetY()[j] > 1.4*clsTarget || fLimitPlot->GetY()[j] < 0.6*clsTarget) continue;
987  xmin = std::min(fLimitPlot->GetX()[j], xmin);
988  xmax = std::max(fLimitPlot->GetX()[j], xmax);
989  }
990  fLimitPlot->GetXaxis()->SetRangeUser(xmin,xmax);
991  fLimitPlot->GetYaxis()->SetRangeUser(0.5*clsTarget, 1.5*clsTarget);
992  fLimitPlot->Draw("AP");
993  expoFit.Draw("SAME");
994  TLine line(fLimitPlot->GetX()[0], clsTarget, fLimitPlot->GetX()[fLimitPlot->GetN()-1], clsTarget);
996  line.DrawLine(limit, 0, limit, fLimitPlot->GetY()[0]);
998  line.DrawLine(limit-limitErr, 0, limit-limitErr, fLimitPlot->GetY()[0]);
999  line.DrawLine(limit+limitErr, 0, limit+limitErr, fLimitPlot->GetY()[0]);
1000  //c1->Print(plot_.c_str());
1001  }
1002 
1003  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Result: \n"
1004  << "\tLimit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " @ " << (1-fSize) * 100 << "% CL\n";
1005  if (fVerbose > 1) oocoutI((TObject*)0,Eval) << "Total toys: " << fTotalToysRun << std::endl;
1006 
1007  // set value in results
1008  fResults->fUpperLimit = limit;
1009  fResults->fUpperLimitError = limitErr;
1010  fResults->fFittedUpperLimit = true;
1011  // lower limit are always min of p value
1014  fResults->fFittedLowerLimit = true;
1015 
1016  return true;
1017 }
1018 
1019 ////////////////////////////////////////////////////////////////////////////////
1020 /// get the distribution of lower limit
1021 /// if rebuild = false (default) it will re-use the results of the scan done
1022 /// for obtained the observed limit and no extra toys will be generated
1023 /// if rebuild a new set of B toys will be done and the procedure will be repeated
1024 /// for each toy
1025 
1027  if (!rebuild) {
1028  if (!fResults) {
1029  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::GetLowerLimitDistribution(false) - result not existing\n";
1030  return 0;
1031  }
1033  }
1034 
1035  TList * clsDist = 0;
1036  TList * clsbDist = 0;
1037  if (fUseCLs) clsDist = &fResults->fExpPValues;
1038  else clsbDist = &fResults->fExpPValues;
1039 
1040  return RebuildDistributions(false, nToys,clsDist, clsbDist);
1041 
1042 }
1043 
1044 ////////////////////////////////////////////////////////////////////////////////
1045 /// get the distribution of lower limit
1046 /// if rebuild = false (default) it will re-use the results of the scan done
1047 /// for obtained the observed limit and no extra toys will be generated
1048 /// if rebuild a new set of B toys will be done and the procedure will be repeated
1049 /// for each toy
1050 /// The nuisance parameter value used for rebuild is the current one in the model
1051 /// so it is user responsibility to set to the desired value (nomi
1052 
1054  if (!rebuild) {
1055  if (!fResults) {
1056  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::GetUpperLimitDistribution(false) - result not existing\n";
1057  return 0;
1058  }
1060  }
1061 
1062  TList * clsDist = 0;
1063  TList * clsbDist = 0;
1064  if (fUseCLs) clsDist = &fResults->fExpPValues;
1065  else clsbDist = &fResults->fExpPValues;
1066 
1067  return RebuildDistributions(true, nToys,clsDist, clsbDist);
1068 }
1069 
1070 ////////////////////////////////////////////////////////////////////////////////
1071 
1073  if (fCalculator0) fCalculator0->SetData(data);
1074 }
1075 
1076 ////////////////////////////////////////////////////////////////////////////////
1077 /// rebuild the sampling distributions by
1078 /// generating some toys and find for each of them a new upper limit
1079 /// Return the upper limit distribution and optionally also the pValue distributions for Cls, Clsb and Clbxs
1080 /// as a TList for each scanned point
1081 /// The method uses the present parameter value. It is user responsibility to give the current parameters to rebuild the distributions
1082 /// It returns a upper or lower limit distribution depending on the isUpper flag, however it computes also the lower limit distribution and it is saved in the
1083 /// output file as an histogram
1084 
1085 SamplingDistribution * HypoTestInverter::RebuildDistributions(bool isUpper, int nToys, TList * clsDist, TList * clsbDist, TList * clbDist, const char *outputfile) {
1086 
1087  if (!fScannedVariable || !fCalculator0) return 0;
1088  // get first background snapshot
1089  const ModelConfig * bModel = fCalculator0->GetAlternateModel();
1090  const ModelConfig * sbModel = fCalculator0->GetNullModel();
1091  if (!bModel || ! sbModel) return 0;
1092  RooArgSet paramPoint;
1093  if (!sbModel->GetParametersOfInterest()) return 0;
1094  paramPoint.add(*sbModel->GetParametersOfInterest());
1095 
1096  const RooArgSet * poibkg = bModel->GetSnapshot();
1097  if (!poibkg) {
1098  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - background snapshot not existing"
1099  << " assume is for POI = 0" << std::endl;
1101  paramPoint = RooArgSet(*fScannedVariable);
1102  }
1103  else
1104  paramPoint = *poibkg;
1105  // generate data at bkg parameter point
1106 
1107  ToyMCSampler * toymcSampler = dynamic_cast<ToyMCSampler *>(fCalculator0->GetTestStatSampler() );
1108  if (!toymcSampler) {
1109  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - no toy MC sampler existing" << std::endl;
1110  return 0;
1111  }
1112  // set up test stat sampler in case of asymptotic calculator
1113  if (dynamic_cast<RooStats::AsymptoticCalculator*>(fCalculator0) ) {
1114  toymcSampler->SetObservables(*sbModel->GetObservables() );
1115  toymcSampler->SetParametersForTestStat(*sbModel->GetParametersOfInterest());
1116  toymcSampler->SetPdf(*sbModel->GetPdf());
1117  toymcSampler->SetNuisanceParameters(*sbModel->GetNuisanceParameters());
1118  if (sbModel->GetGlobalObservables() ) toymcSampler->SetGlobalObservables(*sbModel->GetGlobalObservables() );
1119  // set number of events
1120  if (!sbModel->GetPdf()->canBeExtended())
1121  toymcSampler->SetNEventsPerToy(1);
1122  }
1123 
1124  // loop on data to generate
1125  int nPoints = fNBins;
1126 
1127  bool storePValues = clsDist || clsbDist || clbDist;
1128  if (fNBins <=0 && storePValues) {
1129  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - cannot return p values distribution with the auto scan" << std::endl;
1130  storePValues = false;
1131  nPoints = 0;
1132  }
1133 
1134  if (storePValues) {
1135  if (fResults) nPoints = fResults->ArraySize();
1136  if (nPoints <=0) {
1137  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - result is not existing and number of point to scan is not set"
1138  << std::endl;
1139  return 0;
1140  }
1141  }
1142 
1143  if (nToys <= 0) nToys = 100; // default value
1144 
1145  std::vector<std::vector<double> > CLs_values(nPoints);
1146  std::vector<std::vector<double> > CLsb_values(nPoints);
1147  std::vector<std::vector<double> > CLb_values(nPoints);
1148 
1149  if (storePValues) {
1150  for (int i = 0; i < nPoints; ++i) {
1151  CLs_values[i].reserve(nToys);
1152  CLb_values[i].reserve(nToys);
1153  CLsb_values[i].reserve(nToys);
1154  }
1155  }
1156 
1157  std::vector<double> limit_values; limit_values.reserve(nToys);
1158 
1159  oocoutI((TObject*)0,InputArguments) << "HypoTestInverter - rebuilding the p value distributions by generating ntoys = "
1160  << nToys << std::endl;
1161 
1162 
1163  oocoutI((TObject*)0,InputArguments) << "Rebuilding using parameter of interest point: ";
1165  if (sbModel->GetNuisanceParameters() ) {
1166  oocoutI((TObject*)0,InputArguments) << "And using nuisance parameters: ";
1168  }
1169  // save all parameters to restore them later
1170  assert(bModel->GetPdf() );
1171  assert(bModel->GetObservables() );
1172  RooArgSet * allParams = bModel->GetPdf()->getParameters( *bModel->GetObservables() );
1173  RooArgSet saveParams;
1174  allParams->snapshot(saveParams);
1175 
1176  TFile * fileOut = TFile::Open(outputfile,"RECREATE");
1177  if (!fileOut) {
1178  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - RebuildDistributions - Error opening file " << outputfile
1179  << " - the resulting limits will not be stored" << std::endl;
1180  }
1181  // create temporary histograms to store the limit result
1182  TH1D * hL = new TH1D("lowerLimitDist","Rebuilt lower limit distribution",100,1.,0.);
1183  TH1D * hU = new TH1D("upperLimitDist","Rebuilt upper limit distribution",100,1.,0.);
1184  TH1D * hN = new TH1D("nObs","Observed events",100,1.,0.);
1185  hL->SetBuffer(2*nToys);
1186  hU->SetBuffer(2*nToys);
1187  std::vector<TH1*> hCLb;
1188  std::vector<TH1*> hCLsb;
1189  std::vector<TH1*> hCLs;
1190  if (storePValues) {
1191  for (int i = 0; i < nPoints; ++i) {
1192  hCLb.push_back(new TH1D(TString::Format("CLbDist_bin%d",i),"CLb distribution",100,1.,0.));
1193  hCLs.push_back(new TH1D(TString::Format("ClsDist_bin%d",i),"CLs distribution",100,1.,0.));
1194  hCLsb.push_back(new TH1D(TString::Format("CLsbDist_bin%d",i),"CLs+b distribution",100,1.,0.));
1195  }
1196  }
1197 
1198 
1199  // loop now on the toys
1200  for (int itoy = 0; itoy < nToys; ++itoy) {
1201 
1202  oocoutP((TObject*)0,Eval) << "\nHypoTestInverter - RebuildDistributions - running toy # " << itoy << " / "
1203  << nToys << std::endl;
1204 
1205 
1206  printf("\n\nshnapshot of s+b model \n");
1207  sbModel->GetSnapshot()->Print("v");
1208 
1209  // reset parameters to initial values to be sure in case they are not reset
1210  if (itoy> 0) *allParams = saveParams;
1211 
1212  // need to set th epdf to clear the cache in ToyMCSampler
1213  // pdf we must use is background pdf
1214  toymcSampler->SetPdf(*bModel->GetPdf() );
1215 
1216 
1217  RooAbsData * bkgdata = toymcSampler->GenerateToyData(paramPoint);
1218 
1219  double nObs = bkgdata->sumEntries();
1220  // for debugging in case of number counting models
1221  if (bkgdata->numEntries() ==1 && !bModel->GetPdf()->canBeExtended()) {
1222  oocoutP((TObject*)0,Generation) << "Generate observables are : ";
1223  RooArgList genObs(*bkgdata->get(0));
1225  nObs = 0;
1226  for (int i = 0; i < genObs.getSize(); ++i) {
1227  RooRealVar * x = dynamic_cast<RooRealVar*>(&genObs[i]);
1228  if (x) nObs += x->getVal();
1229  }
1230  }
1231  hN->Fill(nObs);
1232 
1233  // by copying I will have the same min/max as previous ones
1234  HypoTestInverter inverter = *this;
1235  inverter.SetData(*bkgdata);
1236 
1237  // print global observables
1238  auto gobs = bModel->GetPdf()->getVariables()->selectCommon(* sbModel->GetGlobalObservables() );
1239  gobs->Print("v");
1240 
1241  HypoTestInverterResult * r = inverter.GetInterval();
1242 
1243  if (r == 0) continue;
1244 
1245  double value = (isUpper) ? r->UpperLimit() : r->LowerLimit();
1246  limit_values.push_back( value );
1247  hU->Fill(r->UpperLimit() );
1248  hL->Fill(r->LowerLimit() );
1249 
1250 
1251  std::cout << "The computed upper limit for toy #" << itoy << " is " << value << std::endl;
1252 
1253  // write every 10 toys
1254  if (itoy%10 == 0 || itoy == nToys-1) {
1255  hU->Write("",TObject::kOverwrite);
1256  hL->Write("",TObject::kOverwrite);
1257  hN->Write("",TObject::kOverwrite);
1258  }
1259 
1260  if (!storePValues) continue;
1261 
1262  if (nPoints < r->ArraySize()) {
1263  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: skip extra points" << std::endl;
1264  }
1265  else if (nPoints > r->ArraySize()) {
1266  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: missing some points" << std::endl;
1267  }
1268 
1269 
1270  for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1271  HypoTestResult * hr = r->GetResult(ipoint);
1272  if (hr) {
1273  CLs_values[ipoint].push_back( hr->CLs() );
1274  CLsb_values[ipoint].push_back( hr->CLsplusb() );
1275  CLb_values[ipoint].push_back( hr->CLb() );
1276  hCLs[ipoint]->Fill( hr->CLs() );
1277  hCLb[ipoint]->Fill( hr->CLb() );
1278  hCLsb[ipoint]->Fill( hr->CLsplusb() );
1279  }
1280  else {
1281  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: missing result for point: x = "
1282  << fResults->GetXValue(ipoint) << std::endl;
1283  }
1284  }
1285  // write every 10 toys
1286  if (itoy%10 == 0 || itoy == nToys-1) {
1287  for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1288  hCLs[ipoint]->Write("",TObject::kOverwrite);
1289  hCLb[ipoint]->Write("",TObject::kOverwrite);
1290  hCLsb[ipoint]->Write("",TObject::kOverwrite);
1291  }
1292  }
1293 
1294 
1295  delete r;
1296  delete bkgdata;
1297  }
1298 
1299 
1300  if (storePValues) {
1301  if (clsDist) clsDist->SetOwner(true);
1302  if (clbDist) clbDist->SetOwner(true);
1303  if (clsbDist) clsbDist->SetOwner(true);
1304 
1305  oocoutI((TObject*)0,InputArguments) << "HypoTestInverter: storing rebuilt p values " << std::endl;
1306 
1307  for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1308  if (clsDist) {
1309  TString name = TString::Format("CLs_distrib_%d",ipoint);
1310  clsDist->Add( new SamplingDistribution(name,name,CLs_values[ipoint] ) );
1311  }
1312  if (clbDist) {
1313  TString name = TString::Format("CLb_distrib_%d",ipoint);
1314  clbDist->Add( new SamplingDistribution(name,name,CLb_values[ipoint] ) );
1315  }
1316  if (clsbDist) {
1317  TString name = TString::Format("CLsb_distrib_%d",ipoint);
1318  clsbDist->Add( new SamplingDistribution(name,name,CLsb_values[ipoint] ) );
1319  }
1320  }
1321  }
1322 
1323  if (fileOut) {
1324  fileOut->Close();
1325  }
1326  else {
1327  // delete all the histograms
1328  delete hL;
1329  delete hU;
1330  for (int i = 0; i < nPoints && storePValues; ++i) {
1331  delete hCLs[i];
1332  delete hCLb[i];
1333  delete hCLsb[i];
1334  }
1335  }
1336 
1337  const char * disName = (isUpper) ? "upperLimit_dist" : "lowerLimit_dist";
1338  return new SamplingDistribution(disName, disName, limit_values);
1339 }
virtual Double_t sumEntries() const =0
virtual ~HypoTestInverter()
destructor (delete the HypoTestInverterResult)
void CreateResults() const
create a new HypoTestInverterResult to hold all computed results
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
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:778
bool RunOnePoint(double thisX, bool adaptive=false, double clTarget=-1) const
run only one point at the given POI value
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
void SetBackgroundAsAlt(Bool_t l=kTRUE)
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3126
float xmin
Definition: THbookFile.cxx:93
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
const ModelConfig * GetAlternateModel(void) const
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:237
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
HypoTestInverter class for performing an hypothesis test inversion by scanning the hypothesis test re...
IntervalCalculator is an interface class for a tools which produce RooStats ConfIntervals.
virtual void SetParametersForTestStat(const RooArgSet &nullpoi)
Definition: ToyMCSampler.h:155
TLine * line
HypoTestResult * Eval(HypoTestCalculatorGeneric &hc, bool adaptive, double clsTarget) const
Run the Hypothesis test at a previous configured point (internal function called by RunOnePoint) ...
virtual const RooArgSet * get() const
Definition: RooAbsData.h:77
virtual Double_t CLb() const
Convert NullPValue into a "confidence level".
virtual TestStatistic * GetTestStatistic() const =0
virtual void SetGlobalObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:175
RooAbsCollection * selectCommon(const RooAbsCollection &refColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
Definition: Rtypes.h:56
bool RunFixedScan(int nBins, double xMin, double xMax, bool scanLog=false) const
Run a Fixed scan in npoints between min and max.
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
#define oocoutF(o, a)
Definition: RooMsgService.h:48
std::unique_ptr< TGraphErrors > fLimitPlot
TestStatSampler is an interface class for a tools which produce RooStats SamplingDistributions.
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:131
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
TH1 * h
Definition: legend2.C:5
#define oocoutI(o, a)
Definition: RooMsgService.h:44
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results ...
HypoTestResult is a base class for results from hypothesis tests.
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF1.cxx:3288
virtual TLine * DrawLine(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Draw this line with new coordinates.
Definition: TLine.cxx:93
SamplingDistribution * GetUpperLimitDistribution(bool rebuild=false, int nToys=100)
get the distribution of lower limit if rebuild = false (default) it will re-use the results of the sc...
virtual void SetNEventsPerToy(const Int_t nevents)
Definition: ToyMCSampler.h:147
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
static void CheckInputModels(const HypoTestCalculatorGeneric &hc, const RooRealVar &scanVar)
check the model given the given hypotestcalculator
static unsigned int fgNToys
bool Bool_t
Definition: RtypesCore.h:59
SamplingDistribution * GetAltDistribution(void) const
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1087
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:202
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition: TF1.cxx:1656
STL namespace.
Double_t CLbError() const
The error on the "confidence level" of the null hypothesis.
HypoTestInverterResult * fResults
void setMax(const char *name, Double_t value)
Set maximum of name range to given value.
Definition: RooRealVar.cxx:418
virtual void SetBuffer(Int_t buffersize, Option_t *option="")
Set the maximum number of entries to be kept in the buffer.
Definition: TH1.cxx:7582
overwrite existing object with same name
Definition: TObject.h:79
HypoTestInverter & operator=(const HypoTestInverter &rhs)
assignment operator NOTE: this class does not copy the contained result and the HypoTestCalculator, but only the pointers It requires the original HTI to be alive
virtual void SetTestStatistic(TestStatistic *testStatistic)=0
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=1, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3909
Common base class for the Hypothesis Test Calculators.
#define oocoutP(o, a)
Definition: RooMsgService.h:45
Double_t x[n]
Definition: legend1.C:17
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:2345
Double_t CLsError() const
The error on the ratio .
#define oocoutE(o, a)
Definition: RooMsgService.h:47
Bool_t GetPValueIsRightTail(void) const
SamplingDistribution * GetLowerLimitDistribution() const
get expected lower limit distributions implemented using interpolation The size for the sampling dist...
void UseCLs(bool on=true)
flag to switch between using CLsb (default) or CLs as confidence level
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
Double_t LowerLimit()
lower and upper bound of the confidence interval (to get upper/lower limits, multiply the size( = 1-c...
static void SetCloseProof(Bool_t flag)
set flag to close proof for every new run
virtual void SetData(RooAbsData &data)
HypoTestInverter()
default constructor (doesn&#39;t do anything)
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
A doubly linked list.
Definition: TList.h:43
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition: TMath.h:318
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:322
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual HypoTestResult * GetHypoTest() const
inherited methods from HypoTestCalculator interface
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
TestStatistic * GetTestStatistic() const
return the test statistic which is or will be used by the class
TGraphErrors * MakePlot(Option_t *opt="")
return a TGraphErrors with the obtained observed p-values resultinf from the scan By default (Option ...
TRandom2 r(17)
RooAbsArg * first() const
void setMin(const char *name, Double_t value)
Set minimum of name range to given value.
Definition: RooRealVar.cxx:388
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:679
Double_t GetTestStatisticData(void) const
int fTotalToysRun
plot of limits
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
virtual void FixParameter(Int_t ipar, Double_t value)
Fix the value of a parameter The specified value will be used in a fit operation. ...
Definition: TF1.cxx:1289
A simple line.
Definition: TLine.h:23
TList fExpPValues
list of HypoTestResult for each point
Int_t GetSize() const
size of samples
RooRealVar * fScannedVariable
pointer to the generic hypotest calculator used
SamplingDistribution * GetNullDistribution(void) const
HypoTestCalculatorGeneric * fCalculator0
SamplingDistribution * GetUpperLimitDistribution() const
get expected upper limit distributions implemented using interpolation
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
float xmax
Definition: THbookFile.cxx:93
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:71
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:594
void SetTestStatisticData(const Double_t tsd)
This class simply holds a sampling distribution of some test statistic.
virtual void Append(const HypoTestResult *other)
add values from another HypoTestResult
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:216
Same purpose as HybridCalculatorOriginal, but different implementation.
constexpr Double_t E()
Definition: TMath.h:74
virtual Double_t CLs() const
is simply (not a method, but a quantity)
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
virtual Double_t ConfidenceLevel() const
Get the Confidence level for the test.
virtual Double_t CLsplusb() const
Convert AlternatePValue into a "confidence level".
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:222
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:225
#define ClassImp(name)
Definition: Rtypes.h:336
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:84
RooAbsArg * find(const char *name) const
Find object with given name in list.
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
virtual HypoTestInverterResult * GetInterval() const
Run a fixed scan or the automatic scan depending on the configuration Return if needed a copy of the ...
virtual void SetObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:173
Class to plot an HypoTestInverterResult, result of the HypoTestInverter calculator.
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 type
Definition: TGX11.cxx:120
virtual void SetPdf(RooAbsPdf &pdf)
Definition: ToyMCSampler.h:160
Double_t CLsplusbError() const
The error on the "confidence level" of the alternative hypothesis.
virtual RooAbsData * GenerateToyData(RooArgSet &paramPoint, RooAbsPdf &pdf) const
Definition: ToyMCSampler.h:110
const ModelConfig * GetNullModel(void) const
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TAttLine.h:42
bool SetTestStatistic(TestStatistic &stat)
set the test statistic to use
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TNamed.cxx:65
#define oocoutW(o, a)
Definition: RooMsgService.h:46
Mother of all ROOT objects.
Definition: TObject.h:37
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:243
virtual void Add(TObject *obj)
Definition: TList.h:77
double fLowerLimitError
interpolation option (linear or spline)
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:466
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:228
virtual void GetRange(Double_t *xmin, Double_t *xmax) const
Return range of a generic N-D function.
Definition: TF1.cxx:2012
1-Dim function class
Definition: TF1.h:150
Int_t IsNaN(Double_t x)
Definition: TMath.h:778
std::unique_ptr< HypoTestCalculatorGeneric > fHC
Does a frequentist hypothesis test.
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:310
void Clear()
delete contained result and graph
static RooRealVar * GetVariableToScan(const HypoTestCalculatorGeneric &hc)
get the variable to scan try first with null model if not go to alternate model
double result[121]
Py_ssize_t ArraySize(const std::string &name)
Extract size from an array type, if available.
Definition: Utility.cxx:682
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:578
bool RunLimit(double &limit, double &limitErr, double absTol=0, double relTol=0, const double *hint=0) const
run an automatic scan until the desired accuracy is reached Start by default from the full interval (...
const RooArgSet * GetSnapshot() const
get RooArgSet for parameters for a particular hypothesis (return NULL if not existing) ...
SamplingDistribution * RebuildDistributions(bool isUpper=true, int nToys=100, TList *clsDist=0, TList *clsbDist=0, TList *clbDist=0, const char *outputfile="HypoTestInverterRebuiltDist.root")
rebuild the sampling distributions by generating some toys and find for each of them a new upper limi...
double exp(double)
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:155
virtual void SetNuisanceParameters(const RooArgSet &np)
Definition: ToyMCSampler.h:171
virtual void SetData(RooAbsData &)
Set the DataSet ( add to the the workspace if not already there ?)
static void CloseProof(Option_t *option="s")
close all proof connections
Definition: ProofConfig.h:92
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
Definition: TestStatistic.h:31
int ArraySize() const
number of entries in the results array
SamplingDistribution * GetLowerLimitDistribution(bool rebuild=false, int nToys=100)
get the distribution of lower limit if rebuild = false (default) it will re-use the results of the sc...
Hypothesis Test Calculator based on the asymptotic formulae for the profile likelihood ratio...
virtual void SetConfidenceLevel(Double_t cl)
set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval) ...
double log(double)
TestStatSampler * GetTestStatSampler(void) const
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:269
void setError(Double_t value)
Definition: RooRealVar.h:55
std::vector< double > fXValues
number of points used to build expected p-values