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