Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
36
48
49#include "RooAbsData.h"
50#include "RooRealVar.h"
51#include "RooArgSet.h"
52#include "RooAbsPdf.h"
53#include "RooRandom.h"
54#include "RooConstVar.h"
55#include "RooMsgService.h"
56
57#include "TMath.h"
58#include "TF1.h"
59#include "TFile.h"
60#include "TH1.h"
61#include "TLine.h"
62#include "TCanvas.h"
63#include "TGraphErrors.h"
64
66
67#include <cassert>
68#include <cmath>
69#include <memory>
70
72
73using namespace RooStats;
74using namespace std;
75
76// static variable definitions
78unsigned int HypoTestInverter::fgNToys = 500;
79
82std::string HypoTestInverter::fgAlgo = "logSecant";
83
85
86// helper class to wrap the functionality of the various HypoTestCalculators
87
88template<class HypoTestType>
90
91 static void SetToys(HypoTestType * h, int toyNull, int toyAlt) { h->SetToys(toyNull,toyAlt); }
92
93};
94
95////////////////////////////////////////////////////////////////////////////////
96/// set flag to close proof for every new run
97
99 fgCloseProof = flag;
100}
101
102////////////////////////////////////////////////////////////////////////////////
103/// get the variable to scan
104/// try first with null model if not go to alternate model
105
107
108 RooRealVar * varToScan = nullptr;
109 const ModelConfig * mc = hc.GetNullModel();
110 if (mc) {
111 const RooArgSet * poi = mc->GetParametersOfInterest();
112 if (poi) varToScan = dynamic_cast<RooRealVar*> (poi->first() );
113 }
114 if (!varToScan) {
115 mc = hc.GetAlternateModel();
116 if (mc) {
117 const RooArgSet * poi = mc->GetParametersOfInterest();
118 if (poi) varToScan = dynamic_cast<RooRealVar*> (poi->first() );
119 }
120 }
121 return varToScan;
122}
123
124////////////////////////////////////////////////////////////////////////////////
125/// check the model given the given hypotestcalculator
126
128 const ModelConfig * modelSB = hc.GetNullModel();
129 const ModelConfig * modelB = hc.GetAlternateModel();
130 if (!modelSB || ! modelB)
131 oocoutF(nullptr,InputArguments) << "HypoTestInverter - model are not existing" << std::endl;
132 assert(modelSB && modelB);
133
134 oocoutI(nullptr,InputArguments) << "HypoTestInverter ---- Input models: \n"
135 << "\t\t using as S+B (null) model : "
136 << modelSB->GetName() << "\n"
137 << "\t\t using as B (alternate) model : "
138 << modelB->GetName() << "\n" << std::endl;
139
140 // check if scanVariable is included in B model pdf
141 RooAbsPdf * bPdf = modelB->GetPdf();
142 const RooArgSet * bObs = modelB->GetObservables();
143 if (!bPdf || !bObs) {
144 oocoutE(nullptr,InputArguments) << "HypoTestInverter - B model has no pdf or observables defined" << std::endl;
145 return;
146 }
147 std::unique_ptr<RooArgSet> bParams{bPdf->getParameters(*bObs)};
148 if (!bParams) {
149 oocoutE(nullptr,InputArguments) << "HypoTestInverter - pdf of B model has no parameters" << std::endl;
150 return;
151 }
152 if (bParams->find(scanVariable.GetName() ) ) {
153 const RooArgSet * poiB = modelB->GetSnapshot();
154 if (!poiB || !poiB->find(scanVariable.GetName()) ||
155 (static_cast<RooRealVar *>(poiB->find(scanVariable.GetName())))->getVal() != 0) {
156 oocoutW(nullptr, InputArguments)
157 << "HypoTestInverter - using a B model with POI " << scanVariable.GetName() << " not equal to zero "
158 << " user must check input model configurations " << endl;
159 }
160 if (poiB) delete poiB;
161 }
162}
163
164////////////////////////////////////////////////////////////////////////////////
165/// default constructor (doesn't do anything)
166
168 fTotalToysRun(0),
169 fMaxToys(0),
170 fCalculator0(nullptr),
171 fScannedVariable(nullptr),
172 fResults(nullptr),
173 fUseCLs(false),
174 fScanLog(false),
175 fSize(0),
176 fVerbose(0),
177 fCalcType(kUndefined),
178 fNBins(0), fXmin(1), fXmax(1),
179 fNumErr(0)
180{
181}
182
183////////////////////////////////////////////////////////////////////////////////
184/// Constructor from a HypoTestCalculatorGeneric
185/// The HypoTest calculator must be a FrequentistCalculator or HybridCalculator type
186/// Other type of calculators are not supported.
187/// The calculator must be created before by using the S+B model for the null and
188/// the B model for the alt
189/// If no variable to scan are given they are assumed to be the first variable
190/// from the parameter of interests of the null model
191
193 RooRealVar* scannedVariable, double size ) :
194 fTotalToysRun(0),
195 fMaxToys(0),
196 fCalculator0(nullptr),
197 fScannedVariable(scannedVariable),
198 fResults(nullptr),
199 fUseCLs(false),
200 fScanLog(false),
201 fSize(size),
202 fVerbose(0),
203 fCalcType(kUndefined),
204 fNBins(0), fXmin(1), fXmax(1),
205 fNumErr(0)
206{
207
208 if (!fScannedVariable) {
210 }
211 if (!fScannedVariable) {
212 oocoutE(nullptr,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
213 } else {
215 }
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(nullptr,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(nullptr),
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(nullptr,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(nullptr),
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(nullptr,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
300 } else {
302 }
303}
304
305////////////////////////////////////////////////////////////////////////////////
306/// Constructor from a reference to a AsymptoticCalculator
307/// The calculator must be created before by using the S+B model for the null and
308/// the B model for the alt
309/// If no variable to scan are given they are assumed to be the first variable
310/// from the parameter of interests of the null model
311
313 RooRealVar* scannedVariable, double size ) :
314 fTotalToysRun(0),
315 fMaxToys(0),
316 fCalculator0(&hc),
317 fScannedVariable(scannedVariable),
318 fResults(nullptr),
319 fUseCLs(false),
320 fScanLog(false),
321 fSize(size),
322 fVerbose(0),
323 fCalcType(kAsymptotic),
324 fNBins(0), fXmin(1), fXmax(1),
325 fNumErr(0)
326{
327
328 if (!fScannedVariable) {
330 }
331 if (!fScannedVariable) {
332 oocoutE(nullptr,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
333 } else {
335 }
336}
337
338////////////////////////////////////////////////////////////////////////////////
339/// Constructor from a model for B model and a model for S+B.
340/// An HypoTestCalculator (Hybrid of Frequentis) will be created using the
341/// S+B model as the null and the B model as the alternate
342/// If no variable to scan are given they are assumed to be the first variable
343/// from the parameter of interests of the null model
344
346 RooRealVar * scannedVariable, ECalculatorType type, double size) :
347 fTotalToysRun(0),
348 fMaxToys(0),
349 fCalculator0(nullptr),
350 fScannedVariable(scannedVariable),
351 fResults(nullptr),
352 fUseCLs(false),
353 fScanLog(false),
354 fSize(size),
355 fVerbose(0),
356 fCalcType(type),
357 fNBins(0), fXmin(1), fXmax(1),
358 fNumErr(0)
359{
360 if(fCalcType==kFrequentist) fHC = std::make_unique<FrequentistCalculator>(data, bModel, sbModel);
361 if(fCalcType==kHybrid) fHC = std::make_unique<HybridCalculator>(data, bModel, sbModel);
362 if(fCalcType==kAsymptotic) fHC = std::make_unique<AsymptoticCalculator>(data, bModel, sbModel);
363 fCalculator0 = fHC.get();
364 // get scanned variable
365 if (!fScannedVariable) {
367 }
368 if (!fScannedVariable) {
369 oocoutE(nullptr,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
370 } else {
372 }
373}
374
375////////////////////////////////////////////////////////////////////////////////
376/// copy-constructor
377/// NOTE: this class does not copy the contained result and
378/// the HypoTestCalculator, but only the pointers
379/// It requires the original HTI to be alive
380
383 fTotalToysRun(0),
384 fCalculator0(nullptr), fScannedVariable(nullptr), // add these for Coverity
385 fResults(nullptr)
386{
387 (*this) = rhs;
388}
389
390////////////////////////////////////////////////////////////////////////////////
391/// assignment operator
392/// NOTE: this class does not copy the contained result and
393/// the HypoTestCalculator, but only the pointers
394/// It requires the original HTI to be alive
395
397 if (this == &rhs) return *this;
398 fTotalToysRun = 0;
399 fMaxToys = rhs.fMaxToys;
402 fUseCLs = rhs.fUseCLs;
403 fScanLog = rhs.fScanLog;
404 fSize = rhs.fSize;
405 fVerbose = rhs.fVerbose;
406 fCalcType = rhs.fCalcType;
407 fNBins = rhs.fNBins;
408 fXmin = rhs.fXmin;
409 fXmax = rhs.fXmax;
410 fNumErr = rhs.fNumErr;
411
412 return *this;
413}
414
415////////////////////////////////////////////////////////////////////////////////
416/// destructor (delete the HypoTestInverterResult)
417
419{
420 if (fResults) delete fResults;
421 fCalculator0 = nullptr;
422}
423
424////////////////////////////////////////////////////////////////////////////////
425/// return the test statistic which is or will be used by the class
426
428{
431 }
432 else
433 return nullptr;
434}
435
436////////////////////////////////////////////////////////////////////////////////
437/// set the test statistic to use
438
440{
443 return true;
444 }
445 else return false;
446}
447
448////////////////////////////////////////////////////////////////////////////////
449/// delete contained result and graph
450
452 if (fResults) delete fResults;
453 fResults = nullptr;
454 fLimitPlot.reset();
455}
456
457////////////////////////////////////////////////////////////////////////////////
458/// create a new HypoTestInverterResult to hold all computed results
459
461 if (fResults == nullptr) {
462 TString results_name = "result_";
463 results_name += fScannedVariable->GetName();
465 TString title = "HypoTestInverter Result For ";
466 title += fScannedVariable->GetName();
467 fResults->SetTitle(title);
468 }
471 // check if one or two sided scan
472 if (fCalculator0) {
473 // if asymptotic calculator
475 if (ac) {
477 } else {
478 // in case of the other calculators
480 if (sampler) {
482 if (pl && pl->IsTwoSided() ) fResults->fIsTwoSided = true;
483 }
484 }
485 }
486}
487
488////////////////////////////////////////////////////////////////////////////////
489/// Run a fixed scan or the automatic scan depending on the configuration.
490/// Return if needed a copy of the result object which will be managed by the user.
491
493
494 // if having a result with at least one point return it
495 if (fResults && fResults->ArraySize() >= 1) {
496 oocoutI(nullptr,Eval) << "HypoTestInverter::GetInterval - return an already existing interval " << std::endl;
497 return static_cast<HypoTestInverterResult*>(fResults->Clone());
498 }
499
500 if (fNBins > 0) {
501 oocoutI(nullptr,Eval) << "HypoTestInverter::GetInterval - run a fixed scan" << std::endl;
502 bool ret = RunFixedScan(fNBins, fXmin, fXmax, fScanLog);
503 if (!ret)
504 oocoutE(nullptr,Eval) << "HypoTestInverter::GetInterval - error running a fixed scan " << std::endl;
505 }
506 else {
507 oocoutI(nullptr,Eval) << "HypoTestInverter::GetInterval - run an automatic scan" << std::endl;
508 double limit(0);
509 double err(0);
510 bool ret = RunLimit(limit,err);
511 if (!ret)
512 oocoutE(nullptr,Eval) << "HypoTestInverter::GetInterval - error running an auto scan " << std::endl;
513 }
514
516
517 return static_cast<HypoTestInverterResult*> (fResults->Clone());
518}
519
520////////////////////////////////////////////////////////////////////////////////
521/// Run the Hypothesis test at a previous configured point
522/// (internal function called by RunOnePoint)
523
524HypoTestResult * HypoTestInverter::Eval(HypoTestCalculatorGeneric &hc, bool adaptive, double clsTarget) const {
525 //for debug
526 // std::cout << ">>>>>>>>>>> " << std::endl;
527 // std::cout << "alternate model " << std::endl;
528 // hc.GetAlternateModel()->GetNuisanceParameters()->Print("V");
529 // hc.GetAlternateModel()->GetParametersOfInterest()->Print("V");
530 // std::cout << "Null model " << std::endl;
531 // hc.GetNullModel()->GetNuisanceParameters()->Print("V");
532 // hc.GetNullModel()->GetParametersOfInterest()->Print("V");
533 // std::cout << "<<<<<<<<<<<<<<< " << std::endl;
534
535 // run the hypothesis test
536 HypoTestResult * hcResult = hc.GetHypoTest();
537 if (hcResult == nullptr) {
538 oocoutE(nullptr,Eval) << "HypoTestInverter::Eval - HypoTest failed" << std::endl;
539 return hcResult;
540 }
541
542 // since the b model is the alt need to set the flag
543 hcResult->SetBackgroundAsAlt(true);
544
545
546 // bool flipPvalue = false;
547 // if (flipPValues)
548 // hcResult->SetPValueIsRightTail(!hcResult->GetPValueIsRightTail());
549
550 // adjust for some numerical error in discrete models and == is not anymore
551 if (hcResult->GetPValueIsRightTail()) {
552 hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()-fNumErr); // issue with < vs <= in discrete models
553 } else {
554 hcResult->SetTestStatisticData(hcResult->GetTestStatisticData() +
555 fNumErr); // issue with < vs <= in discrete models
556 }
557
558 double clsMid = (fUseCLs ? hcResult->CLs() : hcResult->CLsplusb());
559 double clsMidErr = (fUseCLs ? hcResult->CLsError() : hcResult->CLsplusbError());
560
561 //if (fVerbose) std::cout << (fUseCLs ? "\tCLs = " : "\tCLsplusb = ") << clsMid << " +/- " << clsMidErr << std::endl;
562
563 if (adaptive) {
564
567
568 while (clsMidErr >= fgCLAccuracy && (clsTarget == -1 || std::abs(clsMid-clsTarget) < 3*clsMidErr) ) {
569 std::unique_ptr<HypoTestResult> more(hc.GetHypoTest());
570
571 // if (flipPValues)
572 // more->SetPValueIsRightTail(!more->GetPValueIsRightTail());
573
574 hcResult->Append(more.get());
575 clsMid = (fUseCLs ? hcResult->CLs() : hcResult->CLsplusb());
576 clsMidErr = (fUseCLs ? hcResult->CLsError() : hcResult->CLsplusbError());
577 if (fVerbose) std::cout << (fUseCLs ? "\tCLs = " : "\tCLsplusb = ") << clsMid << " +/- " << clsMidErr << std::endl;
578 }
579
580 }
581 if (fVerbose ) {
582 oocoutP(nullptr,Eval) << "P values for " << fScannedVariable->GetName() << " = " <<
583 fScannedVariable->getVal() << "\n" <<
584 "\tCLs = " << hcResult->CLs() << " +/- " << hcResult->CLsError() << "\n" <<
585 "\tCLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << "\n" <<
586 "\tCLsplusb = " << hcResult->CLsplusb() << " +/- " << hcResult->CLsplusbError() << "\n" <<
587 std::endl;
588 }
589
591 fTotalToysRun += (hcResult->GetAltDistribution()->GetSize() + hcResult->GetNullDistribution()->GetSize());
592
593 // set sampling distribution name
594 TString nullDistName = TString::Format("%s_%s_%4.2f",hcResult->GetNullDistribution()->GetName(),
596 TString altDistName = TString::Format("%s_%s_%4.2f",hcResult->GetAltDistribution()->GetName(),
598
599 hcResult->GetNullDistribution()->SetName(nullDistName);
600 hcResult->GetAltDistribution()->SetName(altDistName);
601 }
602
603 return hcResult;
604}
605
606////////////////////////////////////////////////////////////////////////////////
607/// Run a Fixed scan in npoints between min and max
608
609bool HypoTestInverter::RunFixedScan( int nBins, double xMin, double xMax, bool scanLog ) const
610{
611
613 // interpolate the limits
616
617 // safety checks
618 if ( nBins<=0 ) {
619 oocoutE(nullptr,InputArguments) << "HypoTestInverter::RunFixedScan - Please provide nBins>0\n";
620 return false;
621 }
622 if ( nBins==1 && xMin!=xMax ) {
623 oocoutW(nullptr,InputArguments) << "HypoTestInverter::RunFixedScan - nBins==1 -> I will run for xMin (" << xMin << ")\n";
624 }
625 if ( xMin==xMax && nBins>1 ) {
626 oocoutW(nullptr,InputArguments) << "HypoTestInverter::RunFixedScan - xMin==xMax -> I will enforce nBins==1\n";
627 nBins = 1;
628 }
629 if ( xMin>xMax ) {
630 oocoutE(nullptr,InputArguments) << "HypoTestInverter::RunFixedScan - Please provide xMin ("
631 << xMin << ") smaller than xMax (" << xMax << ")\n";
632 return false;
633 }
634
635 if (xMin < fScannedVariable->getMin()) {
636 xMin = fScannedVariable->getMin();
637 oocoutW(nullptr,InputArguments) << "HypoTestInverter::RunFixedScan - xMin < lower bound, using xmin = "
638 << xMin << std::endl;
639 }
640 if (xMax > fScannedVariable->getMax()) {
641 xMax = fScannedVariable->getMax();
642 oocoutW(nullptr,InputArguments) << "HypoTestInverter::RunFixedScan - xMax > upper bound, using xmax = "
643 << xMax << std::endl;
644 }
645
646 if (xMin <= 0. && scanLog) {
647 oocoutE(nullptr, InputArguments) << "HypoTestInverter::RunFixedScan - cannot go in log steps if xMin <= 0" << std::endl;
648 return false;
649 }
650
651 double thisX = xMin;
652 for (int i=0; i<nBins; i++) {
653
654 if (i > 0) { // avoids case of nBins = 1
655 if (scanLog) {
656 thisX = exp( log(xMin) + i*(log(xMax)-log(xMin))/(nBins-1) ); // scan in log x
657 } else {
658 thisX = xMin + i * (xMax - xMin) / (nBins - 1); // linear scan in x
659 }
660 }
661
662 const bool status = RunOnePoint(thisX);
663
664 // check if failed status
665 if ( status==false ) {
666 oocoutW(nullptr,Eval) << "HypoTestInverter::RunFixedScan - The hypo test for point " << thisX << " failed. Skipping." << std::endl;
667 }
668 }
669
670 return true;
671}
672
673////////////////////////////////////////////////////////////////////////////////
674/// run only one point at the given POI value
675
676bool HypoTestInverter::RunOnePoint( double rVal, bool adaptive, double clTarget) const
677{
678
680
681 // check if rVal is in the range specified for fScannedVariable
682 if ( rVal < fScannedVariable->getMin() ) {
683 oocoutE(nullptr,InputArguments) << "HypoTestInverter::RunOnePoint - Out of range: using the lower bound "
685 << " on the scanned variable rather than " << rVal<< "\n";
686 rVal = fScannedVariable->getMin();
687 }
688 if ( rVal > fScannedVariable->getMax() ) {
689 // print a message when you have a significative difference since rval is computed
690 if (rVal > fScannedVariable->getMax() * (1. + 1.E-12)) {
691 oocoutE(nullptr, InputArguments) << "HypoTestInverter::RunOnePoint - Out of range: using the upper bound "
692 << fScannedVariable->getMax() << " on the scanned variable rather than "
693 << rVal << "\n";
694 }
695 rVal = fScannedVariable->getMax();
696 }
697
698 // save old value
699 double oldValue = fScannedVariable->getVal();
700
701 // evaluate hybrid calculator at a single point
703 // need to set value of rval in hybridcalculator
704 // assume null model is S+B and alternate is B only
705 const ModelConfig * sbModel = fCalculator0->GetNullModel();
706 RooArgSet poi; poi.add(*sbModel->GetParametersOfInterest());
707 // set poi to right values
709 const_cast<ModelConfig*>(sbModel)->SetSnapshot(poi);
710
711 if (fVerbose > 0)
712 oocoutP(nullptr,Eval) << "Running for " << fScannedVariable->GetName() << " = " << fScannedVariable->getVal() << endl;
713
714 // compute the results
715 std::unique_ptr<HypoTestResult> result( Eval(*fCalculator0,adaptive,clTarget) );
716 if (!result) {
717 oocoutE(nullptr,Eval) << "HypoTestInverter - Error running point " << fScannedVariable->GetName() << " = " <<
718 fScannedVariable->getVal() << endl;
719 return false;
720 }
721 // in case of a dummy result
722 const double nullPV = result->NullPValue();
723 const double altPV = result->AlternatePValue();
724 if (!std::isfinite(nullPV) || nullPV < 0. || nullPV > 1. || !std::isfinite(altPV) || altPV < 0. || altPV > 1.) {
725 oocoutW(nullptr,Eval) << "HypoTestInverter - Skipping invalid result for point " << fScannedVariable->GetName() << " = " <<
726 fScannedVariable->getVal() << ". null p-value=" << nullPV << ", alternate p-value=" << altPV << endl;
727 return false;
728 }
729
730 double lastXtested;
731 if ( fResults->ArraySize()!=0 ) lastXtested = fResults->GetXValue(fResults->ArraySize()-1);
732 else lastXtested = -999;
733
734 if ( (std::abs(rVal) < 1 && TMath::AreEqualAbs(rVal, lastXtested,1.E-12) ) ||
735 (std::abs(rVal) >= 1 && TMath::AreEqualRel(rVal, lastXtested,1.E-12) ) ) {
736
737 oocoutI(nullptr,Eval) << "HypoTestInverter::RunOnePoint - Merge with previous result for "
738 << fScannedVariable->GetName() << " = " << rVal << std::endl;
740 if (prevResult && prevResult->GetNullDistribution() && prevResult->GetAltDistribution()) {
741 prevResult->Append(result.get());
742 }
743 else {
744 // if it was empty we re-use it
745 oocoutI(nullptr,Eval) << "HypoTestInverter::RunOnePoint - replace previous empty result\n";
746 auto oldObj = fResults->fYObjects.Remove(prevResult);
747 delete oldObj;
748
749 fResults->fYObjects.Add(result.release());
750 }
751
752 } else {
753
754 // fill the results in the HypoTestInverterResult array
755 fResults->fXValues.push_back(rVal);
756 fResults->fYObjects.Add(result.release());
757
758 }
759
760 fScannedVariable->setVal(oldValue);
761
762 return true;
763}
764
765////////////////////////////////////////////////////////////////////////////////
766/// Run an automatic scan until the desired accuracy is reached.
767/// Start by default from the full interval (min,max) of the POI and then via bisection find the line crossing
768/// the target line.
769/// Optionally, a hint can be provided and the scan will be done closer to that value.
770/// If by bisection the desired accuracy will not be reached, a fit to the points is performed.
771/// \param[out] limit The limit.
772/// \param[out] limitErr The error of the limit.
773/// \param[in] absAccuracy Desired absolute accuracy.
774/// \param[in] relAccuracy Desired relative accuracy.
775/// \param[in] hint Hint to start from or nullptr for no hint.
776
777bool HypoTestInverter::RunLimit(double &limit, double &limitErr, double absAccuracy, double relAccuracy, const double*hint) const {
778
779
780 // routine from G. Petrucciani (from HiggsCombination CMS package)
781
783
784 if ((hint != nullptr) && (*hint > r->getMin())) {
785 r->setMax(std::min<double>(3.0 * (*hint), r->getMax()));
786 r->setMin(std::max<double>(0.3 * (*hint), r->getMin()));
787 oocoutI(nullptr,InputArguments) << "HypoTestInverter::RunLimit - Use hint value " << *hint
788 << " search in interval " << r->getMin() << " , " << r->getMax() << std::endl;
789 }
790
791 // if not specified use the default values for rel and absolute accuracy
792 if (absAccuracy <= 0) absAccuracy = fgAbsAccuracy;
793 if (relAccuracy <= 0) relAccuracy = fgRelAccuracy;
794
795 typedef std::pair<double,double> CLs_t;
796 double clsTarget = fSize;
797 CLs_t clsMin(1, 0);
798 CLs_t clsMax(0, 0);
799 CLs_t clsMid(0, 0);
800 double rMin = r->getMin();
801 double rMax = r->getMax();
802 limit = 0.5*(rMax + rMin);
803 limitErr = 0.5*(rMax - rMin);
804 bool done = false;
805
806 TF1 expoFit("expoFit","[0]*exp([1]*(x-[2]))", rMin, rMax);
807
808 fLimitPlot = std::make_unique<TGraphErrors>();
809
810 if (fVerbose > 0) std::cout << "Search for upper limit to the limit" << std::endl;
811 for (int tries = 0; tries < 6; ++tries) {
812 if (! RunOnePoint(rMax) ) {
813 oocoutE(nullptr,Eval) << "HypoTestInverter::RunLimit - Hypotest failed at upper limit of scan range: " << rMax << std::endl;
814 rMax *= 0.95;
815 continue;
816 }
817 clsMax = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
818 if (clsMax.first == 0 || clsMax.first + 3 * std::abs(clsMax.second) < clsTarget ) break;
819 rMax += rMax;
820 if (tries == 5) {
821 oocoutE(nullptr,Eval) << "HypoTestInverter::RunLimit - Cannot determine upper limit of scan range. At " << r->GetName()
822 << " = " << rMax << " still getting "
823 << (fUseCLs ? "CLs" : "CLsplusb") << " = " << clsMax.first << std::endl;
824 return false;
825 }
826 }
827 if (fVerbose > 0) {
828 oocoutI(nullptr,Eval) << "HypoTestInverter::RunLimit - Search for lower limit to the limit" << std::endl;
829 }
830
831 if ( fUseCLs && rMin == 0 ) {
832 clsMin = CLs_t(1,0);
833 }
834 else {
835 if (! RunOnePoint(rMin) ) {
836 oocoutE(nullptr,Eval) << "HypoTestInverter::RunLimit - Hypotest failed at lower limit of scan range: " << rMin << std::endl;
837 return false;
838 }
839 clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
840 }
841 if (clsMin.first != 1 && clsMin.first - 3 * std::abs(clsMin.second) < clsTarget) {
842 if (fUseCLs) {
843 rMin = 0;
844 clsMin = CLs_t(1,0); // this is always true for CLs
845 } else {
846 rMin = -rMax / 4;
847 for (int tries = 0; tries < 6; ++tries) {
848 if (! RunOnePoint(rMin) ) {
849 oocoutE(nullptr,Eval) << "HypoTestInverter::RunLimit - Hypotest failed at lower limit of scan range: " << rMin << std::endl;
850 rMin = rMin == 0. ? 0.1 : rMin * 1.1;
851 continue;
852 }
853 clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
854 if (clsMin.first == 1 || clsMin.first - 3 * std::abs(clsMin.second) > clsTarget) break;
855 rMin += rMin;
856 if (tries == 5) {
857 oocoutE(nullptr,Eval) << "HypoTestInverter::RunLimit - Cannot determine lower limit of scan range. At " << r->GetName()
858 << " = " << rMin << " still get " << (fUseCLs ? "CLs" : "CLsplusb")
859 << " = " << clsMin.first << std::endl;
860 return false;
861 }
862 }
863 }
864 }
865
866 if (fVerbose > 0)
867 oocoutI(nullptr,Eval) << "HypoTestInverter::RunLimit - Now doing proper bracketing & bisection" << std::endl;
868 do {
869
870 // break loop in case max toys is reached
871 if (fMaxToys > 0 && fTotalToysRun > fMaxToys ) {
872 oocoutW(nullptr,Eval) << "HypoTestInverter::RunLimit - maximum number of toys reached " << std::endl;
873 done = false; break;
874 }
875
876
877 // determine point by bisection or interpolation
878 limit = 0.5*(rMin+rMax); limitErr = 0.5*(rMax-rMin);
879 if (fgAlgo == "logSecant" && clsMax.first != 0) {
880 double logMin = log(clsMin.first);
881 double logMax = log(clsMax.first);
882 double logTarget = log(clsTarget);
883 limit = rMin + (rMax-rMin) * (logTarget - logMin)/(logMax - logMin);
884 if (clsMax.second != 0 && clsMin.second != 0) {
885 limitErr = hypot((logTarget-logMax) * (clsMin.second/clsMin.first), (logTarget-logMin) * (clsMax.second/clsMax.first));
886 limitErr *= (rMax-rMin)/((logMax-logMin)*(logMax-logMin));
887 }
888 }
889 r->setError(limitErr);
890
891 // exit if reached accuracy on r
892 if (limitErr < std::max(absAccuracy, relAccuracy * limit)) {
893 if (fVerbose > 1) {
894 oocoutI(nullptr, Eval) << "HypoTestInverter::RunLimit - reached accuracy " << limitErr << " below "
895 << std::max(absAccuracy, relAccuracy * limit) << std::endl;
896 }
897 done = true; break;
898 }
899
900 // evaluate point
901 if (! RunOnePoint(limit, true, clsTarget) ) {
902 oocoutE(nullptr,Eval) << "HypoTestInverter::RunLimit - Hypo test failed at x=" << limit << " when trying to find limit." << std::endl;
903 return false;
904 }
905 clsMid = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
906
907 if (clsMid.second == -1) {
908 std::cerr << "Hypotest failed" << std::endl;
909 return false;
910 }
911
912 // if sufficiently far away, drop one of the points
913 if (std::abs(clsMid.first-clsTarget) >= 2*clsMid.second) {
914 if ((clsMid.first>clsTarget) == (clsMax.first>clsTarget)) {
915 rMax = limit; clsMax = clsMid;
916 } else {
917 rMin = limit; clsMin = clsMid;
918 }
919 } else {
920 if (fVerbose > 0) std::cout << "Trying to move the interval edges closer" << std::endl;
921 double rMinBound = rMin;
922 double rMaxBound = rMax;
923 // try to reduce the size of the interval
924 while (clsMin.second == 0 || std::abs(rMin-limit) > std::max(absAccuracy, relAccuracy * limit)) {
925 rMin = 0.5*(rMin+limit);
926 if (!RunOnePoint(rMin,true, clsTarget) ) {
927 oocoutE(nullptr,Eval) << "HypoTestInverter::RunLimit - Hypo test failed at x=" << rMin << " when trying to find limit from below." << std::endl;
928 return false;
929 }
930 clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
931 if (std::abs(clsMin.first-clsTarget) <= 2*clsMin.second) break;
932 rMinBound = rMin;
933 }
934 while (clsMax.second == 0 || std::abs(rMax-limit) > std::max(absAccuracy, relAccuracy * limit)) {
935 rMax = 0.5*(rMax+limit);
936 if (!RunOnePoint(rMax,true,clsTarget) ) {
937 oocoutE(nullptr,Eval) << "HypoTestInverter::RunLimit - Hypo test failed at x=" << rMin << " when trying to find limit from above." << std::endl;
938 return false;
939 }
940 clsMax = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
941 if (std::abs(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(nullptr,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;
959 double rMaxBound;
960 expoFit.GetRange(rMinBound, rMaxBound);
961 limitErr = std::max(std::abs(rMinBound-limit), std::abs(rMaxBound-limit));
962 int npoints = 0;
963
964 HypoTestInverterPlot plot("plot","plot",fResults);
965 fLimitPlot.reset(plot.MakePlot() );
966
967
968 for (int j = 0; j < fLimitPlot->GetN(); ++j) {
969 if (fLimitPlot->GetX()[j] >= rMinBound && fLimitPlot->GetX()[j] <= rMaxBound) npoints++;
970 }
971 for (int i = 0, imax = /*(readHybridResults_ ? 0 : */ 8; i <= imax; ++i, ++npoints) {
972 fLimitPlot->Sort();
973 fLimitPlot->Fit(&expoFit,(fVerbose <= 1 ? "QNR EX0" : "NR EXO"));
974 if (fVerbose) {
975 oocoutI(nullptr,Eval) << "Fit to " << npoints << " points: " << expoFit.GetParameter(2) << " +/- " << expoFit.GetParError(2) << std::endl;
976 }
977 if ((rMin < expoFit.GetParameter(2)) && (expoFit.GetParameter(2) < rMax) && (expoFit.GetParError(2) < 0.5*(rMaxBound-rMinBound))) {
978 // sanity check fit result
979 limit = expoFit.GetParameter(2);
980 limitErr = expoFit.GetParError(2);
981 if (limitErr < std::max(absAccuracy, relAccuracy * limit)) break;
982 }
983 // add one point in the interval.
984 double rTry = RooRandom::uniform()*(rMaxBound-rMinBound)+rMinBound;
985 if (i != imax) {
986 if (!RunOnePoint(rTry,true,clsTarget) ) return false;
987 //eval(w, mc_s, mc_b, data, rTry, true, clsTarget);
988 }
989
990 }
991 }
992
993//if (!plot_.empty() && fLimitPlot.get()) {
994 if (fLimitPlot.get() && fLimitPlot->GetN() > 0) {
995 //new TCanvas("c1","c1");
996 fLimitPlot->Sort();
997 fLimitPlot->SetLineWidth(2);
998 double xmin = r->getMin();
999 double xmax = r->getMax();
1000 for (int j = 0; j < fLimitPlot->GetN(); ++j) {
1001 if (fLimitPlot->GetY()[j] > 1.4*clsTarget || fLimitPlot->GetY()[j] < 0.6*clsTarget) continue;
1002 xmin = std::min(fLimitPlot->GetX()[j], xmin);
1003 xmax = std::max(fLimitPlot->GetX()[j], xmax);
1004 }
1005 fLimitPlot->GetXaxis()->SetRangeUser(xmin,xmax);
1006 fLimitPlot->GetYaxis()->SetRangeUser(0.5*clsTarget, 1.5*clsTarget);
1007 fLimitPlot->Draw("AP");
1008 expoFit.Draw("SAME");
1009 TLine line(fLimitPlot->GetX()[0], clsTarget, fLimitPlot->GetX()[fLimitPlot->GetN()-1], clsTarget);
1011 line.DrawLine(limit, 0, limit, fLimitPlot->GetY()[0]);
1013 line.DrawLine(limit-limitErr, 0, limit-limitErr, fLimitPlot->GetY()[0]);
1014 line.DrawLine(limit+limitErr, 0, limit+limitErr, fLimitPlot->GetY()[0]);
1015 //c1->Print(plot_.c_str());
1016 }
1017
1018 oocoutI(nullptr,Eval) << "HypoTestInverter::RunLimit - Result: \n"
1019 << "\tLimit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " @ " << (1-fSize) * 100 << "% CL\n";
1020 if (fVerbose > 1) oocoutI(nullptr,Eval) << "Total toys: " << fTotalToysRun << std::endl;
1021
1022 // set value in results
1023 fResults->fUpperLimit = limit;
1024 fResults->fUpperLimitError = limitErr;
1026 // lower limit are always min of p value
1030
1031 return true;
1032}
1033
1034////////////////////////////////////////////////////////////////////////////////
1035/// get the distribution of lower limit
1036/// if rebuild = false (default) it will re-use the results of the scan done
1037/// for obtained the observed limit and no extra toys will be generated
1038/// if rebuild a new set of B toys will be done and the procedure will be repeated
1039/// for each toy
1040
1042 if (!rebuild) {
1043 if (!fResults) {
1044 oocoutE(nullptr,InputArguments) << "HypoTestInverter::GetLowerLimitDistribution(false) - result not existing\n";
1045 return nullptr;
1046 }
1048 }
1049
1050 TList * clsDist = nullptr;
1051 TList * clsbDist = nullptr;
1052 if (fUseCLs) clsDist = &fResults->fExpPValues;
1053 else clsbDist = &fResults->fExpPValues;
1054
1055 return RebuildDistributions(false, nToys,clsDist, clsbDist);
1056
1057}
1058
1059////////////////////////////////////////////////////////////////////////////////
1060/// get the distribution of lower limit
1061/// if rebuild = false (default) it will re-use the results of the scan done
1062/// for obtained the observed limit and no extra toys will be generated
1063/// if rebuild a new set of B toys will be done and the procedure will be repeated
1064/// for each toy
1065/// The nuisance parameter value used for rebuild is the current one in the model
1066/// so it is user responsibility to set to the desired value (nomi
1067
1069 if (!rebuild) {
1070 if (!fResults) {
1071 oocoutE(nullptr,InputArguments) << "HypoTestInverter::GetUpperLimitDistribution(false) - result not existing\n";
1072 return nullptr;
1073 }
1075 }
1076
1077 TList * clsDist = nullptr;
1078 TList * clsbDist = nullptr;
1079 if (fUseCLs) clsDist = &fResults->fExpPValues;
1080 else clsbDist = &fResults->fExpPValues;
1081
1082 return RebuildDistributions(true, nToys,clsDist, clsbDist);
1083}
1084
1085////////////////////////////////////////////////////////////////////////////////
1086
1089}
1090
1091////////////////////////////////////////////////////////////////////////////////
1092/// rebuild the sampling distributions by
1093/// generating some toys and find for each of them a new upper limit
1094/// Return the upper limit distribution and optionally also the pValue distributions for Cls, Clsb and Clbxs
1095/// as a TList for each scanned point
1096/// The method uses the present parameter value. It is user responsibility to give the current parameters to rebuild the distributions
1097/// 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
1098/// output file as an histogram
1099
1100SamplingDistribution * HypoTestInverter::RebuildDistributions(bool isUpper, int nToys, TList * clsDist, TList * clsbDist, TList * clbDist, const char *outputfile) {
1101
1102 if (!fScannedVariable || !fCalculator0) return nullptr;
1103 // get first background snapshot
1104 const ModelConfig * bModel = fCalculator0->GetAlternateModel();
1105 const ModelConfig * sbModel = fCalculator0->GetNullModel();
1106 if (!bModel || ! sbModel) return nullptr;
1107 RooArgSet paramPoint;
1108 if (!sbModel->GetParametersOfInterest()) return nullptr;
1109 paramPoint.add(*sbModel->GetParametersOfInterest());
1110
1111 const RooArgSet * poibkg = bModel->GetSnapshot();
1112 if (!poibkg) {
1113 oocoutW(nullptr,InputArguments) << "HypoTestInverter::RebuildDistribution - background snapshot not existing"
1114 << " assume is for POI = 0" << std::endl;
1116 paramPoint.assign(RooArgSet(*fScannedVariable));
1117 }
1118 else
1119 paramPoint.assign(*poibkg);
1120 // generate data at bkg parameter point
1121
1122 ToyMCSampler * toymcSampler = dynamic_cast<ToyMCSampler *>(fCalculator0->GetTestStatSampler() );
1123 if (!toymcSampler) {
1124 oocoutE(nullptr,InputArguments) << "HypoTestInverter::RebuildDistribution - no toy MC sampler existing" << std::endl;
1125 return nullptr;
1126 }
1127 // set up test stat sampler in case of asymptotic calculator
1128 if (dynamic_cast<RooStats::AsymptoticCalculator*>(fCalculator0) ) {
1129 toymcSampler->SetObservables(*sbModel->GetObservables() );
1130 toymcSampler->SetParametersForTestStat(*sbModel->GetParametersOfInterest());
1131 toymcSampler->SetPdf(*sbModel->GetPdf());
1132 toymcSampler->SetNuisanceParameters(*sbModel->GetNuisanceParameters());
1133 if (sbModel->GetGlobalObservables() ) toymcSampler->SetGlobalObservables(*sbModel->GetGlobalObservables() );
1134 // set number of events
1135 if (!sbModel->GetPdf()->canBeExtended())
1136 toymcSampler->SetNEventsPerToy(1);
1137 }
1138
1139 // loop on data to generate
1140 int nPoints = fNBins;
1141
1142 bool storePValues = clsDist || clsbDist || clbDist;
1143 if (fNBins <=0 && storePValues) {
1144 oocoutW(nullptr,InputArguments) << "HypoTestInverter::RebuildDistribution - cannot return p values distribution with the auto scan" << std::endl;
1145 storePValues = false;
1146 nPoints = 0;
1147 }
1148
1149 if (storePValues) {
1150 if (fResults) nPoints = fResults->ArraySize();
1151 if (nPoints <=0) {
1152 oocoutE(nullptr,InputArguments) << "HypoTestInverter - result is not existing and number of point to scan is not set"
1153 << std::endl;
1154 return nullptr;
1155 }
1156 }
1157
1158 if (nToys <= 0) nToys = 100; // default value
1159
1160 std::vector<std::vector<double> > CLs_values(nPoints);
1161 std::vector<std::vector<double> > CLsb_values(nPoints);
1162 std::vector<std::vector<double> > CLb_values(nPoints);
1163
1164 if (storePValues) {
1165 for (int i = 0; i < nPoints; ++i) {
1166 CLs_values[i].reserve(nToys);
1167 CLb_values[i].reserve(nToys);
1168 CLsb_values[i].reserve(nToys);
1169 }
1170 }
1171
1172 std::vector<double> limit_values; limit_values.reserve(nToys);
1173
1174 oocoutI(nullptr,InputArguments) << "HypoTestInverter - rebuilding the p value distributions by generating ntoys = "
1175 << nToys << std::endl;
1176
1177
1178 oocoutI(nullptr,InputArguments) << "Rebuilding using parameter of interest point: ";
1179 RooStats::PrintListContent(paramPoint, oocoutI(nullptr,InputArguments) );
1180 if (sbModel->GetNuisanceParameters() ) {
1181 oocoutI(nullptr,InputArguments) << "And using nuisance parameters: ";
1182 RooStats::PrintListContent(*sbModel->GetNuisanceParameters(), oocoutI(nullptr,InputArguments) );
1183 }
1184 // save all parameters to restore them later
1185 assert(bModel->GetPdf() );
1186 assert(bModel->GetObservables() );
1187 std::unique_ptr<RooArgSet> allParams{bModel->GetPdf()->getParameters( *bModel->GetObservables() )};
1188 RooArgSet saveParams;
1189 allParams->snapshot(saveParams);
1190
1191 std::unique_ptr<TFile> fileOut{TFile::Open(outputfile,"RECREATE")};
1192 if (!fileOut) {
1193 oocoutE(nullptr,InputArguments) << "HypoTestInverter - RebuildDistributions - Error opening file " << outputfile
1194 << " - the resulting limits will not be stored" << std::endl;
1195 }
1196 // create temporary histograms to store the limit result
1197 TH1D * hL = new TH1D("lowerLimitDist","Rebuilt lower limit distribution",100,1.,0.);
1198 TH1D * hU = new TH1D("upperLimitDist","Rebuilt upper limit distribution",100,1.,0.);
1199 TH1D * hN = new TH1D("nObs","Observed events",100,1.,0.);
1200 hL->SetBuffer(2*nToys);
1201 hU->SetBuffer(2*nToys);
1202 std::vector<TH1*> hCLb;
1203 std::vector<TH1*> hCLsb;
1204 std::vector<TH1*> hCLs;
1205 if (storePValues) {
1206 for (int i = 0; i < nPoints; ++i) {
1207 hCLb.push_back(new TH1D(TString::Format("CLbDist_bin%d",i),"CLb distribution",100,1.,0.));
1208 hCLs.push_back(new TH1D(TString::Format("ClsDist_bin%d",i),"CLs distribution",100,1.,0.));
1209 hCLsb.push_back(new TH1D(TString::Format("CLsbDist_bin%d",i),"CLs+b distribution",100,1.,0.));
1210 }
1211 }
1212
1213
1214 // loop now on the toys
1215 for (int itoy = 0; itoy < nToys; ++itoy) {
1216
1217 oocoutP(nullptr,Eval) << "\nHypoTestInverter - RebuildDistributions - running toy # " << itoy << " / "
1218 << nToys << std::endl;
1219
1220
1221 printf("\n\nshnapshot of s+b model \n");
1222 sbModel->GetSnapshot()->Print("v");
1223
1224 // reset parameters to initial values to be sure in case they are not reset
1225 if (itoy> 0) {
1226 allParams->assign(saveParams);
1227 }
1228
1229 // need to set th epdf to clear the cache in ToyMCSampler
1230 // pdf we must use is background pdf
1231 toymcSampler->SetPdf(*bModel->GetPdf() );
1232
1233
1234 RooAbsData * bkgdata = toymcSampler->GenerateToyData(paramPoint);
1235
1236 double nObs = bkgdata->sumEntries();
1237 // for debugging in case of number counting models
1238 if (bkgdata->numEntries() ==1 && !bModel->GetPdf()->canBeExtended()) {
1239 oocoutP(nullptr,Generation) << "Generate observables are : ";
1240 RooArgList genObs(*bkgdata->get(0));
1241 RooStats::PrintListContent(genObs, oocoutP(nullptr,Generation) );
1242 nObs = 0;
1243 for (std::size_t i = 0; i < genObs.size(); ++i) {
1244 RooRealVar * x = dynamic_cast<RooRealVar*>(&genObs[i]);
1245 if (x) nObs += x->getVal();
1246 }
1247 }
1248 hN->Fill(nObs);
1249
1250 // by copying I will have the same min/max as previous ones
1251 HypoTestInverter inverter = *this;
1252 inverter.SetData(*bkgdata);
1253
1254 // print global observables
1255 auto gobs = bModel->GetPdf()->getVariables()->selectCommon(* sbModel->GetGlobalObservables() );
1256 gobs->Print("v");
1257
1258 HypoTestInverterResult * r = inverter.GetInterval();
1259
1260 if (r == nullptr) continue;
1261
1262 double value = (isUpper) ? r->UpperLimit() : r->LowerLimit();
1263 limit_values.push_back( value );
1264 hU->Fill(r->UpperLimit() );
1265 hL->Fill(r->LowerLimit() );
1266
1267
1268 std::cout << "The computed upper limit for toy #" << itoy << " is " << value << std::endl;
1269
1270 // write every 10 toys
1271 if (itoy%10 == 0 || itoy == nToys-1) {
1272 hU->Write("",TObject::kOverwrite);
1273 hL->Write("",TObject::kOverwrite);
1274 hN->Write("",TObject::kOverwrite);
1275 }
1276
1277 if (!storePValues) continue;
1278
1279 if (nPoints < r->ArraySize()) {
1280 oocoutW(nullptr,InputArguments) << "HypoTestInverter: skip extra points" << std::endl;
1281 }
1282 else if (nPoints > r->ArraySize()) {
1283 oocoutW(nullptr,InputArguments) << "HypoTestInverter: missing some points" << std::endl;
1284 }
1285
1286
1287 for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1288 HypoTestResult * hr = r->GetResult(ipoint);
1289 if (hr) {
1290 CLs_values[ipoint].push_back( hr->CLs() );
1291 CLsb_values[ipoint].push_back( hr->CLsplusb() );
1292 CLb_values[ipoint].push_back( hr->CLb() );
1293 hCLs[ipoint]->Fill( hr->CLs() );
1294 hCLb[ipoint]->Fill( hr->CLb() );
1295 hCLsb[ipoint]->Fill( hr->CLsplusb() );
1296 }
1297 else {
1298 oocoutW(nullptr,InputArguments) << "HypoTestInverter: missing result for point: x = "
1299 << fResults->GetXValue(ipoint) << std::endl;
1300 }
1301 }
1302 // write every 10 toys
1303 if (itoy%10 == 0 || itoy == nToys-1) {
1304 for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1305 hCLs[ipoint]->Write("",TObject::kOverwrite);
1306 hCLb[ipoint]->Write("",TObject::kOverwrite);
1307 hCLsb[ipoint]->Write("",TObject::kOverwrite);
1308 }
1309 }
1310
1311
1312 delete r;
1313 delete bkgdata;
1314 }
1315
1316
1317 if (storePValues) {
1318 if (clsDist) clsDist->SetOwner(true);
1319 if (clbDist) clbDist->SetOwner(true);
1320 if (clsbDist) clsbDist->SetOwner(true);
1321
1322 oocoutI(nullptr,InputArguments) << "HypoTestInverter: storing rebuilt p values " << std::endl;
1323
1324 for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1325 if (clsDist) {
1326 TString name = TString::Format("CLs_distrib_%d",ipoint);
1327 clsDist->Add( new SamplingDistribution(name,name,CLs_values[ipoint] ) );
1328 }
1329 if (clbDist) {
1330 TString name = TString::Format("CLb_distrib_%d",ipoint);
1331 clbDist->Add( new SamplingDistribution(name,name,CLb_values[ipoint] ) );
1332 }
1333 if (clsbDist) {
1334 TString name = TString::Format("CLsb_distrib_%d",ipoint);
1335 clsbDist->Add( new SamplingDistribution(name,name,CLsb_values[ipoint] ) );
1336 }
1337 }
1338 }
1339
1340 else {
1341 // delete all the histograms
1342 delete hL;
1343 delete hU;
1344 for (int i = 0; i < nPoints && storePValues; ++i) {
1345 delete hCLs[i];
1346 delete hCLb[i];
1347 delete hCLsb[i];
1348 }
1349 }
1350
1351 const char * disName = (isUpper) ? "upperLimit_dist" : "lowerLimit_dist";
1352 return new SamplingDistribution(disName, disName, limit_values);
1353}
size_t fSize
#define h(i)
Definition RSha256.hxx:106
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define oocoutW(o, a)
#define oocoutE(o, a)
#define oocoutF(o, a)
#define oocoutI(o, a)
#define oocoutP(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
@ kRed
Definition Rtypes.h:66
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expression tree)
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
Storage_t::size_type size() const
RooAbsArg * first() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual double sumEntries() const =0
Return effective number of entries in dataset, i.e., sum all weights.
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
bool canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition RooAbsPdf.h:219
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:178
static double uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition RooRandom.cxx:81
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
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
void SetData(RooAbsData &data) override
Set the DataSet.
HypoTestResult * GetHypoTest() const override
inherited methods from HypoTestCalculator interface
const ModelConfig * GetAlternateModel(void) const
TestStatSampler * GetTestStatSampler(void) const
Returns instance of TestStatSampler.
Class to plot a HypoTestInverterResult, the output of the HypoTestInverter calculator.
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
int ArraySize() const
number of entries in the results array
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 expected sampling distribution for each point
bool fIsTwoSided
two sided scan (look for lower/upper limit)
TList fYObjects
list of HypoTestResult for each point
void SetConfidenceLevel(double cl) override
set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
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...
static void CheckInputModels(const HypoTestCalculatorGeneric &hc, const RooRealVar &scanVar)
check the model given the given hypotestcalculator
RooRealVar * fScannedVariable
pointer to the constrained variable
std::unique_ptr< HypoTestCalculatorGeneric > fHC
! pointer to the generic hypotest calculator used
bool RunLimit(double &limit, double &limitErr, double absTol=0, double relTol=0, const double *hint=nullptr) const
Run an automatic scan until the desired accuracy is reached.
void SetData(RooAbsData &) override
Set the DataSet ( add to the workspace if not already there ?)
HypoTestInverterResult * fResults
pointer to the result
SamplingDistribution * RebuildDistributions(bool isUpper=true, int nToys=100, TList *clsDist=nullptr, TList *clsbDist=nullptr, TList *clbDist=nullptr, const char *outputfile="HypoTestInverterRebuiltDist.root")
function to rebuild the distributions
int fMaxToys
maximum number of toys to run
HypoTestInverter()
default constructor (used only for I/O)
~HypoTestInverter() override
destructor
double ConfidenceLevel() const override
Get the Confidence level for the test.
bool SetTestStatistic(TestStatistic &stat)
set the test statistic
HypoTestInverterResult * GetInterval() const override
Run a fixed scan or the automatic scan depending on the configuration.
static void SetCloseProof(bool flag)
set flag to close proof for every new run
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
get the test statistic
std::unique_ptr< TGraphErrors > fLimitPlot
! plot of limits
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 upper/lower limit distribution
HypoTestInverter & operator=(const HypoTestInverter &rhs)
assignment
HypoTestCalculatorGeneric * fCalculator0
pointer to the calculator passed in the constructor
void CreateResults() const
create a new HypoTestInverterResult to hold all computed results
static RooRealVar * GetVariableToScan(const HypoTestCalculatorGeneric &hc)
helper functions
HypoTestResult * Eval(HypoTestCalculatorGeneric &hc, bool adaptive, double clsTarget) const
run the hybrid at a single point
HypoTestResult is a base class for results from hypothesis tests.
virtual double CLsplusb() const
Convert AlternatePValue into a "confidence level".
double GetTestStatisticData(void) const
virtual void Append(const HypoTestResult *other)
add values from another HypoTestResult
double CLsError() const
The error on the ratio .
void SetBackgroundAsAlt(bool l=true)
bool GetPValueIsRightTail(void) const
double CLbError() const
The error on the "confidence level" of the null hypothesis.
void SetTestStatisticData(const double tsd)
double CLsplusbError() const
The error on the "confidence level" of the alternative hypothesis.
SamplingDistribution * GetNullDistribution(void) const
virtual double CLs() const
is simply (not a method, but a quantity)
virtual double CLb() const
Convert NullPValue into a "confidence level".
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:35
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return nullptr if not existing)
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
const RooArgSet * GetObservables() const
get RooArgSet for observables (return nullptr if not existing)
const RooArgSet * GetSnapshot() const
get RooArgSet for parameters for a particular hypothesis (return nullptr if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
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:91
This class simply holds a sampling distribution of some test statistic.
Int_t GetSize() const
size of samples
double fUpperLimit
upper interval limit
double fLowerLimit
lower interval limit
TestStatSampler is an interface class for a tools which produce RooStats SamplingDistributions.
virtual TestStatistic * GetTestStatistic() const =0
Get the TestStatistic.
virtual void SetTestStatistic(TestStatistic *testStatistic)=0
Set the TestStatistic (want the argument to be a function of the data & parameter points.
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
ToyMCSampler is an implementation of the TestStatSampler interface.
void SetParametersForTestStat(const RooArgSet &nullpoi) override
Set the Pdf, add to the workspace if not already there.
void SetObservables(const RooArgSet &o) override
specify the observables in the dataset (needed to evaluate the test statistic)
virtual RooAbsData * GenerateToyData(RooArgSet &paramPoint, RooAbsPdf &pdf) const
generates toy data without weight
void SetNuisanceParameters(const RooArgSet &np) override
specify the nuisance parameters (eg. the rest of the parameters)
void SetPdf(RooAbsPdf &pdf) override
Set the Pdf, add to the workspace if not already there.
void SetGlobalObservables(const RooArgSet &o) override
specify the conditional observables
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
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:233
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition TF1.cxx:1932
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition TF1.cxx:3528
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition TF1.cxx:1335
virtual void GetRange(Double_t *xmin, Double_t *xmax) const
Return range of a generic N-D function.
Definition TF1.cxx:2281
virtual void SetParameter(Int_t param, Double_t value)
Definition TF1.h:659
virtual void FixParameter(Int_t ipar, Double_t value)
Fix the value of a parameter for a fit operation The specified value will be used in the fit and the ...
Definition TF1.cxx:1559
virtual Double_t GetParameter(Int_t ipar) const
Definition TF1.h:537
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:4070
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:621
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3325
virtual void SetBuffer(Int_t buffersize, Option_t *option="")
Set the maximum number of entries to be kept in the buffer.
Definition TH1.cxx:8329
Use the TLine constructor to create a simple line.
Definition TLine.h:22
virtual TLine * DrawLine(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Draw this line with new coordinates.
Definition TLine.cxx:103
A doubly linked list.
Definition TList.h:38
void Add(TObject *obj) override
Definition TList.h:81
TObject * Remove(TObject *obj) override
Remove object from the list.
Definition TList.cxx:822
TObject * Clone(const char *newname="") const override
Make a clone of an object using the Streamer facility.
Definition TNamed.cxx:74
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition TNamed.cxx:140
@ kOverwrite
overwrite existing object with same name
Definition TObject.h:92
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition TObject.cxx:880
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
Basic string class.
Definition TString.h:139
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:2356
TLine * line
Double_t x[n]
Definition legend1.C:17
Namespace for the RooStats classes.
Definition Asimov.h:19
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
useful function to print in one line the content of a set with their values
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Comparing floating points.
Definition TMath.h:426
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Comparing floating points.
Definition TMath.h:418
static void SetToys(HypoTestType *h, int toyNull, int toyAlt)