Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
54
55////////////////////////////////////////////////////////////////////////////////
56/// default constructor
57
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/// Remove problematic points from this result.
193///
194/// This function can be used to clean up a result that has failed fits, spiking CLs
195/// or similar problems. It removes
196/// - Points where CLs is not falling monotonously. These may result from a lack of numerical precision.
197/// - Points where CLs spikes to more than 0.999.
198/// - Points with very low CLs. These are not needed to run the inverter, which speeds up the process.
199/// - Points where CLs < 0. These occur when fits fail.
201{
202 const int nEntries = ArraySize();
203
204 // initialization
205 double nsig1(1.0);
206 double nsig2(2.0);
207 double p[5];
208 double q[5];
209
210 p[0] = ROOT::Math::normal_cdf(-nsig2);
211 p[1] = ROOT::Math::normal_cdf(-nsig1);
212 p[2] = 0.5;
213 p[3] = ROOT::Math::normal_cdf(nsig1);
214 p[4] = ROOT::Math::normal_cdf(nsig2);
215
216 bool resultIsAsymptotic(false);
217 if (nEntries>=1) {
218 HypoTestResult* r = dynamic_cast<HypoTestResult *> ( GetResult(0) );
219 assert(r!=0);
220 if ( !r->GetNullDistribution() && !r->GetAltDistribution() ) {
221 resultIsAsymptotic = true;
222 }
223 }
224
225 int nPointsRemoved(0);
226
227 double CLsobsprev(1.0);
228
229 for (auto itr = fXValues.begin(); itr != fXValues.end(); ++itr) {
230 const double x = *itr;
231 const int i = FindIndex(x);
232
234 if (!s) break;
235
236 /////////////////////////////////////////////////////////////////////////////////////////
237
238 const std::vector<double> & values = s->GetSamplingDistribution();
239 if ((int) values.size() != fgAsymptoticNumPoints) {
240 oocoutE(this,Eval) << "HypoTestInverterResult::ExclusionCleanup - invalid size of sampling distribution" << std::endl;
241 delete s;
242 break;
243 }
244
245 /// expected p-values
246 // special case for asymptotic results (cannot use TMath::quantile in that case)
247 if (resultIsAsymptotic) {
248 double maxSigma = fgAsymptoticMaxSigma;
249 double dsig = 2.*maxSigma / (values.size() -1) ;
250 int i0 = (int) TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5 );
251 int i1 = (int) TMath::Floor ( ( -nsig1 + maxSigma )/dsig + 0.5 );
252 int i2 = (int) TMath::Floor ( ( maxSigma )/dsig + 0.5 );
253 int i3 = (int) TMath::Floor ( ( nsig1 + maxSigma )/dsig + 0.5 );
254 int i4 = (int) TMath::Floor ( ( nsig2 + maxSigma )/dsig + 0.5 );
255 //
256 q[0] = values[i0];
257 q[1] = values[i1];
258 q[2] = values[i2];
259 q[3] = values[i3];
260 q[4] = values[i4];
261 } else {
262 double * z = const_cast<double *>( &values[0] ); // need to change TMath::Quantiles
263 TMath::Quantiles(values.size(), 5, z, q, p, false);
264 }
265
266 delete s;
267
268 const double CLsobs = CLs(i);
269
270 /////////////////////////////////////////////////////////////////////////////////////////
271
272 bool removeThisPoint(false);
273
274 // 1. CLs should drop, else skip this point
275 if (resultIsAsymptotic && i>=1 && CLsobs>CLsobsprev) {
276 removeThisPoint = true;
277 } else if (CLsobs >= 0.) {
278 CLsobsprev = CLsobs;
279 }
280
281 // 2. CLs should not spike, else skip this point
282 removeThisPoint |= i>=1 && CLsobs >= 0.9999;
283
284 // 3. Not interested in CLs values that become too low.
285 removeThisPoint |= i>=1 && q[4] < fCLsCleanupThreshold;
286
287 // 4. Negative CLs indicate failed fits
288 removeThisPoint |= CLsobs < 0.;
289
290 // to remove or not to remove
291 if (removeThisPoint) {
292 itr = fXValues.erase(itr)--;
295 nPointsRemoved++;
296 continue;
297 } else { // keep
298 CLsobsprev = CLsobs;
299 }
300 }
301
302 // after cleanup, reset existing limits
303 fFittedUpperLimit = false;
304 fFittedLowerLimit = false;
306
307 return nPointsRemoved;
308}
309
310////////////////////////////////////////////////////////////////////////////////
311/// Merge this HypoTestInverterResult with another
312/// HypoTestInverterResult passed as argument
313/// The merge is done by combining the HypoTestResult when the same point value exist in both results.
314/// If results exist at different points these are added in the new result
315/// NOTE: Merging of the expected p-values obtained with pseudo-data.
316/// When expected p-values exist in the result (i.e. when rebuild option is used when getting the expected
317/// limit distribution in the HYpoTestInverter) then the expected p-values are also merged. This is equivalent
318/// at merging the pseudo-data. However there can be an inconsistency if the expected p-values have been
319/// obtained with different toys. In this case the merge is done but a warning message is printed.
320
322{
323 int nThis = ArraySize();
324 int nOther = otherResult.ArraySize();
325 if (nOther == 0) return true;
326 if (nOther != otherResult.fYObjects.GetSize() ) return false;
327 if (nThis != fYObjects.GetSize() ) return false;
328
329 // cannot merge in case of inconsistent members
330 if (fExpPValues.GetSize() > 0 && fExpPValues.GetSize() != nThis ) return false;
331 if (otherResult.fExpPValues.GetSize() > 0 && otherResult.fExpPValues.GetSize() != nOther ) return false;
332
333 oocoutI(this,Eval) << "HypoTestInverterResult::Add - merging result from " << otherResult.GetName()
334 << " in " << GetName() << std::endl;
335
336 bool addExpPValues = (fExpPValues.GetSize() == 0 && otherResult.fExpPValues.GetSize() > 0);
337 bool mergeExpPValues = (fExpPValues.GetSize() > 0 && otherResult.fExpPValues.GetSize() > 0);
338
339 if (addExpPValues || mergeExpPValues)
340 oocoutI(this,Eval) << "HypoTestInverterResult::Add - merging also the expected p-values from pseudo-data" << std::endl;
341
342
343 // case current result is empty
344 // just make a simple copy of the other result
345 if (nThis == 0) {
346 fXValues = otherResult.fXValues;
347 for (int i = 0; i < nOther; ++i)
348 fYObjects.Add( otherResult.fYObjects.At(i)->Clone() );
349 for (int i = 0; i < fExpPValues.GetSize() ; ++i)
350 fExpPValues.Add( otherResult.fExpPValues.At(i)->Clone() );
351 }
352 // now do the real merge combining point with same value or adding extra ones
353 else {
354 for (int i = 0; i < nOther; ++i) {
355 double otherVal = otherResult.fXValues[i];
356 HypoTestResult * otherHTR = (HypoTestResult*) otherResult.fYObjects.At(i);
357 if (otherHTR == 0) continue;
358 bool sameXFound = false;
359 for (int j = 0; j < nThis; ++j) {
360 double thisVal = fXValues[j];
361
362 // if same value merge the result
363 if ( (std::abs(otherVal) < 1 && TMath::AreEqualAbs(otherVal, thisVal,1.E-12) ) ||
364 (std::abs(otherVal) >= 1 && TMath::AreEqualRel(otherVal, thisVal,1.E-12) ) ) {
365 HypoTestResult * thisHTR = (HypoTestResult*) fYObjects.At(j);
366 thisHTR->Append(otherHTR);
367 sameXFound = true;
368 if (mergeExpPValues) {
370 // check if same toys have been used for the test statistic distribution
371 int thisNToys = (thisHTR->GetNullDistribution() ) ? thisHTR->GetNullDistribution()->GetSize() : 0;
372 int otherNToys = (otherHTR->GetNullDistribution() ) ? otherHTR->GetNullDistribution()->GetSize() : 0;
373 if (thisNToys != otherNToys )
374 oocoutW(this,Eval) << "HypoTestInverterResult::Add expected p values have been generated with different toys " << thisNToys << " , " << otherNToys << std::endl;
375 }
376 break;
377 }
378 }
379 if (!sameXFound) {
380 // add the new result
381 fYObjects.Add(otherHTR->Clone() );
382 fXValues.push_back( otherVal);
383 }
384 // add in any case also when same x found
385 if (addExpPValues)
386 fExpPValues.Add( otherResult.fExpPValues.At(i)->Clone() );
387
388
389 }
390 }
391
392 if (ArraySize() > nThis)
393 oocoutI(this,Eval) << "HypoTestInverterResult::Add - new number of points is " << fXValues.size()
394 << std::endl;
395 else
396 oocoutI(this,Eval) << "HypoTestInverterResult::Add - new toys/point is "
397 << ((HypoTestResult*) fYObjects.At(0))->GetNullDistribution()->GetSize()
398 << std::endl;
399
400 // reset cached limit values
403
404 return true;
405}
406
407////////////////////////////////////////////////////////////////////////////////
408/// Add a single point result (an HypoTestResult)
409
411{
412 int i= FindIndex(x);
413 if (i<0) {
414 fXValues.push_back(x);
415 fYObjects.Add(res.Clone());
416 } else {
418 if (!r) return false;
419 r->Append(&res);
420 }
421
422 // reset cached limit values
425
426 return true;
427}
428
429////////////////////////////////////////////////////////////////////////////////
430/// function to return the value of the parameter of interest for the i^th entry in the results
431
432double HypoTestInverterResult::GetXValue( int index ) const
433{
434 if ( index >= ArraySize() || index<0 ) {
435 oocoutE(this,InputArguments) << "Problem: You are asking for an impossible array index value\n";
436 return -999;
437 }
438
439 return fXValues[index];
440}
441
442////////////////////////////////////////////////////////////////////////////////
443/// function to return the value of the confidence level for the i^th entry in the results
444
445double HypoTestInverterResult::GetYValue( int index ) const
446{
447 auto result = GetResult(index);
448 if ( !result ) {
449 return -999;
450 }
451
452 if (fUseCLs) {
453 return result->CLs();
454 } else {
455 return result->CLsplusb(); // CLs+b
456 }
457}
458
459////////////////////////////////////////////////////////////////////////////////
460/// function to return the estimated error on the value of the confidence level for the i^th entry in the results
461
462double HypoTestInverterResult::GetYError( int index ) const
463{
464 auto result = GetResult(index);
465 if ( !result ) {
466 return -999;
467 }
468
469 if (fUseCLs) {
470 return result->CLsError();
471 } else {
472 return result->CLsplusbError();
473 }
474}
475
476////////////////////////////////////////////////////////////////////////////////
477/// function to return the observed CLb value for the i-th entry
478
479double HypoTestInverterResult::CLb( int index ) const
480{
481 auto result = GetResult(index);
482 if ( !result ) {
483 return -999;
484 }
485 return result->CLb();
486}
487
488////////////////////////////////////////////////////////////////////////////////
489/// function to return the observed CLs+b value for the i-th entry
490
491double HypoTestInverterResult::CLsplusb( int index ) const
492{
493 auto result = GetResult(index);
494 if ( !result) {
495 return -999;
496 }
497 return result->CLsplusb();
498}
499
500////////////////////////////////////////////////////////////////////////////////
501/// function to return the observed CLs value for the i-th entry
502
503double HypoTestInverterResult::CLs( int index ) const
504{
505 auto result = GetResult(index);
506 if ( !result ) {
507 return -999;
508 }
509 return result->CLs();
510}
511
512////////////////////////////////////////////////////////////////////////////////
513/// function to return the error on the observed CLb value for the i-th entry
514
515double HypoTestInverterResult::CLbError( int index ) const
516{
517 auto result = GetResult(index);
518 if ( !result ) {
519 return -999;
520 }
521 return result->CLbError();
522}
523
524////////////////////////////////////////////////////////////////////////////////
525/// function to return the error on the observed CLs+b value for the i-th entry
526
528{
529 auto result = GetResult(index);
530 if ( ! result ) {
531 return -999;
532 }
533 return result->CLsplusbError();
534}
535
536////////////////////////////////////////////////////////////////////////////////
537/// function to return the error on the observed CLs value for the i-th entry
538
539double HypoTestInverterResult::CLsError( int index ) const
540{
541 auto result = GetResult(index);
542 if ( ! result ){
543 return -999;
544 }
545 return result->CLsError();
546}
547
548////////////////////////////////////////////////////////////////////////////////
549/// get the HypoTestResult object at the given index point
550
552{
553 if ( index >= ArraySize() || index<0 ) {
554 oocoutE(this,InputArguments) << "Problem: You are asking for an impossible array index value\n";
555 return 0;
556 }
557
558 return ((HypoTestResult*) fYObjects.At(index));
559}
560
561////////////////////////////////////////////////////////////////////////////////
562/// find the index corresponding at the poi value xvalue
563/// If no points is found return -1
564/// Note that a tolerance is used of 10^-12 to find the closest point
565
566int HypoTestInverterResult::FindIndex(double xvalue) const
567{
568 const double tol = 1.E-12;
569 for (int i=0; i<ArraySize(); i++) {
570 double xpoint = fXValues[i];
571 if ( (std::abs(xvalue) > 1 && TMath::AreEqualRel( xvalue, xpoint, tol) ) ||
572 (std::abs(xvalue) <= 1 && TMath::AreEqualAbs( xvalue, xpoint, tol) ) )
573 return i;
574 }
575 return -1;
576}
577
578////////////////////////////////////////////////////////////////////////////////
579/// return the X value of the given graph for the target value y0
580/// the graph is evaluated using linear interpolation by default.
581/// if option = "S" a TSpline3 is used
582
583double HypoTestInverterResult::GetGraphX(const TGraph & graph, double y0, bool lowSearch, double &axmin, double &axmax) const {
584
585//#define DO_DEBUG
586#ifdef DO_DEBUG
587 std::cout << "using graph for search " << lowSearch << " min " << axmin << " max " << axmax << std::endl;
588#endif
589
590
591 // find reasonable xmin and xmax if not given
592 const double * y = graph.GetY();
593 int n = graph.GetN();
594 if (n < 2) {
595 ooccoutE(this,Eval) << "HypoTestInverterResult::GetGraphX - need at least 2 points for interpolation (n=" << n << ")\n";
596 return (n>0) ? y[0] : 0;
597 }
598
599 double varmin = - TMath::Infinity();
600 double varmax = TMath::Infinity();
601 const RooRealVar* var = dynamic_cast<RooRealVar*>( fParameters.first() );
602 if (var) {
603 varmin = var->getMin();
604 varmax = var->getMax();
605 }
606
607
608 // find ymin and ymax and corresponding values
609 double ymin = TMath::MinElement(n,y);
610 double ymax = TMath::MaxElement(n,y);
611 // cannot find intercept in the full range - return min or max valie
612 if (ymax < y0) {
613 return (lowSearch) ? varmax : varmin;
614 }
615 if (ymin > y0) {
616 return (lowSearch) ? varmin : varmax;
617 }
618
619 double xmin = axmin;
620 double xmax = axmax;
621
622 // case no range is given check if need to extrapolate to lower/upper values
623 if (axmin >= axmax ) {
624
625#ifdef DO_DEBUG
626 std::cout << "No rage given - check if extrapolation is needed " << std::endl;
627#endif
628
629 xmin = graph.GetX()[0];
630 xmax = graph.GetX()[n-1];
631
632 double yfirst = graph.GetY()[0];
633 double ylast = graph.GetY()[n-1];
634
635 // distinguish the case we have lower /upper limits
636 // check if a possible crossing exists otherwise return variable min/max
637
638 // do lower extrapolation
639 if ( (ymax < y0 && !lowSearch) || ( yfirst > y0 && lowSearch) ) {
640 xmin = varmin;
641 }
642 // do upper extrapolation
643 if ( (ymax < y0 && lowSearch) || ( ylast > y0 && !lowSearch) ) {
644 xmax = varmax;
645 }
646 }
647
648 auto func = [&](double x) {
649 return (fInterpolOption == kSpline) ? graph.Eval(x, nullptr, "S") - y0 : graph.Eval(x) - y0;
650 };
651 ROOT::Math::Functor1D f1d(func);
652
654 brf.SetFunction(f1d,xmin,xmax);
655 brf.SetNpx(TMath::Max(graph.GetN()*2,100) );
656#ifdef DO_DEBUG
657 std::cout << "findind root for " << xmin << " , "<< xmax << "f(x) : " << graph.Eval(xmin) << " , " << graph.Eval(0.5*(xmax+xmin))
658 << " , " << graph.Eval(xmax) << " target " << y0 << std::endl;
659#endif
660
661 bool ret = brf.Solve(100, 1.E-16, 1.E-6);
662 if (!ret) {
663 ooccoutE(this,Eval) << "HypoTestInverterResult - interpolation failed for interval [" << xmin << "," << xmax
664 << " ] g(xmin,xmax) =" << graph.Eval(xmin) << "," << graph.Eval(xmax)
665 << " target=" << y0 << " return inf" << std::endl
666 << "One may try to clean up invalid points using HypoTestInverterResult::ExclusionCleanup()." << std::endl;
667 return TMath::Infinity();
668 }
669 double limit = brf.Root();
670
671#ifdef DO_DEBUG
672 if (lowSearch) std::cout << "lower limit search : ";
673 else std::cout << "Upper limit search : ";
674 std::cout << "interpolation done between " << xmin << " and " << xmax
675 << "\n Found limit using RootFinder is " << limit << std::endl;
676
677 TString fname = "graph_upper.root";
678 if (lowSearch) fname = "graph_lower.root";
679 auto file = TFile::Open(fname,"RECREATE");
680 graph.Write("graph");
681 file->Close();
682#endif
683
684 // look in case if a new intersection exists
685 // only when boundaries are not given
686 if (axmin >= axmax) {
687 int index = TMath::BinarySearch(n, graph.GetX(), limit);
688#ifdef DO_DEBUG
689 std::cout << "do new interpolation dividing from " << index << " and " << y[index] << std::endl;
690#endif
691
692 if (lowSearch && index >= 1 && (y[0] - y0) * ( y[index]- y0) < 0) {
693 //search if another intersection exists at a lower value
694 limit = GetGraphX(graph, y0, lowSearch, graph.GetX()[0], graph.GetX()[index] );
695 }
696 else if (!lowSearch && index < n-2 && (y[n-1] - y0) * ( y[index+1]- y0) < 0) {
697 // another intersection exists at an higher value
698 limit = GetGraphX(graph, y0, lowSearch, graph.GetX()[index+1], graph.GetX()[n-1] );
699 }
700 }
701 // return also xmin, xmax values
702 axmin = xmin;
703 axmax = xmax;
704
705 return limit;
706}
707
708////////////////////////////////////////////////////////////////////////////////
709/// interpolate to find a limit value
710/// Use a linear or a spline interpolation depending on the interpolation option
711
712double HypoTestInverterResult::FindInterpolatedLimit(double target, bool lowSearch, double xmin, double xmax )
713{
714
715 // variable minimum and maximum
716 double varmin = - TMath::Infinity();
717 double varmax = TMath::Infinity();
718 const RooRealVar* var = dynamic_cast<RooRealVar*>( fParameters.first() );
719 if (var) {
720 varmin = var->getMin();
721 varmax = var->getMax();
722 }
723
724 if (ArraySize()<2) {
725 double val = (lowSearch) ? xmin : xmax;
726 oocoutW(this,Eval) << "HypoTestInverterResult::FindInterpolatedLimit"
727 << " - not enough points to get the inverted interval - return "
728 << val << std::endl;
729 fLowerLimit = varmin;
730 fUpperLimit = varmax;
731 return (lowSearch) ? fLowerLimit : fUpperLimit;
732 }
733
734 // sort the values in x
735 int n = ArraySize();
736 std::vector<unsigned int> index(n );
737 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
738 // make a graph with the sorted point
739 TGraph graph(n);
740 for (int i = 0; i < n; ++i)
741 graph.SetPoint(i, GetXValue(index[i]), GetYValue(index[i] ) );
742
743
744 //std::cout << " search for " << lowSearch << " xmin = " << xmin << " xmax " << xmax << std::endl;
745
746
747 // search first for min/max in the given range
748 if (xmin >= xmax) {
749
750
751 // search for maximum between the point
752 double * itrmax = std::max_element(graph.GetY() , graph.GetY() +n);
753 double ymax = *itrmax;
754 int iymax = itrmax - graph.GetY();
755 double xwithymax = graph.GetX()[iymax];
756
757#ifdef DO_DEBUG
758 std::cout << " max of y " << iymax << " " << xwithymax << " " << ymax << " target is " << target << std::endl;
759#endif
760 // look if maximum is above/below target
761 if (ymax > target) {
762 if (lowSearch) {
763 if ( iymax > 0) {
764 // low search (minimum is first point or minimum range)
765 xmin = ( graph.GetY()[0] <= target ) ? graph.GetX()[0] : varmin;
766 xmax = xwithymax;
767 }
768 else {
769 // no room for lower limit
770 fLowerLimit = varmin;
771 return fLowerLimit;
772 }
773 }
774 if (!lowSearch ) {
775 // up search
776 if ( iymax < n-1 ) {
777 xmin = xwithymax;
778 xmax = ( graph.GetY()[n-1] <= target ) ? graph.GetX()[n-1] : varmax;
779 }
780 else {
781 // no room for upper limit
782 fUpperLimit = varmax;
783 return fUpperLimit;
784 }
785 }
786 }
787 else {
788 // in case is below the target
789 // find out if is a lower or upper search
790 if (iymax <= (n-1)/2 ) {
791 lowSearch = false;
792 fLowerLimit = varmin;
793 }
794 else {
795 lowSearch = true;
796 fUpperLimit = varmax;
797 }
798 }
799
800#ifdef DO_DEBUG
801 std::cout << " found xmin, xmax = " << xmin << " " << xmax << " for search " << lowSearch << std::endl;
802#endif
803
804 // now come here if I have already found a lower/upper limit
805 // i.e. I am calling routine for the second time
806#ifdef ISNEEDED
807 // should not really come here
808 if (lowSearch && fUpperLimit < varmax) {
809 xmin = fXValues[ index.front() ];
810 // find xmax (is first point before upper limit)
811 int upI = FindClosestPointIndex(target, 2, fUpperLimit);
812 if (upI < 1) return xmin;
813 xmax = GetXValue(upI);
814 }
815 else if (!lowSearch && fLowerLimit > varmin ) {
816 // find xmin (is first point after lower limit)
817 int lowI = FindClosestPointIndex(target, 3, fLowerLimit);
818 if (lowI >= n-1) return xmax;
819 xmin = GetXValue(lowI);
820 xmax = fXValues[ index.back() ];
821 }
822#endif
823 }
824
825#ifdef DO_DEBUG
826 std::cout << "finding " << lowSearch << " limit between " << xmin << " " << xmax << endl;
827#endif
828
829 // compute noe the limit using the TGraph interpolations routine
830 double limit = GetGraphX(graph, target, lowSearch, xmin, xmax);
831 if (lowSearch) fLowerLimit = limit;
832 else fUpperLimit = limit;
833 // estimate the error
834 double error = CalculateEstimatedError( target, lowSearch, xmin, xmax);
835
836 TString limitType = (lowSearch) ? "lower" : "upper";
837 ooccoutD(this,Eval) << "HypoTestInverterResult::FindInterpolateLimit "
838 << "the computed " << limitType << " limit is " << limit << " +/- " << error << std::endl;
839
840#ifdef DO_DEBUG
841 std::cout << "Found limit is " << limit << " +/- " << error << std::endl;
842#endif
843
844
845 if (lowSearch) return fLowerLimit;
846 return fUpperLimit;
847
848
849// if (lowSearch && !TMath::IsNaN(fUpperLimit)) return fLowerLimit;
850// if (!lowSearch && !TMath::IsNaN(fLowerLimit)) return fUpperLimit;
851// // is this needed ?
852// // we call again the function for the upper limits
853
854// // now perform the opposite search on the complement interval
855// if (lowSearch) {
856// xmin = xmax;
857// xmax = varmax;
858// } else {
859// xmax = xmin;
860// xmin = varmin;
861// }
862// double limit2 = GetGraphX(graph, target, !lowSearch, xmin, xmax);
863// if (!lowSearch) fLowerLimit = limit2;
864// else fUpperLimit = limit2;
865
866// CalculateEstimatedError( target, !lowSearch, xmin, xmax);
867
868// #ifdef DO_DEBUG
869// std::cout << "other limit is " << limit2 << std::endl;
870// #endif
871
872// return (lowSearch) ? fLowerLimit : fUpperLimit;
873
874}
875
876////////////////////////////////////////////////////////////////////////////////
877/// - if mode = 0
878/// find closest point to target in Y, the object closest to the target which is 3 sigma from the target
879/// and has smaller error
880/// - if mode = 1
881/// find 2 closest point to target in X and between these two take the one closer to the target
882/// - if mode = 2 as in mode = 1 but return the lower point not the closest one
883/// - if mode = 3 as in mode = 1 but return the upper point not the closest one
884
885int HypoTestInverterResult::FindClosestPointIndex(double target, int mode, double xtarget)
886{
887
888 int bestIndex = -1;
889 int closestIndex = -1;
890 if (mode == 0) {
891 double smallestError = 2; // error must be < 1
892 double bestValue = 2;
893 for (int i=0; i<ArraySize(); i++) {
894 double dist = fabs(GetYValue(i)-target);
895 if ( dist <3 *GetYError(i) ) { // less than 1 sigma from target CL
896 if (GetYError(i) < smallestError ) {
897 smallestError = GetYError(i);
898 bestIndex = i;
899 }
900 }
901 if ( dist < bestValue) {
902 bestValue = dist;
903 closestIndex = i;
904 }
905 }
906 if (bestIndex >=0) return bestIndex;
907 // if no points found just return the closest one to the target
908 return closestIndex;
909 }
910 // else mode = 1,2,3
911 // find the two closest points to limit value
912 // sort the array first
913 int n = fXValues.size();
914 std::vector<unsigned int> indx(n);
915 TMath::SortItr(fXValues.begin(), fXValues.end(), indx.begin(), false);
916 std::vector<double> xsorted( n);
917 for (int i = 0; i < n; ++i) xsorted[i] = fXValues[indx[i] ];
918 int index1 = TMath::BinarySearch( n, &xsorted[0], xtarget);
919
920#ifdef DO_DEBUG
921 std::cout << "finding closest point to " << xtarget << " is " << index1 << " " << indx[index1] << std::endl;
922#endif
923
924 // case xtarget is outside the range (before or afterwards)
925 if (index1 < 0) return indx[0];
926 if (index1 >= n-1) return indx[n-1];
927 int index2 = index1 +1;
928
929 if (mode == 2) return (GetXValue(indx[index1]) < GetXValue(indx[index2])) ? indx[index1] : indx[index2];
930 if (mode == 3) return (GetXValue(indx[index1]) > GetXValue(indx[index2])) ? indx[index1] : indx[index2];
931 // get smaller point of the two (mode == 1)
932 if (fabs(GetYValue(indx[index1])-target) <= fabs(GetYValue(indx[index2])-target) )
933 return indx[index1];
934 return indx[index2];
935
936}
937
938////////////////////////////////////////////////////////////////////////////////
939
941{
942 if (fFittedLowerLimit) return fLowerLimit;
943 //std::cout << "finding point with cl = " << 1-(1-ConfidenceLevel())/2 << endl;
945 // find both lower/upper limit
947 } else {
948 //LM: I think this is never called
950 }
951 return fLowerLimit;
952}
953
954////////////////////////////////////////////////////////////////////////////////
955
957{
958 //std::cout << "finding point with cl = " << (1-ConfidenceLevel())/2 << endl;
959 if (fFittedUpperLimit) return fUpperLimit;
962 } else {
963 //LM: I think this is never called
965 }
966 return fUpperLimit;
967}
968
969////////////////////////////////////////////////////////////////////////////////
970/// Return an error estimate on the upper(lower) limit. This is the error on
971/// either CLs or CLsplusb divided by an estimate of the slope at this
972/// point.
973
974Double_t HypoTestInverterResult::CalculateEstimatedError(double target, bool lower, double xmin, double xmax)
975{
976
977 if (ArraySize()==0) {
978 oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimateError"
979 << "Empty result \n";
980 return 0;
981 }
982
983 if (ArraySize()<2) {
984 oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimateError"
985 << " only points - return its error\n";
986 return GetYError(0);
987 }
988
989 // it does not make sense in case of asymptotic which do not have point errors
990 if (!GetNullTestStatDist(0) ) return 0;
991
992 TString type = (!lower) ? "upper" : "lower";
993
994#ifdef DO_DEBUG
995 std::cout << "calculate estimate error " << type << " between " << xmin << " and " << xmax << std::endl;
996 std::cout << "computed limit is " << ( (lower) ? fLowerLimit : fUpperLimit ) << std::endl;
997#endif
998
999 // make a TGraph Errors with the sorted points
1000 std::vector<unsigned int> indx(fXValues.size());
1001 TMath::SortItr(fXValues.begin(), fXValues.end(), indx.begin(), false);
1002 // make a graph with the sorted point
1004 int ip = 0, np = 0;
1005 for (int i = 0; i < ArraySize(); ++i) {
1006 if ( (xmin < xmax) && ( GetXValue(indx[i]) >= xmin && GetXValue(indx[i]) <= xmax) ) {
1007 np++;
1008 // exclude points with zero or very small errors
1009 if (GetYError(indx[i] ) > 1.E-6) {
1010 graph.SetPoint(ip, GetXValue(indx[i]), GetYValue(indx[i] ) );
1011 graph.SetPointError(ip, 0., GetYError(indx[i]) );
1012 ip++;
1013 }
1014 }
1015 }
1016 if (graph.GetN() < 2) {
1017 if (np >= 2) oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimatedError - no valid points - cannot estimate the " << type << " limit error " << std::endl;
1018 return 0;
1019 }
1020
1021 double minX = xmin;
1022 double maxX = xmax;
1023 if (xmin >= xmax) {
1024 minX = fXValues[ indx.front() ];
1025 maxX = fXValues[ indx.back() ];
1026 }
1027
1028 TF1 fct("fct", "exp([0] * (x - [2] ) + [1] * (x-[2])**2)", minX, maxX);
1029 double scale = maxX-minX;
1030 if (lower) {
1031 fct.SetParameters( 2./scale, 0.1/scale, graph.GetX()[0] );
1032 fct.SetParLimits(0,0,100./scale);
1033 fct.SetParLimits(1,0, 10./scale); }
1034 else {
1035 fct.SetParameters( -2./scale, -0.1/scale );
1036 fct.SetParLimits(0,-100./scale, 0);
1037 fct.SetParLimits(1,-100./scale, 0); }
1038
1039 if (graph.GetN() < 3) fct.FixParameter(1,0.);
1040
1041 // find the point closest to the limit
1042 double limit = (!lower) ? fUpperLimit : fLowerLimit;
1043 if (TMath::IsNaN(limit)) return 0; // cannot do if limit not computed
1044
1045
1046#ifdef DO_DEBUG
1047 TCanvas * c1 = new TCanvas();
1048 std::cout << "fitting for limit " << type << "between " << minX << " , " << maxX << " points considered " << graph.GetN() << std::endl;
1049 int fitstat = graph.Fit(&fct," EX0");
1050 graph.SetMarkerStyle(20);
1051 graph.Draw("AP");
1052 graph.Print();
1053 c1->SaveAs(TString::Format("graphFit_%s.pdf",type.Data()) );
1054 delete c1;
1055#else
1056 int fitstat = graph.Fit(&fct,"Q EX0");
1057#endif
1058
1059 int index = FindClosestPointIndex(target, 1, limit);
1060 double theError = 0;
1061 if (fitstat == 0) {
1062 double errY = GetYError(index);
1063 if (errY > 0) {
1064 double m = fct.Derivative( GetXValue(index) );
1065 theError = std::min(fabs( GetYError(index) / m), maxX-minX);
1066 }
1067 }
1068 else {
1069 oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimatedError - cannot estimate the " << type << " limit error " << std::endl;
1070 theError = 0;
1071 }
1072 if (lower)
1073 fLowerLimitError = theError;
1074 else
1075 fUpperLimitError = theError;
1076
1077#ifdef DO_DEBUG
1078 std::cout << "closes point to the limit is " << index << " " << GetXValue(index) << " and has error " << GetYError(index) << std::endl;
1079#endif
1080
1081 return theError;
1082}
1083
1084////////////////////////////////////////////////////////////////////////////////
1085/// need to have compute first lower limit
1086
1088{
1090 if (fLowerLimitError >= 0) return fLowerLimitError;
1091 // try to recompute the error
1092 return CalculateEstimatedError(1-ConfidenceLevel(), true);
1093}
1094
1095////////////////////////////////////////////////////////////////////////////////
1096
1098{
1100 if (fUpperLimitError >= 0) return fUpperLimitError;
1101 // try to recompute the error
1102 return CalculateEstimatedError(1-ConfidenceLevel(), false);
1103}
1104
1105////////////////////////////////////////////////////////////////////////////////
1106/// get the background test statistic distribution
1107
1109
1110 HypoTestResult * firstResult = (HypoTestResult*) fYObjects.At(index);
1111 if (!firstResult) return 0;
1112 return firstResult->GetBackGroundIsAlt() ? firstResult->GetAltDistribution() : firstResult->GetNullDistribution();
1113}
1114
1115////////////////////////////////////////////////////////////////////////////////
1116/// get the signal and background test statistic distribution
1117
1119 HypoTestResult * result = (HypoTestResult*) fYObjects.At(index);
1120 if (!result) return 0;
1121 return !result->GetBackGroundIsAlt() ? result->GetAltDistribution() : result->GetNullDistribution();
1122}
1123
1124////////////////////////////////////////////////////////////////////////////////
1125/// get the expected p-value distribution at the scanned point index
1126
1128
1129 if (index < 0 || index >= ArraySize() ) return 0;
1130
1131 if (fExpPValues.GetSize() == ArraySize() ) {
1132 return (SamplingDistribution*) fExpPValues.At(index)->Clone();
1133 }
1134
1135 static bool useFirstB = false;
1136 // get the S+B distribution
1137 int bIndex = (useFirstB) ? 0 : index;
1138
1139 SamplingDistribution * bDistribution = GetBackgroundTestStatDist(bIndex);
1141
1142 HypoTestResult * result = (HypoTestResult*) fYObjects.At(index);
1143
1144 if (bDistribution && sbDistribution) {
1145
1146 // create a new HypoTestResult
1147 HypoTestResult tempResult;
1148 tempResult.SetPValueIsRightTail( result->GetPValueIsRightTail() );
1149 tempResult.SetBackgroundAsAlt( true);
1150 // ownership of SamplingDistribution is in HypoTestResult class now
1151 tempResult.SetNullDistribution( new SamplingDistribution(*sbDistribution) );
1152 tempResult.SetAltDistribution( new SamplingDistribution(*bDistribution ) );
1153
1154 std::vector<double> values(bDistribution->GetSize());
1155 for (int i = 0; i < bDistribution->GetSize(); ++i) {
1156 tempResult.SetTestStatisticData( bDistribution->GetSamplingDistribution()[i] );
1157 values[i] = (fUseCLs) ? tempResult.CLs() : tempResult.CLsplusb();
1158 }
1159 return new SamplingDistribution("expected values","expected values",values);
1160 }
1161 // in case b abs sbDistribution are null assume is coming from the asymptotic calculator
1162 // hard -coded this value (no really needed to be used by user)
1165 const double smax = fgAsymptoticMaxSigma;
1166 const int npoints = fgAsymptoticNumPoints;
1167 const double dsig = 2* smax/ (npoints-1) ;
1168 std::vector<double> values(npoints);
1169 for (int i = 0; i < npoints; ++i) {
1170 double nsig = -smax + dsig*i;
1171 double pval = AsymptoticCalculator::GetExpectedPValues( result->NullPValue(), result->AlternatePValue(), nsig, fUseCLs, !fIsTwoSided);
1172 if (pval < 0) { return 0;}
1173 values[i] = pval;
1174 }
1175 return new SamplingDistribution("Asymptotic expected values","Asymptotic expected values",values);
1176
1177}
1178
1179////////////////////////////////////////////////////////////////////////////////
1180/// get the limit distribution (lower/upper depending on the flag)
1181/// by interpolating the expected p values for each point
1182
1184 if (ArraySize()<2) {
1185 oocoutE(this,Eval) << "HypoTestInverterResult::GetLimitDistribution"
1186 << " not enough points - return 0 " << std::endl;
1187 return 0;
1188 }
1189
1190 ooccoutD(this,Eval) << "HypoTestInverterResult - computing limit distribution...." << std::endl;
1191
1192
1193 // find optimal size by looking at the PValue distribution obtained
1194 int npoints = ArraySize();
1195 std::vector<SamplingDistribution*> distVec( npoints );
1196 double sum = 0;
1197 for (unsigned int i = 0; i < distVec.size(); ++i) {
1198 distVec[i] = GetExpectedPValueDist(i);
1199 // sort the distributions
1200 // hack (by calling InverseCDF(0) will sort the sampling distribution
1201 if (distVec[i] ) {
1202 distVec[i]->InverseCDF(0);
1203 sum += distVec[i]->GetSize();
1204 }
1205 }
1206 int size = int( sum/ npoints);
1207
1208 if (size < 10) {
1209 ooccoutW(this,InputArguments) << "HypoTestInverterResult - set a minimum size of 10 for limit distribution" << std::endl;
1210 size = 10;
1211 }
1212
1213
1214 double target = 1-fConfidenceLevel;
1215
1216 // vector with the quantiles of the p-values for each scanned poi point
1217 std::vector< std::vector<double> > quantVec(npoints );
1218 for (int i = 0; i < npoints; ++i) {
1219
1220 if (!distVec[i]) continue;
1221
1222 // make quantiles from the sampling distributions of the expected p values
1223 std::vector<double> pvalues = distVec[i]->GetSamplingDistribution();
1224 delete distVec[i]; distVec[i] = 0;
1225 std::sort(pvalues.begin(), pvalues.end());
1226 // find the quantiles of the distribution
1227 double p[1] = {0};
1228 double q[1] = {0};
1229
1230 quantVec[i] = std::vector<double>(size);
1231 for (int ibin = 0; ibin < size; ++ibin) {
1232 // exclude for a bug in TMath::Quantiles last item
1233 p[0] = std::min( (ibin+1) * 1./double(size), 1.0);
1234 // use the type 1 which give the point value
1235 TMath::Quantiles(pvalues.size(), 1, &pvalues[0], q, p, true, 0, 1 );
1236 (quantVec[i])[ibin] = q[0];
1237 }
1238
1239 }
1240
1241 // sort the values in x
1242 std::vector<unsigned int> index( npoints );
1243 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1244
1245 // SamplingDistribution * dist0 = distVec[index[0]];
1246 // SamplingDistribution * dist1 = distVec[index[1]];
1247
1248
1249 std::vector<double> limits(size);
1250 // loop on the p values and find the limit for each expected point in the quantiles vector
1251 for (int j = 0; j < size; ++j ) {
1252
1253 TGraph g;
1254 for (int k = 0; k < npoints ; ++k) {
1255 if (quantVec[index[k]].size() > 0 )
1256 g.SetPoint(g.GetN(), GetXValue(index[k]), (quantVec[index[k]])[j] );
1257 }
1258
1259 limits[j] = GetGraphX(g, target, lower);
1260
1261 }
1262
1263
1264 if (lower)
1265 return new SamplingDistribution("Expected lower Limit","Expected lower limits",limits);
1266 else
1267 return new SamplingDistribution("Expected upper Limit","Expected upper limits",limits);
1268
1269}
1270
1271////////////////////////////////////////////////////////////////////////////////
1272/// Get the expected lower limit
1273/// nsig is used to specify which expected value of the UpperLimitDistribution
1274/// For example
1275/// - nsig = 0 (default value) returns the expected value
1276/// - nsig = -1 returns the lower band value at -1 sigma
1277/// - nsig + 1 return the upper value
1278/// - opt = "" (default) : compute limit by interpolating all the p values, find the corresponding limit distribution
1279/// and then find the quantiles in the limit distribution
1280/// ioption = "P" is the method used for plotting. One Finds the corresponding nsig quantile in the p values and then
1281/// interpolates them
1282
1283double HypoTestInverterResult::GetExpectedLowerLimit(double nsig, const char * opt ) const {
1284
1285 return GetExpectedLimit(nsig, true, opt);
1286}
1287
1288////////////////////////////////////////////////////////////////////////////////
1289/// Get the expected upper limit
1290/// nsig is used to specify which expected value of the UpperLimitDistribution
1291/// For example
1292/// - nsig = 0 (default value) returns the expected value
1293/// - nsig = -1 returns the lower band value at -1 sigma
1294/// - nsig + 1 return the upper value
1295/// - opt is an option specifying the type of method used for computing the upper limit
1296/// - opt = "" (default) : compute limit by interpolating all the p values, find the corresponding limit distribution
1297/// and then find the quantiles in the limit distribution
1298/// ioption = "P" is the method used for plotting. One Finds the corresponding nsig quantile in the p values and then
1299/// interpolates them
1300
1301double HypoTestInverterResult::GetExpectedUpperLimit(double nsig, const char * opt ) const {
1302
1303 return GetExpectedLimit(nsig, false, opt);
1304}
1305
1306////////////////////////////////////////////////////////////////////////////////
1307/// get expected limit (lower/upper) depending on the flag
1308/// for asymptotic is a special case (the distribution is generated an step in sigma values)
1309/// distinguish asymptotic looking at the hypotest results
1310/// if option = "P" get expected limit using directly quantiles of p value distribution
1311/// else (default) find expected limit by obtaining first a full limit distributions
1312/// The last one is in general more correct
1313
1314double HypoTestInverterResult::GetExpectedLimit(double nsig, bool lower, const char * opt ) const {
1315
1316 const int nEntries = ArraySize();
1317 if (nEntries <= 0) return (lower) ? 1 : 0; // return 1 for lower, 0 for upper
1318
1319 HypoTestResult * r = dynamic_cast<HypoTestResult *> (fYObjects.First() );
1320 assert(r != 0);
1321 if (!r->GetNullDistribution() && !r->GetAltDistribution() ) {
1322 // we are in the asymptotic case
1323 // get the limits obtained at the different sigma values
1324 SamplingDistribution * limitDist = GetLimitDistribution(lower);
1325 if (!limitDist) return 0;
1326 const std::vector<double> & values = limitDist->GetSamplingDistribution();
1327 if (values.size() <= 1) return 0;
1328 double dsig = 2* fgAsymptoticMaxSigma/ (values.size() -1) ;
1329 int i = (int) TMath::Floor ( (nsig + fgAsymptoticMaxSigma)/dsig + 0.5);
1330 return values[i];
1331 }
1332
1333 double p[1] = {0};
1334 double q[1] = {0};
1335 p[0] = ROOT::Math::normal_cdf(nsig,1);
1336
1337 // for CLs+b can get the quantiles of p-value distribution and
1338 // interpolate them
1339 // In the case of CLs (since it is not a real p-value anymore but a ratio)
1340 // then it is needed to get first all limit distribution values and then the quantiles
1341
1342 TString option(opt);
1343 option.ToUpper();
1344 if (option.Contains("P")) {
1345
1346 TGraph g;
1347
1348 // sort the arrays based on the x values
1349 std::vector<unsigned int> index(nEntries);
1350 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1351
1352 for (int j=0; j<nEntries; ++j) {
1353 int i = index[j]; // i is the order index
1355 if (!s) {
1356 ooccoutI(this,Eval) << "HypoTestInverterResult - cannot compute expected p value distribution for point, x = "
1357 << GetXValue(i) << " skip it " << std::endl;
1358 continue;
1359 }
1360 const std::vector<double> & values = s->GetSamplingDistribution();
1361 double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
1362 TMath::Quantiles(values.size(), 1, x,q,p,false);
1363 g.SetPoint(g.GetN(), fXValues[i], q[0] );
1364 delete s;
1365 }
1366 if (g.GetN() < 2) {
1367 ooccoutE(this,Eval) << "HypoTestInverterResult - cannot compute limits , not enough points, n = " << g.GetN() << std::endl;
1368 return 0;
1369 }
1370
1371 // interpolate the graph to obtain the limit
1372 double target = 1-fConfidenceLevel;
1373 return GetGraphX(g, target, lower);
1374
1375 }
1376 // here need to use the limit distribution
1377 SamplingDistribution * limitDist = GetLimitDistribution(lower);
1378 if (!limitDist) return 0;
1379 const std::vector<double> & values = limitDist->GetSamplingDistribution();
1380 double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
1381 TMath::Quantiles(values.size(), 1, x,q,p,false);
1382 return q[0];
1383
1384}
ROOT::R::TRInterface & r
Definition Object.C:4
#define g(i)
Definition RSha256.hxx:105
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define oocoutW(o, a)
#define oocoutE(o, a)
#define oocoutI(o, a)
#define ooccoutW(o, a)
#define ooccoutI(o, a)
#define ooccoutD(o, a)
#define ooccoutE(o, a)
#define ClassImp(name)
Definition Rtypes.h:364
char name[80]
Definition TGX11.cxx:110
int type
Definition TGX11.cxx:121
float xmin
float * q
float ymin
float xmax
float ymax
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:506
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Add an argument and transfer the ownership to the collection.
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.
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
virtual TObject * clone(const char *newname) const
Definition RooRealVar.h:51
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:23
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.
1-Dim function class
Definition TF1.h:213
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:1108
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
Definition TF1.cxx:3518
virtual void SetParameters(const Double_t *params)
Definition TF1.h:644
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:1553
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:4025
A TGraphErrors is a TGraph with error bars.
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2298
virtual void Add(TObject *obj)
Definition TList.h:81
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition TList.cxx:357
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition TList.cxx:470
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition TList.cxx:659
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:216
virtual TObject * RemoveAt(Int_t idx)
Basic string class.
Definition TString.h:136
void ToUpper()
Change string to upper case.
Definition TString.cxx:1163
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:2336
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:624
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
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition Common.h:18
@ InputArguments
Namespace for the RooStats classes.
Definition Asimov.h:19
TMath.
Definition TMathBase.h:35
Bool_t IsNaN(Double_t x)
Definition TMath.h:842
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:208
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
Definition TMathBase.h:333
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition TMath.h:851
Double_t Floor(Double_t x)
Definition TMath.h:653
T MinElement(Long64_t n, const T *a)
Return minimum of array a of length n.
Definition TMath.h:902
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition TMath.h:430
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition TMath.h:424
T MaxElement(Long64_t n, const T *a)
Return maximum of array a of length n.
Definition TMath.h:909
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition TMathBase.h:274
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:1183
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition TMath.h:864
Definition file.py:1
Definition graph.py:1
auto * m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345