Logo ROOT   6.18/05
Reference Guide
HypoTestInverterResult.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
3/*************************************************************************
4 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11
12/** \class RooStats::HypoTestInverterResult
13 \ingroup Roostats
14
15HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence interval.
16Based on the RatioFinder code available in the RooStatsCms package developed by Gregory Schott and Danilo Piparo
17Ported and adapted to RooStats by Gregory Schott
18Some contributions to this class have been written by Matthias Wolf (error estimation)
19
20*/
21
22// include header file of this class
24
28#include "RooMsgService.h"
29#include "RooGlobalFunc.h"
30
31#include "TF1.h"
32#include "TGraphErrors.h"
33#include <cmath>
36#include "Math/Functor.h"
37
38#include "TCanvas.h"
39#include "TFile.h"
41
42#include <algorithm>
43
45
46using namespace RooStats;
47using namespace RooFit;
48using namespace std;
49
50
51// initialize static value
52double HypoTestInverterResult::fgAsymptoticMaxSigma = 5;
53int HypoTestInverterResult::fgAsymptoticNumPoints = 11;
54
55////////////////////////////////////////////////////////////////////////////////
56/// default constructor
57
58HypoTestInverterResult::HypoTestInverterResult(const char * name ) :
60 fUseCLs(false),
61 fIsTwoSided(false),
62 fInterpolateLowerLimit(true),
63 fInterpolateUpperLimit(true),
64 fFittedLowerLimit(false),
65 fFittedUpperLimit(false),
66 fInterpolOption(kLinear),
67 fLowerLimitError(-1),
68 fUpperLimitError(-1),
69 fCLsCleanupThreshold(0.005)
70{
73
76}
77
78////////////////////////////////////////////////////////////////////////////////
79/// copy constructor
80
82 SimpleInterval(other,name),
83 fUseCLs(other.fUseCLs),
84 fIsTwoSided(other.fIsTwoSided),
85 fInterpolateLowerLimit(other.fInterpolateLowerLimit),
86 fInterpolateUpperLimit(other.fInterpolateUpperLimit),
87 fFittedLowerLimit(other.fFittedLowerLimit),
88 fFittedUpperLimit(other.fFittedUpperLimit),
89 fInterpolOption(other.fInterpolOption),
90 fLowerLimitError(other.fLowerLimitError),
91 fUpperLimitError(other.fUpperLimitError),
92 fCLsCleanupThreshold(other.fCLsCleanupThreshold)
93{
96 int nOther = other.ArraySize();
97
98 fXValues = other.fXValues;
99 for (int i = 0; i < nOther; ++i)
100 fYObjects.Add( other.fYObjects.At(i)->Clone() );
101 for (int i = 0; i < fExpPValues.GetSize() ; ++i)
102 fExpPValues.Add( other.fExpPValues.At(i)->Clone() );
103
106}
107
108////////////////////////////////////////////////////////////////////////////////
109
112{
113 if (&other==this) {
114 return *this ;
115 }
116
118 fLowerLimit = other.fLowerLimit;
119 fUpperLimit = other.fUpperLimit;
120 fUseCLs = other.fUseCLs;
121 fIsTwoSided = other.fIsTwoSided;
130
131 int nOther = other.ArraySize();
132 fXValues = other.fXValues;
133
135 for (int i=0; i < nOther; ++i) {
136 fYObjects.Add( other.fYObjects.At(i)->Clone() );
137 }
139 for (int i=0; i < fExpPValues.GetSize() ; ++i) {
140 fExpPValues.Add( other.fExpPValues.At(i)->Clone() );
141 }
142
145
146 return *this;
147}
148
149////////////////////////////////////////////////////////////////////////////////
150/// constructor
151
153 const RooRealVar& scannedVariable,
154 double cl ) :
155 SimpleInterval(name,scannedVariable,TMath::QuietNaN(),TMath::QuietNaN(),cl),
156 fUseCLs(false),
157 fIsTwoSided(false),
158 fInterpolateLowerLimit(true),
159 fInterpolateUpperLimit(true),
160 fFittedLowerLimit(false),
161 fFittedUpperLimit(false),
162 fInterpolOption(kLinear),
163 fLowerLimitError(-1),
164 fUpperLimitError(-1),
165 fCLsCleanupThreshold(0.005)
166{
169
170 // put a cloned copy of scanned variable to set in the interval
171 // to avoid I/O problem of the Result class -
172 // make the set owning the cloned copy (use clone instead of Clone to not copying all links)
175 fParameters.addOwned(*((RooRealVar *) scannedVariable.clone(scannedVariable.GetName()) ));
176}
177
178////////////////////////////////////////////////////////////////////////////////
179/// destructor
180
182{
183 // explicitly empty the TLists - these contain pointers, not objects
186
189}
190
191////////////////////////////////////////////////////////////////////////////////
192
194{
195 const int nEntries = ArraySize();
196
197 // initialization
198 double nsig1(1.0);
199 double nsig2(2.0);
200 double p[5];
201 double q[5];
202 std::vector<double> qv;
203 qv.resize(11,-1.0);
204
205 p[0] = ROOT::Math::normal_cdf(-nsig2);
206 p[1] = ROOT::Math::normal_cdf(-nsig1);
207 p[2] = 0.5;
208 p[3] = ROOT::Math::normal_cdf(nsig1);
209 p[4] = ROOT::Math::normal_cdf(nsig2);
210
211 bool resultIsAsymptotic(false);
212 if (nEntries>=1) {
213 HypoTestResult* r = dynamic_cast<HypoTestResult *> ( GetResult(0) );
214 assert(r!=0);
215 if ( !r->GetNullDistribution() && !r->GetAltDistribution() ) {
216 resultIsAsymptotic = true;
217 }
218 }
219
220 int nPointsRemoved(0);
221
222 double CLsobsprev(1.0);
223 std::vector<double>::iterator itr = fXValues.begin();
224
225 for (; itr!=fXValues.end();) {
226
227 double x = (*itr);
228 int i = FindIndex(x);
229 //HypoTestResult* oneresult = GetResult(i);
230
232 if (!s) break;
233
234 /////////////////////////////////////////////////////////////////////////////////////////
235
236 const std::vector<double> & values = s->GetSamplingDistribution();
237 if ((int) values.size() != fgAsymptoticNumPoints) {
238 oocoutE(this,Eval) << "HypoTestInverterResult::ExclusionCleanup - invalid size of sampling distribution" << std::endl;
239 delete s;
240 break;
241 }
242
243 /// expected p-values
244 // special case for asymptotic results (cannot use TMath::quantile in that case)
245 if (resultIsAsymptotic) {
246 double maxSigma = fgAsymptoticMaxSigma;
247 double dsig = 2.*maxSigma / (values.size() -1) ;
248 int i0 = (int) TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5 );
249 int i1 = (int) TMath::Floor ( ( -nsig1 + maxSigma )/dsig + 0.5 );
250 int i2 = (int) TMath::Floor ( ( maxSigma )/dsig + 0.5 );
251 int i3 = (int) TMath::Floor ( ( nsig1 + maxSigma )/dsig + 0.5 );
252 int i4 = (int) TMath::Floor ( ( nsig2 + maxSigma )/dsig + 0.5 );
253 //
254 q[0] = values[i0];
255 q[1] = values[i1];
256 q[2] = values[i2];
257 q[3] = values[i3];
258 q[4] = values[i4];
259 } else {
260 double * z = const_cast<double *>( &values[0] ); // need to change TMath::Quantiles
261 TMath::Quantiles(values.size(), 5, z, q, p, false);
262 }
263
264 delete s;
265
266 /// store useful quantities for reuse later ...
267 /// http://root.cern.ch/root/html532/src/RooStats__HypoTestInverterPlot.cxx.html#197
268 for (int j=0; j<5; ++j) { qv[j]=q[j]; }
269 qv[5] = CLs(i) ; //
270 qv[6] = CLsError(i) ; //
271 qv[7] = CLb(i) ; //
272 qv[8] = CLbError(i) ; //
273 qv[9] = CLsplusb(i) ; //
274 qv[10] = CLsplusbError(i) ; //
275 double CLsobs = qv[5];
276
277 /////////////////////////////////////////////////////////////////////////////////////////
278
279 bool removeThisPoint(false);
280
281 // 1. CLs should drop, else skip this point
282 if (!removeThisPoint && resultIsAsymptotic && i>=1 && CLsobs>CLsobsprev) {
283 //StatToolsLogger << kERROR << "Asymptotic. CLs not dropping: " << CLsobs << ". Remove this point." << GEndl;
284 removeThisPoint = true;
285 } else { CLsobsprev = CLsobs; }
286
287 // 2. CLs should not spike, else skip this point
288 if (!removeThisPoint && i>=1 && CLsobs>=0.9999) {
289 //StatToolsLogger << kERROR << "CLs spiking at 1.0: " << CLsobs << ". Remove this point." << GEndl;
290 removeThisPoint = true;
291 }
292 // 3. Not interested in CLs values that become too low.
293 if (!removeThisPoint && i>=1 && qv[4]<fCLsCleanupThreshold) { removeThisPoint = true; }
294
295 // to remove or not to remove
296 if (removeThisPoint) {
297 itr = fXValues.erase(itr); // returned itr has been updated.
300 nPointsRemoved++;
301 continue;
302 } else { // keep
303 CLsobsprev = CLsobs;
304 ++itr;
305 }
306 }
307
308 // after cleanup, reset existing limits
309 fFittedUpperLimit = false;
310 fFittedLowerLimit = false;
312
313 return nPointsRemoved;
314}
315
316////////////////////////////////////////////////////////////////////////////////
317/// Merge this HypoTestInverterResult with another
318/// HypoTestInverterResult passed as argument
319/// The merge is done by combining the HypoTestResult when the same point value exist in both results.
320/// If results exist at different points these are added in the new result
321/// NOTE: Merging of the expected p-values obtained with pseudo-data.
322/// When expected p-values exist in the result (i.e. when rebuild option is used when getting the expected
323/// limit distribution in the HYpoTestInverter) then the expected p-values are also merged. This is equivalent
324/// at merging the pseudo-data. However there can be an inconsistency if the expected p-values have been
325/// obtained with different toys. In this case the merge is done but a warning message is printed.
326
328{
329 int nThis = ArraySize();
330 int nOther = otherResult.ArraySize();
331 if (nOther == 0) return true;
332 if (nOther != otherResult.fYObjects.GetSize() ) return false;
333 if (nThis != fYObjects.GetSize() ) return false;
334
335 // cannot merge in case of inconsistent members
336 if (fExpPValues.GetSize() > 0 && fExpPValues.GetSize() != nThis ) return false;
337 if (otherResult.fExpPValues.GetSize() > 0 && otherResult.fExpPValues.GetSize() != nOther ) return false;
338
339 oocoutI(this,Eval) << "HypoTestInverterResult::Add - merging result from " << otherResult.GetName()
340 << " in " << GetName() << std::endl;
341
342 bool addExpPValues = (fExpPValues.GetSize() == 0 && otherResult.fExpPValues.GetSize() > 0);
343 bool mergeExpPValues = (fExpPValues.GetSize() > 0 && otherResult.fExpPValues.GetSize() > 0);
344
345 if (addExpPValues || mergeExpPValues)
346 oocoutI(this,Eval) << "HypoTestInverterResult::Add - merging also the expected p-values from pseudo-data" << std::endl;
347
348
349 // case current result is empty
350 // just make a simple copy of the other result
351 if (nThis == 0) {
352 fXValues = otherResult.fXValues;
353 for (int i = 0; i < nOther; ++i)
354 fYObjects.Add( otherResult.fYObjects.At(i)->Clone() );
355 for (int i = 0; i < fExpPValues.GetSize() ; ++i)
356 fExpPValues.Add( otherResult.fExpPValues.At(i)->Clone() );
357 }
358 // now do the real merge combining point with same value or adding extra ones
359 else {
360 for (int i = 0; i < nOther; ++i) {
361 double otherVal = otherResult.fXValues[i];
362 HypoTestResult * otherHTR = (HypoTestResult*) otherResult.fYObjects.At(i);
363 if (otherHTR == 0) continue;
364 bool sameXFound = false;
365 for (int j = 0; j < nThis; ++j) {
366 double thisVal = fXValues[j];
367
368 // if same value merge the result
369 if ( (std::abs(otherVal) < 1 && TMath::AreEqualAbs(otherVal, thisVal,1.E-12) ) ||
370 (std::abs(otherVal) >= 1 && TMath::AreEqualRel(otherVal, thisVal,1.E-12) ) ) {
371 HypoTestResult * thisHTR = (HypoTestResult*) fYObjects.At(j);
372 thisHTR->Append(otherHTR);
373 sameXFound = true;
374 if (mergeExpPValues) {
376 // check if same toys have been used for the test statistic distribution
377 int thisNToys = (thisHTR->GetNullDistribution() ) ? thisHTR->GetNullDistribution()->GetSize() : 0;
378 int otherNToys = (otherHTR->GetNullDistribution() ) ? otherHTR->GetNullDistribution()->GetSize() : 0;
379 if (thisNToys != otherNToys )
380 oocoutW(this,Eval) << "HypoTestInverterResult::Add expected p values have been generated with different toys " << thisNToys << " , " << otherNToys << std::endl;
381 }
382 break;
383 }
384 }
385 if (!sameXFound) {
386 // add the new result
387 fYObjects.Add(otherHTR->Clone() );
388 fXValues.push_back( otherVal);
389 }
390 // add in any case also when same x found
391 if (addExpPValues)
392 fExpPValues.Add( otherResult.fExpPValues.At(i)->Clone() );
393
394
395 }
396 }
397
398 if (ArraySize() > nThis)
399 oocoutI(this,Eval) << "HypoTestInverterResult::Add - new number of points is " << fXValues.size()
400 << std::endl;
401 else
402 oocoutI(this,Eval) << "HypoTestInverterResult::Add - new toys/point is "
403 << ((HypoTestResult*) fYObjects.At(0))->GetNullDistribution()->GetSize()
404 << std::endl;
405
406 // reset cached limit values
409
410 return true;
411}
412
413////////////////////////////////////////////////////////////////////////////////
414/// Add a single point result (an HypoTestResult)
415
417{
418 int i= FindIndex(x);
419 if (i<0) {
420 fXValues.push_back(x);
421 fYObjects.Add(res.Clone());
422 } else {
424 if (!r) return false;
425 r->Append(&res);
426 }
427
428 // reset cached limit values
431
432 return true;
433}
434
435////////////////////////////////////////////////////////////////////////////////
436/// function to return the value of the parameter of interest for the i^th entry in the results
437
438double HypoTestInverterResult::GetXValue( int index ) const
439{
440 if ( index >= ArraySize() || index<0 ) {
441 oocoutE(this,InputArguments) << "Problem: You are asking for an impossible array index value\n";
442 return -999;
443 }
444
445 return fXValues[index];
446}
447
448////////////////////////////////////////////////////////////////////////////////
449/// function to return the value of the confidence level for the i^th entry in the results
450
451double HypoTestInverterResult::GetYValue( int index ) const
452{
453 auto result = GetResult(index);
454 if ( !result ) {
455 return -999;
456 }
457
458 if (fUseCLs) {
459 return result->CLs();
460 } else {
461 return result->CLsplusb(); // CLs+b
462 }
463}
464
465////////////////////////////////////////////////////////////////////////////////
466/// function to return the estimated error on the value of the confidence level for the i^th entry in the results
467
468double HypoTestInverterResult::GetYError( int index ) const
469{
470 auto result = GetResult(index);
471 if ( !result ) {
472 return -999;
473 }
474
475 if (fUseCLs) {
476 return result->CLsError();
477 } else {
478 return result->CLsplusbError();
479 }
480}
481
482////////////////////////////////////////////////////////////////////////////////
483/// function to return the observed CLb value for the i-th entry
484
485double HypoTestInverterResult::CLb( int index ) const
486{
487 auto result = GetResult(index);
488 if ( !result ) {
489 return -999;
490 }
491 return result->CLb();
492}
493
494////////////////////////////////////////////////////////////////////////////////
495/// function to return the observed CLs+b value for the i-th entry
496
497double HypoTestInverterResult::CLsplusb( int index ) const
498{
499 auto result = GetResult(index);
500 if ( !result) {
501 return -999;
502 }
503 return result->CLsplusb();
504}
505
506////////////////////////////////////////////////////////////////////////////////
507/// function to return the observed CLs value for the i-th entry
508
509double HypoTestInverterResult::CLs( int index ) const
510{
511 auto result = GetResult(index);
512 if ( !result ) {
513 return -999;
514 }
515 return result->CLs();
516}
517
518////////////////////////////////////////////////////////////////////////////////
519/// function to return the error on the observed CLb value for the i-th entry
520
521double HypoTestInverterResult::CLbError( int index ) const
522{
523 auto result = GetResult(index);
524 if ( !result ) {
525 return -999;
526 }
527 return result->CLbError();
528}
529
530////////////////////////////////////////////////////////////////////////////////
531/// function to return the error on the observed CLs+b value for the i-th entry
532
534{
535 auto result = GetResult(index);
536 if ( ! result ) {
537 return -999;
538 }
539 return result->CLsplusbError();
540}
541
542////////////////////////////////////////////////////////////////////////////////
543/// function to return the error on the observed CLs value for the i-th entry
544
545double HypoTestInverterResult::CLsError( int index ) const
546{
547 auto result = GetResult(index);
548 if ( ! result ){
549 return -999;
550 }
551 return result->CLsError();
552}
553
554////////////////////////////////////////////////////////////////////////////////
555/// get the HypoTestResult object at the given index point
556
558{
559 if ( index >= ArraySize() || index<0 ) {
560 oocoutE(this,InputArguments) << "Problem: You are asking for an impossible array index value\n";
561 return 0;
562 }
563
564 return ((HypoTestResult*) fYObjects.At(index));
565}
566
567////////////////////////////////////////////////////////////////////////////////
568/// find the index corresponding at the poi value xvalue
569/// If no points is found return -1
570/// Note that a tolerance is used of 10^-12 to find the closest point
571
572int HypoTestInverterResult::FindIndex(double xvalue) const
573{
574 const double tol = 1.E-12;
575 for (int i=0; i<ArraySize(); i++) {
576 double xpoint = fXValues[i];
577 if ( (std::abs(xvalue) > 1 && TMath::AreEqualRel( xvalue, xpoint, tol) ) ||
578 (std::abs(xvalue) < 1 && TMath::AreEqualAbs( xvalue, xpoint, tol) ) )
579 return i;
580 }
581 return -1;
582}
583
584////////////////////////////////////////////////////////////////////////////////
585/// return the X value of the given graph for the target value y0
586/// the graph is evaluated using linear interpolation by default.
587/// if option = "S" a TSpline3 is used
588
589double HypoTestInverterResult::GetGraphX(const TGraph & graph, double y0, bool lowSearch, double &axmin, double &axmax) const {
590
591//#define DO_DEBUG
592#ifdef DO_DEBUG
593 std::cout << "using graph for search " << lowSearch << " min " << axmin << " max " << axmax << std::endl;
594#endif
595
596
597 // find reasonable xmin and xmax if not given
598 const double * y = graph.GetY();
599 int n = graph.GetN();
600 if (n < 2) {
601 ooccoutE(this,Eval) << "HypoTestInverterResult::GetGraphX - need at least 2 points for interpolation (n=" << n << ")\n";
602 return (n>0) ? y[0] : 0;
603 }
604
605 double varmin = - TMath::Infinity();
606 double varmax = TMath::Infinity();
607 const RooRealVar* var = dynamic_cast<RooRealVar*>( fParameters.first() );
608 if (var) {
609 varmin = var->getMin();
610 varmax = var->getMax();
611 }
612
613
614 // find ymin and ymax and corresponding values
615 double ymin = TMath::MinElement(n,y);
616 double ymax = TMath::MaxElement(n,y);
617 // cannot find intercept in the full range - return min or max valie
618 if (ymax < y0) {
619 return (lowSearch) ? varmax : varmin;
620 }
621 if (ymin > y0) {
622 return (lowSearch) ? varmin : varmax;
623 }
624
625 double xmin = axmin;
626 double xmax = axmax;
627
628 // case no range is given check if need to extrapolate to lower/upper values
629 if (axmin >= axmax ) {
630
631#ifdef DO_DEBUG
632 std::cout << "No rage given - check if extrapolation is needed " << std::endl;
633#endif
634
635 xmin = graph.GetX()[0];
636 xmax = graph.GetX()[n-1];
637
638 double yfirst = graph.GetY()[0];
639 double ylast = graph.GetY()[n-1];
640
641 // distinguish the case we have lower /upper limits
642 // check if a possible crossing exists otherwise return variable min/max
643
644 // do lower extrapolation
645 if ( (ymax < y0 && !lowSearch) || ( yfirst > y0 && lowSearch) ) {
646 xmin = varmin;
647 }
648 // do upper extrapolation
649 if ( (ymax < y0 && lowSearch) || ( ylast > y0 && !lowSearch) ) {
650 xmax = varmax;
651 }
652 }
653
654 auto func = [&](double x) {
655 return (fInterpolOption == kSpline) ? graph.Eval(x, nullptr, "S") - y0 : graph.Eval(x) - y0;
656 };
657 ROOT::Math::Functor1D f1d(func);
658
660 brf.SetFunction(f1d,xmin,xmax);
661 brf.SetNpx(TMath::Max(graph.GetN()*2,100) );
662#ifdef DO_DEBUG
663 std::cout << "findind root for " << xmin << " , "<< xmax << "f(x) : " << graph.Eval(xmin) << " , " << graph.Eval(0.5*(xmax+xmin))
664 << " , " << graph.Eval(xmax) << " target " << y0 << std::endl;
665#endif
666
667 bool ret = brf.Solve(100, 1.E-16, 1.E-6);
668 if (!ret) {
669 ooccoutE(this,Eval) << "HypoTestInverterResult - interpolation failed for interval [" << xmin << "," << xmax
670 << " ] g(xmin,xmax) =" << graph.Eval(xmin) << "," << graph.Eval(xmax)
671 << " target=" << y0 << " return inf" << std::endl;
672 return TMath::Infinity();
673 }
674 double limit = brf.Root();
675
676 // auto grfunc = [&](double * x, double *) {
677 // return (fInterpolOption == kSpline) ? graph.Eval(*x, nullptr, "S") - y0 : graph.Eval(*x) - y0;
678 // };
679 // TF1 tgrfunc("tgrfunc",grfunc,xmin,xmax,0);
680 // double limit = tgrfunc.GetX(0,xmin,xmax);
681
682#ifdef DO_DEBUG
683 if (lowSearch) std::cout << "lower limit search : ";
684 else std::cout << "Upper limit search : ";
685 std::cout << "interpolation done between " << xmin << " and " << xmax
686 << "\n Found limit using RootFinder is " << limit << std::endl;
687
688 TString fname = "graph_upper.root";
689 if (lowSearch) fname = "graph_lower.root";
690 auto file = TFile::Open(fname,"RECREATE");
691 graph.Write("graph");
692 file->Close();
693#endif
694
695 // look in case if a new intersection exists
696 // only when boundaries are not given
697 if (axmin >= axmax) {
698 int index = TMath::BinarySearch(n, graph.GetX(), limit);
699#ifdef DO_DEBUG
700 std::cout << "do new interpolation dividing from " << index << " and " << y[index] << std::endl;
701#endif
702
703 if (lowSearch && index >= 1 && (y[0] - y0) * ( y[index]- y0) < 0) {
704 //search if another intersection exists at a lower value
705 limit = GetGraphX(graph, y0, lowSearch, graph.GetX()[0], graph.GetX()[index] );
706 }
707 else if (!lowSearch && index < n-2 && (y[n-1] - y0) * ( y[index+1]- y0) < 0) {
708 // another intersection exists at an higher value
709 limit = GetGraphX(graph, y0, lowSearch, graph.GetX()[index+1], graph.GetX()[n-1] );
710 }
711 }
712 // return also xmin, xmax values
713 axmin = xmin;
714 axmax = xmax;
715
716 return limit;
717}
718
719////////////////////////////////////////////////////////////////////////////////
720/// interpolate to find a limit value
721/// Use a linear or a spline interpolation depending on the interpolation option
722
723double HypoTestInverterResult::FindInterpolatedLimit(double target, bool lowSearch, double xmin, double xmax )
724{
725
726 // variable minimum and maximum
727 double varmin = - TMath::Infinity();
728 double varmax = TMath::Infinity();
729 const RooRealVar* var = dynamic_cast<RooRealVar*>( fParameters.first() );
730 if (var) {
731 varmin = var->getMin();
732 varmax = var->getMax();
733 }
734
735 if (ArraySize()<2) {
736 double val = (lowSearch) ? xmin : xmax;
737 oocoutW(this,Eval) << "HypoTestInverterResult::FindInterpolatedLimit"
738 << " - not enough points to get the inverted interval - return "
739 << val << std::endl;
740 fLowerLimit = varmin;
741 fUpperLimit = varmax;
742 return (lowSearch) ? fLowerLimit : fUpperLimit;
743 }
744
745 // sort the values in x
746 int n = ArraySize();
747 std::vector<unsigned int> index(n );
748 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
749 // make a graph with the sorted point
750 TGraph graph(n);
751 for (int i = 0; i < n; ++i)
752 graph.SetPoint(i, GetXValue(index[i]), GetYValue(index[i] ) );
753
754
755 //std::cout << " search for " << lowSearch << " xmin = " << xmin << " xmax " << xmax << std::endl;
756
757
758 // search first for min/max in the given range
759 if (xmin >= xmax) {
760
761
762 // search for maximum between the point
763 double * itrmax = std::max_element(graph.GetY() , graph.GetY() +n);
764 double ymax = *itrmax;
765 int iymax = itrmax - graph.GetY();
766 double xwithymax = graph.GetX()[iymax];
767
768#ifdef DO_DEBUG
769 std::cout << " max of y " << iymax << " " << xwithymax << " " << ymax << " target is " << target << std::endl;
770#endif
771 // look if maximum is above/below target
772 if (ymax > target) {
773 if (lowSearch) {
774 if ( iymax > 0) {
775 // low search (minimum is first point or minimum range)
776 xmin = ( graph.GetY()[0] <= target ) ? graph.GetX()[0] : varmin;
777 xmax = xwithymax;
778 }
779 else {
780 // no room for lower limit
781 fLowerLimit = varmin;
782 return fLowerLimit;
783 }
784 }
785 if (!lowSearch ) {
786 // up search
787 if ( iymax < n-1 ) {
788 xmin = xwithymax;
789 xmax = ( graph.GetY()[n-1] <= target ) ? graph.GetX()[n-1] : varmax;
790 }
791 else {
792 // no room for upper limit
793 fUpperLimit = varmax;
794 return fUpperLimit;
795 }
796 }
797 }
798 else {
799 // in case is below the target
800 // find out if is a lower or upper search
801 if (iymax <= (n-1)/2 ) {
802 lowSearch = false;
803 fLowerLimit = varmin;
804 }
805 else {
806 lowSearch = true;
807 fUpperLimit = varmax;
808 }
809 }
810
811#ifdef DO_DEBUG
812 std::cout << " found xmin, xmax = " << xmin << " " << xmax << " for search " << lowSearch << std::endl;
813#endif
814
815 // now come here if I have already found a lower/upper limit
816 // i.e. I am calling routine for the second time
817#ifdef ISNEEDED
818 // should not really come here
819 if (lowSearch && fUpperLimit < varmax) {
820 xmin = fXValues[ index.front() ];
821 // find xmax (is first point before upper limit)
822 int upI = FindClosestPointIndex(target, 2, fUpperLimit);
823 if (upI < 1) return xmin;
824 xmax = GetXValue(upI);
825 }
826 else if (!lowSearch && fLowerLimit > varmin ) {
827 // find xmin (is first point after lower limit)
828 int lowI = FindClosestPointIndex(target, 3, fLowerLimit);
829 if (lowI >= n-1) return xmax;
830 xmin = GetXValue(lowI);
831 xmax = fXValues[ index.back() ];
832 }
833#endif
834 }
835
836#ifdef DO_DEBUG
837 std::cout << "finding " << lowSearch << " limit between " << xmin << " " << xmax << endl;
838#endif
839
840 // compute noe the limit using the TGraph interpolations routine
841 double limit = GetGraphX(graph, target, lowSearch, xmin, xmax);
842 if (lowSearch) fLowerLimit = limit;
843 else fUpperLimit = limit;
844 // estimate the error
845 double error = CalculateEstimatedError( target, lowSearch, xmin, xmax);
846
847 TString limitType = (lowSearch) ? "lower" : "upper";
848 ooccoutD(this,Eval) << "HypoTestInverterResult::FindInterpolateLimit "
849 << "the computed " << limitType << " limit is " << limit << " +/- " << error << std::endl;
850
851#ifdef DO_DEBUG
852 std::cout << "Found limit is " << limit << " +/- " << error << std::endl;
853#endif
854
855
856 if (lowSearch) return fLowerLimit;
857 return fUpperLimit;
858
859
860// if (lowSearch && !TMath::IsNaN(fUpperLimit)) return fLowerLimit;
861// if (!lowSearch && !TMath::IsNaN(fLowerLimit)) return fUpperLimit;
862// // is this needed ?
863// // we call again the function for the upper limits
864
865// // now perform the opposite search on the complement interval
866// if (lowSearch) {
867// xmin = xmax;
868// xmax = varmax;
869// } else {
870// xmax = xmin;
871// xmin = varmin;
872// }
873// double limit2 = GetGraphX(graph, target, !lowSearch, xmin, xmax);
874// if (!lowSearch) fLowerLimit = limit2;
875// else fUpperLimit = limit2;
876
877// CalculateEstimatedError( target, !lowSearch, xmin, xmax);
878
879// #ifdef DO_DEBUG
880// std::cout << "other limit is " << limit2 << std::endl;
881// #endif
882
883// return (lowSearch) ? fLowerLimit : fUpperLimit;
884
885}
886
887////////////////////////////////////////////////////////////////////////////////
888/// - if mode = 0
889/// find closest point to target in Y, the object closest to the target which is 3 sigma from the target
890/// and has smaller error
891/// - if mode = 1
892/// find 2 closest point to target in X and between these two take the one closer to the target
893/// - if mode = 2 as in mode = 1 but return the lower point not the closest one
894/// - if mode = 3 as in mode = 1 but return the upper point not the closest one
895
896int HypoTestInverterResult::FindClosestPointIndex(double target, int mode, double xtarget)
897{
898
899 int bestIndex = -1;
900 int closestIndex = -1;
901 if (mode == 0) {
902 double smallestError = 2; // error must be < 1
903 double bestValue = 2;
904 for (int i=0; i<ArraySize(); i++) {
905 double dist = fabs(GetYValue(i)-target);
906 if ( dist <3 *GetYError(i) ) { // less than 1 sigma from target CL
907 if (GetYError(i) < smallestError ) {
908 smallestError = GetYError(i);
909 bestIndex = i;
910 }
911 }
912 if ( dist < bestValue) {
913 bestValue = dist;
914 closestIndex = i;
915 }
916 }
917 if (bestIndex >=0) return bestIndex;
918 // if no points found just return the closest one to the target
919 return closestIndex;
920 }
921 // else mode = 1,2,3
922 // find the two closest points to limit value
923 // sort the array first
924 int n = fXValues.size();
925 std::vector<unsigned int> indx(n);
926 TMath::SortItr(fXValues.begin(), fXValues.end(), indx.begin(), false);
927 std::vector<double> xsorted( n);
928 for (int i = 0; i < n; ++i) xsorted[i] = fXValues[indx[i] ];
929 int index1 = TMath::BinarySearch( n, &xsorted[0], xtarget);
930
931#ifdef DO_DEBUG
932 std::cout << "finding closest point to " << xtarget << " is " << index1 << " " << indx[index1] << std::endl;
933#endif
934
935 // case xtarget is outside the range (before or afterwards)
936 if (index1 < 0) return indx[0];
937 if (index1 >= n-1) return indx[n-1];
938 int index2 = index1 +1;
939
940 if (mode == 2) return (GetXValue(indx[index1]) < GetXValue(indx[index2])) ? indx[index1] : indx[index2];
941 if (mode == 3) return (GetXValue(indx[index1]) > GetXValue(indx[index2])) ? indx[index1] : indx[index2];
942 // get smaller point of the two (mode == 1)
943 if (fabs(GetYValue(indx[index1])-target) <= fabs(GetYValue(indx[index2])-target) )
944 return indx[index1];
945 return indx[index2];
946
947}
948
949////////////////////////////////////////////////////////////////////////////////
950
952{
953 if (fFittedLowerLimit) return fLowerLimit;
954 //std::cout << "finding point with cl = " << 1-(1-ConfidenceLevel())/2 << endl;
956 // find both lower/upper limit
958 } else {
959 //LM: I think this is never called
961 }
962 return fLowerLimit;
963}
964
965////////////////////////////////////////////////////////////////////////////////
966
968{
969 //std::cout << "finding point with cl = " << (1-ConfidenceLevel())/2 << endl;
970 if (fFittedUpperLimit) return fUpperLimit;
973 } else {
974 //LM: I think this is never called
976 }
977 return fUpperLimit;
978}
979
980////////////////////////////////////////////////////////////////////////////////
981/// Return an error estimate on the upper(lower) limit. This is the error on
982/// either CLs or CLsplusb divided by an estimate of the slope at this
983/// point.
984
985Double_t HypoTestInverterResult::CalculateEstimatedError(double target, bool lower, double xmin, double xmax)
986{
987
988 if (ArraySize()==0) {
989 oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimateError"
990 << "Empty result \n";
991 return 0;
992 }
993
994 if (ArraySize()<2) {
995 oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimateError"
996 << " only points - return its error\n";
997 return GetYError(0);
998 }
999
1000 // it does not make sense in case of asymptotic which do not have point errors
1001 if (!GetNullTestStatDist(0) ) return 0;
1002
1003 TString type = (!lower) ? "upper" : "lower";
1004
1005#ifdef DO_DEBUG
1006 std::cout << "calculate estimate error " << type << " between " << xmin << " and " << xmax << std::endl;
1007 std::cout << "computed limit is " << ( (lower) ? fLowerLimit : fUpperLimit ) << std::endl;
1008#endif
1009
1010 // make a TGraph Errors with the sorted points
1011 std::vector<unsigned int> indx(fXValues.size());
1012 TMath::SortItr(fXValues.begin(), fXValues.end(), indx.begin(), false);
1013 // make a graph with the sorted point
1015 int ip = 0, np = 0;
1016 for (int i = 0; i < ArraySize(); ++i) {
1017 if ( (xmin < xmax) && ( GetXValue(indx[i]) >= xmin && GetXValue(indx[i]) <= xmax) ) {
1018 np++;
1019 // exclude points with zero or very small errors
1020 if (GetYError(indx[i] ) > 1.E-6) {
1021 graph.SetPoint(ip, GetXValue(indx[i]), GetYValue(indx[i] ) );
1022 graph.SetPointError(ip, 0., GetYError(indx[i]) );
1023 ip++;
1024 }
1025 }
1026 }
1027 if (graph.GetN() < 2) {
1028 if (np >= 2) oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimatedError - no valid points - cannot estimate the " << type << " limit error " << std::endl;
1029 return 0;
1030 }
1031
1032 double minX = xmin;
1033 double maxX = xmax;
1034 if (xmin >= xmax) {
1035 minX = fXValues[ indx.front() ];
1036 maxX = fXValues[ indx.back() ];
1037 }
1038
1039 TF1 fct("fct", "exp([0] * (x - [2] ) + [1] * (x-[2])**2)", minX, maxX);
1040 double scale = maxX-minX;
1041 if (lower) {
1042 fct.SetParameters( 2./scale, 0.1/scale, graph.GetX()[0] );
1043 fct.SetParLimits(0,0,100./scale);
1044 fct.SetParLimits(1,0, 10./scale); }
1045 else {
1046 fct.SetParameters( -2./scale, -0.1/scale );
1047 fct.SetParLimits(0,-100./scale, 0);
1048 fct.SetParLimits(1,-100./scale, 0); }
1049
1050 if (graph.GetN() < 3) fct.FixParameter(1,0.);
1051
1052 // find the point closest to the limit
1053 double limit = (!lower) ? fUpperLimit : fLowerLimit;
1054 if (TMath::IsNaN(limit)) return 0; // cannot do if limit not computed
1055
1056
1057#ifdef DO_DEBUG
1058 TCanvas * c1 = new TCanvas();
1059 std::cout << "fitting for limit " << type << "between " << minX << " , " << maxX << " points considered " << graph.GetN() << std::endl;
1060 int fitstat = graph.Fit(&fct," EX0");
1061 graph.SetMarkerStyle(20);
1062 graph.Draw("AP");
1063 graph.Print();
1064 c1->SaveAs(TString::Format("graphFit_%s.pdf",type.Data()) );
1065 delete c1;
1066#else
1067 int fitstat = graph.Fit(&fct,"Q EX0");
1068#endif
1069
1070 int index = FindClosestPointIndex(target, 1, limit);
1071 double theError = 0;
1072 if (fitstat == 0) {
1073 double errY = GetYError(index);
1074 if (errY > 0) {
1075 double m = fct.Derivative( GetXValue(index) );
1076 theError = std::min(fabs( GetYError(index) / m), maxX-minX);
1077 }
1078 }
1079 else {
1080 oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimatedError - cannot estimate the " << type << " limit error " << std::endl;
1081 theError = 0;
1082 }
1083 if (lower)
1084 fLowerLimitError = theError;
1085 else
1086 fUpperLimitError = theError;
1087
1088#ifdef DO_DEBUG
1089 std::cout << "closes point to the limit is " << index << " " << GetXValue(index) << " and has error " << GetYError(index) << std::endl;
1090#endif
1091
1092 return theError;
1093}
1094
1095////////////////////////////////////////////////////////////////////////////////
1096/// need to have compute first lower limit
1097
1099{
1101 if (fLowerLimitError >= 0) return fLowerLimitError;
1102 // try to recompute the error
1103 return CalculateEstimatedError(1-ConfidenceLevel(), true);
1104}
1105
1106////////////////////////////////////////////////////////////////////////////////
1107
1109{
1111 if (fUpperLimitError >= 0) return fUpperLimitError;
1112 // try to recompute the error
1113 return CalculateEstimatedError(1-ConfidenceLevel(), false);
1114}
1115
1116////////////////////////////////////////////////////////////////////////////////
1117/// get the background test statistic distribution
1118
1120
1121 HypoTestResult * firstResult = (HypoTestResult*) fYObjects.At(index);
1122 if (!firstResult) return 0;
1123 return firstResult->GetBackGroundIsAlt() ? firstResult->GetAltDistribution() : firstResult->GetNullDistribution();
1124}
1125
1126////////////////////////////////////////////////////////////////////////////////
1127/// get the signal and background test statistic distribution
1128
1130 HypoTestResult * result = (HypoTestResult*) fYObjects.At(index);
1131 if (!result) return 0;
1132 return !result->GetBackGroundIsAlt() ? result->GetAltDistribution() : result->GetNullDistribution();
1133}
1134
1135////////////////////////////////////////////////////////////////////////////////
1136/// get the expected p-value distribution at the scanned point index
1137
1139
1140 if (index < 0 || index >= ArraySize() ) return 0;
1141
1142 if (fExpPValues.GetSize() == ArraySize() ) {
1143 return (SamplingDistribution*) fExpPValues.At(index)->Clone();
1144 }
1145
1146 static bool useFirstB = false;
1147 // get the S+B distribution
1148 int bIndex = (useFirstB) ? 0 : index;
1149
1150 SamplingDistribution * bDistribution = GetBackgroundTestStatDist(bIndex);
1152
1153 HypoTestResult * result = (HypoTestResult*) fYObjects.At(index);
1154
1155 if (bDistribution && sbDistribution) {
1156
1157 // create a new HypoTestResult
1158 HypoTestResult tempResult;
1159 tempResult.SetPValueIsRightTail( result->GetPValueIsRightTail() );
1160 tempResult.SetBackgroundAsAlt( true);
1161 // ownership of SamplingDistribution is in HypoTestResult class now
1162 tempResult.SetNullDistribution( new SamplingDistribution(*sbDistribution) );
1163 tempResult.SetAltDistribution( new SamplingDistribution(*bDistribution ) );
1164
1165 std::vector<double> values(bDistribution->GetSize());
1166 for (int i = 0; i < bDistribution->GetSize(); ++i) {
1167 tempResult.SetTestStatisticData( bDistribution->GetSamplingDistribution()[i] );
1168 values[i] = (fUseCLs) ? tempResult.CLs() : tempResult.CLsplusb();
1169 }
1170 return new SamplingDistribution("expected values","expected values",values);
1171 }
1172 // in case b abs sbDistribution are null assume is coming from the asymptotic calculator
1173 // hard -coded this value (no really needed to be used by user)
1176 const double smax = fgAsymptoticMaxSigma;
1177 const int npoints = fgAsymptoticNumPoints;
1178 const double dsig = 2* smax/ (npoints-1) ;
1179 std::vector<double> values(npoints);
1180 for (int i = 0; i < npoints; ++i) {
1181 double nsig = -smax + dsig*i;
1182 double pval = AsymptoticCalculator::GetExpectedPValues( result->NullPValue(), result->AlternatePValue(), nsig, fUseCLs, !fIsTwoSided);
1183 if (pval < 0) { return 0;}
1184 values[i] = pval;
1185 }
1186 return new SamplingDistribution("Asymptotic expected values","Asymptotic expected values",values);
1187
1188}
1189
1190////////////////////////////////////////////////////////////////////////////////
1191/// get the limit distribution (lower/upper depending on the flag)
1192/// by interpolating the expected p values for each point
1193
1195 if (ArraySize()<2) {
1196 oocoutE(this,Eval) << "HypoTestInverterResult::GetLimitDistribution"
1197 << " not enough points - return 0 " << std::endl;
1198 return 0;
1199 }
1200
1201 ooccoutD(this,Eval) << "HypoTestInverterResult - computing limit distribution...." << std::endl;
1202
1203
1204 // find optimal size by looking at the PValue distribution obtained
1205 int npoints = ArraySize();
1206 std::vector<SamplingDistribution*> distVec( npoints );
1207 double sum = 0;
1208 for (unsigned int i = 0; i < distVec.size(); ++i) {
1209 distVec[i] = GetExpectedPValueDist(i);
1210 // sort the distributions
1211 // hack (by calling InverseCDF(0) will sort the sampling distribution
1212 if (distVec[i] ) {
1213 distVec[i]->InverseCDF(0);
1214 sum += distVec[i]->GetSize();
1215 }
1216 }
1217 int size = int( sum/ npoints);
1218
1219 if (size < 10) {
1220 ooccoutW(this,InputArguments) << "HypoTestInverterResult - set a minimum size of 10 for limit distribution" << std::endl;
1221 size = 10;
1222 }
1223
1224
1225 double target = 1-fConfidenceLevel;
1226
1227 // vector with the quantiles of the p-values for each scanned poi point
1228 std::vector< std::vector<double> > quantVec(npoints );
1229 for (int i = 0; i < npoints; ++i) {
1230
1231 if (!distVec[i]) continue;
1232
1233 // make quantiles from the sampling distributions of the expected p values
1234 std::vector<double> pvalues = distVec[i]->GetSamplingDistribution();
1235 delete distVec[i]; distVec[i] = 0;
1236 std::sort(pvalues.begin(), pvalues.end());
1237 // find the quantiles of the distribution
1238 double p[1] = {0};
1239 double q[1] = {0};
1240
1241 quantVec[i] = std::vector<double>(size);
1242 for (int ibin = 0; ibin < size; ++ibin) {
1243 // exclude for a bug in TMath::Quantiles last item
1244 p[0] = std::min( (ibin+1) * 1./double(size), 1.0);
1245 // use the type 1 which give the point value
1246 TMath::Quantiles(pvalues.size(), 1, &pvalues[0], q, p, true, 0, 1 );
1247 (quantVec[i])[ibin] = q[0];
1248 }
1249
1250 }
1251
1252 // sort the values in x
1253 std::vector<unsigned int> index( npoints );
1254 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1255
1256 // SamplingDistribution * dist0 = distVec[index[0]];
1257 // SamplingDistribution * dist1 = distVec[index[1]];
1258
1259
1260 std::vector<double> limits(size);
1261 // loop on the p values and find the limit for each expected point in the quantiles vector
1262 for (int j = 0; j < size; ++j ) {
1263
1264 TGraph g;
1265 for (int k = 0; k < npoints ; ++k) {
1266 if (quantVec[index[k]].size() > 0 )
1267 g.SetPoint(g.GetN(), GetXValue(index[k]), (quantVec[index[k]])[j] );
1268 }
1269
1270 limits[j] = GetGraphX(g, target, lower);
1271
1272 }
1273
1274
1275 if (lower)
1276 return new SamplingDistribution("Expected lower Limit","Expected lower limits",limits);
1277 else
1278 return new SamplingDistribution("Expected upper Limit","Expected upper limits",limits);
1279
1280}
1281
1282////////////////////////////////////////////////////////////////////////////////
1283/// Get the expected lower limit
1284/// nsig is used to specify which expected value of the UpperLimitDistribution
1285/// For example
1286/// - nsig = 0 (default value) returns the expected value
1287/// - nsig = -1 returns the lower band value at -1 sigma
1288/// - nsig + 1 return the upper value
1289/// - opt = "" (default) : compute limit by interpolating all the p values, find the corresponding limit distribution
1290/// and then find the quantiles in the limit distribution
1291/// ioption = "P" is the method used for plotting. One Finds the corresponding nsig quantile in the p values and then
1292/// interpolates them
1293
1294double HypoTestInverterResult::GetExpectedLowerLimit(double nsig, const char * opt ) const {
1295
1296 return GetExpectedLimit(nsig, true, opt);
1297}
1298
1299////////////////////////////////////////////////////////////////////////////////
1300/// Get the expected upper limit
1301/// nsig is used to specify which expected value of the UpperLimitDistribution
1302/// For example
1303/// - nsig = 0 (default value) returns the expected value
1304/// - nsig = -1 returns the lower band value at -1 sigma
1305/// - nsig + 1 return the upper value
1306/// - opt is an option specifying the type of method used for computing the upper limit
1307/// - opt = "" (default) : compute limit by interpolating all the p values, find the corresponding limit distribution
1308/// and then find the quantiles in the limit distribution
1309/// ioption = "P" is the method used for plotting. One Finds the corresponding nsig quantile in the p values and then
1310/// interpolates them
1311
1312double HypoTestInverterResult::GetExpectedUpperLimit(double nsig, const char * opt ) const {
1313
1314 return GetExpectedLimit(nsig, false, opt);
1315}
1316
1317////////////////////////////////////////////////////////////////////////////////
1318/// get expected limit (lower/upper) depending on the flag
1319/// for asymptotic is a special case (the distribution is generated an step in sigma values)
1320/// distinguish asymptotic looking at the hypotest results
1321/// if option = "P" get expected limit using directly quantiles of p value distribution
1322/// else (default) find expected limit by obtaining first a full limit distributions
1323/// The last one is in general more correct
1324
1325double HypoTestInverterResult::GetExpectedLimit(double nsig, bool lower, const char * opt ) const {
1326
1327 const int nEntries = ArraySize();
1328 if (nEntries <= 0) return (lower) ? 1 : 0; // return 1 for lower, 0 for upper
1329
1330 HypoTestResult * r = dynamic_cast<HypoTestResult *> (fYObjects.First() );
1331 assert(r != 0);
1332 if (!r->GetNullDistribution() && !r->GetAltDistribution() ) {
1333 // we are in the asymptotic case
1334 // get the limits obtained at the different sigma values
1335 SamplingDistribution * limitDist = GetLimitDistribution(lower);
1336 if (!limitDist) return 0;
1337 const std::vector<double> & values = limitDist->GetSamplingDistribution();
1338 if (values.size() <= 1) return 0;
1339 double dsig = 2* fgAsymptoticMaxSigma/ (values.size() -1) ;
1340 int i = (int) TMath::Floor ( (nsig + fgAsymptoticMaxSigma)/dsig + 0.5);
1341 return values[i];
1342 }
1343
1344 double p[1] = {0};
1345 double q[1] = {0};
1346 p[0] = ROOT::Math::normal_cdf(nsig,1);
1347
1348 // for CLs+b can get the quantiles of p-value distribution and
1349 // interpolate them
1350 // In the case of CLs (since it is not a real p-value anymore but a ratio)
1351 // then it is needed to get first all limit distribution values and then the quantiles
1352
1353 TString option(opt);
1354 option.ToUpper();
1355 if (option.Contains("P")) {
1356
1357 TGraph g;
1358
1359 // sort the arrays based on the x values
1360 std::vector<unsigned int> index(nEntries);
1361 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1362
1363 for (int j=0; j<nEntries; ++j) {
1364 int i = index[j]; // i is the order index
1366 if (!s) {
1367 ooccoutI(this,Eval) << "HypoTestInverterResult - cannot compute expected p value distribution for point, x = "
1368 << GetXValue(i) << " skip it " << std::endl;
1369 continue;
1370 }
1371 const std::vector<double> & values = s->GetSamplingDistribution();
1372 double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
1373 TMath::Quantiles(values.size(), 1, x,q,p,false);
1374 g.SetPoint(g.GetN(), fXValues[i], q[0] );
1375 delete s;
1376 }
1377 if (g.GetN() < 2) {
1378 ooccoutE(this,Eval) << "HypoTestInverterResult - cannot compute limits , not enough points, n = " << g.GetN() << std::endl;
1379 return 0;
1380 }
1381
1382 // interpolate the graph to obtain the limit
1383 double target = 1-fConfidenceLevel;
1384 return GetGraphX(g, target, lower);
1385
1386 }
1387 // here need to use the limit distribution
1388 SamplingDistribution * limitDist = GetLimitDistribution(lower);
1389 if (!limitDist) return 0;
1390 const std::vector<double> & values = limitDist->GetSamplingDistribution();
1391 double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
1392 TMath::Quantiles(values.size(), 1, x,q,p,false);
1393 return q[0];
1394
1395}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define g(i)
Definition: RSha256.hxx:105
#define oocoutW(o, a)
Definition: RooMsgService.h:46
#define oocoutE(o, a)
Definition: RooMsgService.h:47
#define oocoutI(o, a)
Definition: RooMsgService.h:44
#define ooccoutW(o, a)
Definition: RooMsgService.h:53
#define ooccoutI(o, a)
Definition: RooMsgService.h:51
#define ooccoutD(o, a)
Definition: RooMsgService.h:50
#define ooccoutE(o, a)
Definition: RooMsgService.h:54
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
float xmin
Definition: THbookFile.cxx:93
float * q
Definition: THbookFile.cxx:87
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float ymax
Definition: THbookFile.cxx:93
Class for finding the root of a one dimensional function using the Brent algorithm.
double Root() const
Returns root value.
bool Solve(int maxIter=100, double absTol=1E-8, double relTol=1E-10)
Returns the X value corresponding to the function value fy for (xmin<x<xmax).
void SetNpx(int npx)
Set the number of point used to bracket root using a grid.
bool SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup)
Sets the function for the rest of the algorithms.
Functor1D class for one-dimensional functions.
Definition: Functor.h:487
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
RooAbsArg * first() const
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.
virtual Bool_t addOwned(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling addOwned() for each element in the source...
Definition: RooArgSet.h:92
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual TObject * clone(const char *newname) const
Definition: RooRealVar.h:48
static double GetExpectedPValues(double pnull, double palt, double nsigma, bool usecls, bool oneSided=true)
function given the null and the alt p value - return the expected one given the N - sigma value
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
Double_t LowerLimit()
lower and upper bound of the confidence interval (to get upper/lower limits, multiply the size( = 1-c...
Double_t LowerLimitEstimatedError()
rough estimation of the error on the computed bound of the confidence interval Estimate of lower limi...
double GetYValue(int index) const
function to return the value of the confidence level for the i^th entry in the results
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
get the signal and background test statistic distribution
Double_t UpperLimitEstimatedError()
Estimate of lower limit error function evaluates only a rough error on the lower limit.
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
double fLowerLimitError
interpolation option (linear or spline)
HypoTestInverterResult & operator=(const HypoTestInverterResult &other)
operator =
double CalculateEstimatedError(double target, bool lower=true, double xmin=1, double xmax=0)
Return an error estimate on the upper(lower) limit.
HypoTestInverterResult(const char *name=0)
default constructor
int ArraySize() const
number of entries in the results array
double FindInterpolatedLimit(double target, bool lowSearch=false, double xmin=1, double xmax=0)
interpolate to find a limit value Use a linear or a spline interpolation depending on the interpolati...
bool fInterpolateLowerLimit
two sided scan (look for lower/upper limit)
std::vector< double > fXValues
number of points used to build expected p-values
double CLsplusbError(int index) const
return the observed CLsplusb value for the i-th entry
double CLsError(int index) const
return the observed CLb value for the i-th entry
double GetGraphX(const TGraph &g, double y0, bool lowSearch, double &xmin, double &xmax) const
return the X value of the given graph for the target value y0 the graph is evaluated using linear int...
double CLbError(int index) const
return the observed CLb value for the i-th entry
double GetExpectedUpperLimit(double nsig=0, const char *opt="") const
get Limit value corresponding at the desired nsigma level (0) is median -1 sigma is 1 sigma
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results
bool Add(const HypoTestInverterResult &otherResult)
merge with the content of another HypoTestInverterResult object
int FindClosestPointIndex(double target, int mode=0, double xtarget=0)
TList fExpPValues
list of HypoTestResult for each point
double GetExpectedLowerLimit(double nsig=0, const char *opt="") const
get Limit value corresponding at the desired nsigma level (0) is median -1 sigma is 1 sigma
static int fgAsymptoticNumPoints
max sigma value used to scan asymptotic expected p values
int ExclusionCleanup()
remove points that appear to have failed.
double CLs(int index) const
return the observed CLb value for the i-th entry
int FindIndex(double xvalue) const
find the index corresponding at the poi value xvalue If no points is found return -1 Note that a tole...
double GetYError(int index) const
function to return the estimated error on the value of the confidence level for the i^th entry in the...
double CLb(int index) const
return the observed CLb value for the i-th entry
double GetExpectedLimit(double nsig, bool lower, const char *opt="") const
get expected limit (lower/upper) depending on the flag for asymptotic is a special case (the distribu...
SamplingDistribution * GetExpectedPValueDist(int index) const
return expected distribution of p-values (Cls or Clsplusb)
double CLsplusb(int index) const
return the observed CLsplusb value for the i-th entry
SamplingDistribution * GetLimitDistribution(bool lower) const
get the limit distribution (lower/upper depending on the flag) by interpolating the expected p values...
SamplingDistribution * GetNullTestStatDist(int index) const
same in terms of alt and null
SamplingDistribution * GetBackgroundTestStatDist(int index) const
get the background test statistic distribution
HypoTestResult is a base class for results from hypothesis tests.
void SetBackgroundAsAlt(Bool_t l=kTRUE)
Bool_t GetPValueIsRightTail(void) const
virtual void Append(const HypoTestResult *other)
add values from another HypoTestResult
virtual Double_t CLsplusb() const
Convert AlternatePValue into a "confidence level".
void SetPValueIsRightTail(Bool_t pr)
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)
void SetNullDistribution(SamplingDistribution *null)
virtual Double_t CLs() const
is simply (not a method, but a quantity)
Bool_t GetBackGroundIsAlt(void) const
void SetAltDistribution(SamplingDistribution *alt)
SamplingDistribution * GetNullDistribution(void) const
SamplingDistribution * GetAltDistribution(void) const
This class simply holds a sampling distribution of some test statistic.
Int_t GetSize() const
size of samples
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
SimpleInterval is a concrete implementation of the ConfInterval interface.
virtual Double_t ConfidenceLevel() const
return confidence level
SimpleInterval & operator=(const SimpleInterval &other)
The Canvas class.
Definition: TCanvas.h:31
virtual void RemoveAll(TCollection *col)
Remove all objects in collection col from this collection.
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
1-Dim function class
Definition: TF1.h:211
virtual Double_t Derivative(Double_t x, Double_t *params=0, Double_t epsilon=0.001) const
Returns the first derivative of the function at point x, computed by Richardson's extrapolation metho...
Definition: TF1.cxx:1097
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
Definition: TF1.cxx:3497
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:638
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
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseGeneralPurpose, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3980
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
virtual void Add(TObject *obj)
Definition: TList.h:87
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:354
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition: TList.cxx:467
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TList.cxx:656
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
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TObject.cxx:144
virtual TObject * RemoveAt(Int_t idx)
Basic string class.
Definition: TString.h:131
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
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
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
return c1
Definition: legend1.C:41
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Template specialisation used in RooAbsArg:
@ InputArguments
Definition: RooGlobalFunc.h:58
Namespace for the RooStats classes.
Definition: Asimov.h:20
static constexpr double s
TMath.
Definition: TMathBase.h:35
Bool_t IsNaN(Double_t x)
Definition: TMath.h:880
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
Definition: TMathBase.h:337
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:889
Double_t Floor(Double_t x)
Definition: TMath.h:691
T MinElement(Long64_t n, const T *a)
Return minimum of array a of length n.
Definition: TMath.h:940
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:416
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition: TMath.h:412
T MaxElement(Long64_t n, const T *a)
Return maximum of array a of length n.
Definition: TMath.h:947
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMathBase.h:278
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=0, Int_t type=7)
Computes sample quantiles, corresponding to the given probabilities.
Definition: TMath.cxx:1190
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition: TMath.h:902
Definition: file.py:1
Definition: graph.py:1
auto * m
Definition: textangle.C:8
static long int sum(long int i)
Definition: Factory.cxx:2258