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 fXValues(other.fXValues)
94{
97 int nOther = other.ArraySize();
98
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)
174 fParameters.addOwned(std::unique_ptr<RooRealVar>{static_cast<RooRealVar *>(scannedVariable.clone(scannedVariable.GetName()))});
175}
176
177////////////////////////////////////////////////////////////////////////////////
178/// destructor
179
181{
182 // explicitly empty the TLists - these contain pointers, not objects
185
188}
189
190////////////////////////////////////////////////////////////////////////////////
191/// Remove problematic points from this result.
192///
193/// This function can be used to clean up a result that has failed fits, spiking CLs
194/// or similar problems. It removes
195/// - Points where CLs is not falling monotonously. These may result from a lack of numerical precision.
196/// - Points where CLs spikes to more than 0.999.
197/// - Points with very low CLs. These are not needed to run the inverter, which speeds up the process.
198/// - Points where CLs < 0. These occur when fits fail.
200{
201 const int nEntries = ArraySize();
202
203 // initialization
204 double nsig1(1.0);
205 double nsig2(2.0);
206 double p[5];
207 double q[5];
208
209 p[0] = ROOT::Math::normal_cdf(-nsig2);
210 p[1] = ROOT::Math::normal_cdf(-nsig1);
211 p[2] = 0.5;
212 p[3] = ROOT::Math::normal_cdf(nsig1);
213 p[4] = ROOT::Math::normal_cdf(nsig2);
214
215 bool resultIsAsymptotic(false);
216 if (nEntries>=1) {
217 HypoTestResult* r = dynamic_cast<HypoTestResult *> ( GetResult(0) );
218 assert(r!=nullptr);
219 if ( !r->GetNullDistribution() && !r->GetAltDistribution() ) {
220 resultIsAsymptotic = true;
221 }
222 }
223
224 int nPointsRemoved(0);
225
226 double CLsobsprev(1.0);
227
228 for (auto itr = fXValues.begin(); itr != fXValues.end(); ++itr) {
229 const double x = *itr;
230 const int i = FindIndex(x);
231
233 if (!s) break;
234
235 /////////////////////////////////////////////////////////////////////////////////////////
236
237 const std::vector<double> & values = s->GetSamplingDistribution();
238 if ((int) values.size() != fgAsymptoticNumPoints) {
239 coutE(Eval) << "HypoTestInverterResult::ExclusionCleanup - invalid size of sampling distribution" << std::endl;
240 delete s;
241 break;
242 }
243
244 /// expected p-values
245 // special case for asymptotic results (cannot use TMath::quantile in that case)
246 if (resultIsAsymptotic) {
247 double maxSigma = fgAsymptoticMaxSigma;
248 double dsig = 2.*maxSigma / (values.size() -1) ;
249 int i0 = (int) TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5 );
250 int i1 = (int) TMath::Floor ( ( -nsig1 + maxSigma )/dsig + 0.5 );
251 int i2 = (int) TMath::Floor ( ( maxSigma )/dsig + 0.5 );
252 int i3 = (int) TMath::Floor ( ( nsig1 + maxSigma )/dsig + 0.5 );
253 int i4 = (int) TMath::Floor ( ( nsig2 + maxSigma )/dsig + 0.5 );
254 //
255 q[0] = values[i0];
256 q[1] = values[i1];
257 q[2] = values[i2];
258 q[3] = values[i3];
259 q[4] = values[i4];
260 } else {
261 double * z = const_cast<double *>(&values[0] ); // need to change TMath::Quantiles
262 TMath::Quantiles(values.size(), 5, z, q, p, false);
263 }
264
265 delete s;
266
267 const double CLsobs = CLs(i);
268
269 /////////////////////////////////////////////////////////////////////////////////////////
270
271 bool removeThisPoint(false);
272
273 // 1. CLs should drop, else skip this point
274 if (resultIsAsymptotic && i>=1 && CLsobs>CLsobsprev) {
275 removeThisPoint = true;
276 } else if (CLsobs >= 0.) {
277 CLsobsprev = CLsobs;
278 }
279
280 // 2. CLs should not spike, else skip this point
281 removeThisPoint |= i>=1 && CLsobs >= 0.9999;
282
283 // 3. Not interested in CLs values that become too low.
284 removeThisPoint |= i>=1 && q[4] < fCLsCleanupThreshold;
285
286 // 4. Negative CLs indicate failed fits
287 removeThisPoint |= CLsobs < 0.;
288
289 // to remove or not to remove
290 if (removeThisPoint) {
291 itr = fXValues.erase(itr)--;
294 nPointsRemoved++;
295 continue;
296 } else { // keep
297 CLsobsprev = CLsobs;
298 }
299 }
300
301 // after cleanup, reset existing limits
302 fFittedUpperLimit = false;
303 fFittedLowerLimit = false;
305
306 return nPointsRemoved;
307}
308
309////////////////////////////////////////////////////////////////////////////////
310/// Merge this HypoTestInverterResult with another
311/// HypoTestInverterResult passed as argument
312/// The merge is done by combining the HypoTestResult when the same point value exist in both results.
313/// If results exist at different points these are added in the new result
314/// NOTE: Merging of the expected p-values obtained with pseudo-data.
315/// When expected p-values exist in the result (i.e. when rebuild option is used when getting the expected
316/// limit distribution in the HYpoTestInverter) then the expected p-values are also merged. This is equivalent
317/// at merging the pseudo-data. However there can be an inconsistency if the expected p-values have been
318/// obtained with different toys. In this case the merge is done but a warning message is printed.
319
321{
322 int nThis = ArraySize();
323 int nOther = otherResult.ArraySize();
324 if (nOther == 0) return true;
325 if (nOther != otherResult.fYObjects.GetSize() ) return false;
326 if (nThis != fYObjects.GetSize() ) return false;
327
328 // cannot merge in case of inconsistent members
329 if (fExpPValues.GetSize() > 0 && fExpPValues.GetSize() != nThis ) return false;
330 if (otherResult.fExpPValues.GetSize() > 0 && otherResult.fExpPValues.GetSize() != nOther ) return false;
331
332 coutI(Eval) << "HypoTestInverterResult::Add - merging result from " << otherResult.GetName()
333 << " in " << GetName() << std::endl;
334
335 bool addExpPValues = (fExpPValues.GetSize() == 0 && otherResult.fExpPValues.GetSize() > 0);
336 bool mergeExpPValues = (fExpPValues.GetSize() > 0 && otherResult.fExpPValues.GetSize() > 0);
337
338 if (addExpPValues || mergeExpPValues)
339 coutI(Eval) << "HypoTestInverterResult::Add - merging also the expected p-values from pseudo-data" << std::endl;
340
341
342 // case current result is empty
343 // just make a simple copy of the other result
344 if (nThis == 0) {
345 fXValues = otherResult.fXValues;
346 for (int i = 0; i < nOther; ++i)
347 fYObjects.Add( otherResult.fYObjects.At(i)->Clone() );
348 for (int i = 0; i < fExpPValues.GetSize() ; ++i)
349 fExpPValues.Add( otherResult.fExpPValues.At(i)->Clone() );
350 }
351 // now do the real merge combining point with same value or adding extra ones
352 else {
353 for (int i = 0; i < nOther; ++i) {
354 double otherVal = otherResult.fXValues[i];
355 HypoTestResult * otherHTR = static_cast<HypoTestResult*>(otherResult.fYObjects.At(i));
356 if (otherHTR == nullptr) continue;
357 bool sameXFound = false;
358 for (int j = 0; j < nThis; ++j) {
359 double thisVal = fXValues[j];
360
361 // if same value merge the result
362 if ( (std::abs(otherVal) < 1 && TMath::AreEqualAbs(otherVal, thisVal,1.E-12) ) ||
363 (std::abs(otherVal) >= 1 && TMath::AreEqualRel(otherVal, thisVal,1.E-12) ) ) {
364 HypoTestResult * thisHTR = static_cast<HypoTestResult*>(fYObjects.At(j));
365 thisHTR->Append(otherHTR);
366 sameXFound = true;
367 if (mergeExpPValues) {
368 (static_cast<SamplingDistribution*>(fExpPValues.At(j)))->Add( static_cast<SamplingDistribution*>(otherResult.fExpPValues.At(i)) );
369 // check if same toys have been used for the test statistic distribution
370 int thisNToys = (thisHTR->GetNullDistribution() ) ? thisHTR->GetNullDistribution()->GetSize() : 0;
371 int otherNToys = (otherHTR->GetNullDistribution() ) ? otherHTR->GetNullDistribution()->GetSize() : 0;
372 if (thisNToys != otherNToys )
373 coutW(Eval) << "HypoTestInverterResult::Add expected p values have been generated with different toys " << thisNToys << " , " << otherNToys << std::endl;
374 }
375 break;
376 }
377 }
378 if (!sameXFound) {
379 // add the new result
380 fYObjects.Add(otherHTR->Clone() );
381 fXValues.push_back( otherVal);
382 }
383 // add in any case also when same x found
384 if (addExpPValues)
385 fExpPValues.Add( otherResult.fExpPValues.At(i)->Clone() );
386
387
388 }
389 }
390
391 if (ArraySize() > nThis) {
392 coutI(Eval) << "HypoTestInverterResult::Add - new number of points is " << fXValues.size()
393 << std::endl;
394 } else {
395 coutI(Eval) << "HypoTestInverterResult::Add - new toys/point is "
396 << (static_cast<HypoTestResult *>(fYObjects.At(0)))->GetNullDistribution()->GetSize() << std::endl;
397 }
398
399 // reset cached limit values
402
403 return true;
404}
405
406////////////////////////////////////////////////////////////////////////////////
407/// Add a single point result (an HypoTestResult)
408
410{
411 int i= FindIndex(x);
412 if (i<0) {
413 fXValues.push_back(x);
414 fYObjects.Add(res.Clone());
415 } else {
417 if (!r) return false;
418 r->Append(&res);
419 }
420
421 // reset cached limit values
424
425 return true;
426}
427
428////////////////////////////////////////////////////////////////////////////////
429/// function to return the value of the parameter of interest for the i^th entry in the results
430
432{
433 if ( index >= ArraySize() || index<0 ) {
434 coutE(InputArguments) << "Problem: You are asking for an impossible array index value\n";
435 return -999;
436 }
437
438 return fXValues[index];
439}
440
441////////////////////////////////////////////////////////////////////////////////
442/// function to return the value of the confidence level for the i^th entry in the results
443
445{
446 auto result = GetResult(index);
447 if ( !result ) {
448 return -999;
449 }
450
451 if (fUseCLs) {
452 return result->CLs();
453 } else {
454 return result->CLsplusb(); // CLs+b
455 }
456}
457
458////////////////////////////////////////////////////////////////////////////////
459/// function to return the estimated error on the value of the confidence level for the i^th entry in the results
460
462{
463 auto result = GetResult(index);
464 if ( !result ) {
465 return -999;
466 }
467
468 if (fUseCLs) {
469 return result->CLsError();
470 } else {
471 return result->CLsplusbError();
472 }
473}
474
475////////////////////////////////////////////////////////////////////////////////
476/// function to return the observed CLb value for the i-th entry
477
479{
480 auto result = GetResult(index);
481 if ( !result ) {
482 return -999;
483 }
484 return result->CLb();
485}
486
487////////////////////////////////////////////////////////////////////////////////
488/// function to return the observed CLs+b value for the i-th entry
489
491{
492 auto result = GetResult(index);
493 if ( !result) {
494 return -999;
495 }
496 return result->CLsplusb();
497}
498
499////////////////////////////////////////////////////////////////////////////////
500/// function to return the observed CLs value for the i-th entry
501
503{
504 auto result = GetResult(index);
505 if ( !result ) {
506 return -999;
507 }
508 return result->CLs();
509}
510
511////////////////////////////////////////////////////////////////////////////////
512/// function to return the error on the observed CLb value for the i-th entry
513
515{
516 auto result = GetResult(index);
517 if ( !result ) {
518 return -999;
519 }
520 return result->CLbError();
521}
522
523////////////////////////////////////////////////////////////////////////////////
524/// function to return the error on the observed CLs+b value for the i-th entry
525
527{
528 auto result = GetResult(index);
529 if ( ! result ) {
530 return -999;
531 }
532 return result->CLsplusbError();
533}
534
535////////////////////////////////////////////////////////////////////////////////
536/// function to return the error on the observed CLs value for the i-th entry
537
539{
540 auto result = GetResult(index);
541 if ( ! result ){
542 return -999;
543 }
544 return result->CLsError();
545}
546
547////////////////////////////////////////////////////////////////////////////////
548/// get the HypoTestResult object at the given index point
549
551{
552 if ( index >= ArraySize() || index<0 ) {
553 coutE(InputArguments) << "Problem: You are asking for an impossible array index value\n";
554 return nullptr;
555 }
556
557 return (static_cast<HypoTestResult*>(fYObjects.At(index)));
558}
559
560////////////////////////////////////////////////////////////////////////////////
561/// find the index corresponding at the poi value xvalue
562/// If no points is found return -1
563/// Note that a tolerance is used of 10^-12 to find the closest point
564
565int HypoTestInverterResult::FindIndex(double xvalue) const
566{
567 const double tol = 1.E-12;
568 for (int i=0; i<ArraySize(); i++) {
569 double xpoint = fXValues[i];
570 if ( (std::abs(xvalue) > 1 && TMath::AreEqualRel( xvalue, xpoint, tol) ) ||
571 (std::abs(xvalue) <= 1 && TMath::AreEqualAbs( xvalue, xpoint, tol) ) )
572 return i;
573 }
574 return -1;
575}
576
577////////////////////////////////////////////////////////////////////////////////
578/// return the X value of the given graph for the target value y0
579/// the graph is evaluated using linear interpolation by default.
580/// if option = "S" a TSpline3 is used
581
582double HypoTestInverterResult::GetGraphX(const TGraph & graph, double y0, bool lowSearch, double &axmin, double &axmax) const {
583
584//#define DO_DEBUG
585#ifdef DO_DEBUG
586 std::cout << "using graph for search " << lowSearch << " min " << axmin << " max " << axmax << std::endl;
587#endif
588
589
590 // find reasonable xmin and xmax if not given
591 const double * y = graph.GetY();
592 int n = graph.GetN();
593 if (n < 2) {
594 ooccoutE(this,Eval) << "HypoTestInverterResult::GetGraphX - need at least 2 points for interpolation (n=" << n << ")\n";
595 return (n>0) ? y[0] : 0;
596 }
597
598 double varmin = - TMath::Infinity();
599 double varmax = TMath::Infinity();
600 const RooRealVar* var = dynamic_cast<RooRealVar*>(fParameters.first() );
601 if (var) {
602 varmin = var->getMin();
603 varmax = var->getMax();
604 }
605
606
607 // find ymin and ymax and corresponding values
608 double ymin = TMath::MinElement(n,y);
609 double ymax = TMath::MaxElement(n,y);
610 // cannot find intercept in the full range - return min or max valie
611 if (ymax < y0) {
612 return (lowSearch) ? varmax : varmin;
613 }
614 if (ymin > y0) {
615 return (lowSearch) ? varmin : varmax;
616 }
617
618 double xmin = axmin;
619 double xmax = axmax;
620
621 // case no range is given check if need to extrapolate to lower/upper values
622 if (axmin >= axmax ) {
623
624#ifdef DO_DEBUG
625 std::cout << "No range given - check if extrapolation is needed " << std::endl;
626#endif
627
628 xmin = graph.GetX()[0];
629 xmax = graph.GetX()[n-1];
630
631 double yfirst = graph.GetY()[0];
632 double ylast = graph.GetY()[n-1];
633
634 // distinguish the case we have lower /upper limits
635 // check if a possible crossing exists otherwise return variable min/max
636
637 // do lower extrapolation
638 if ( (ymax < y0 && !lowSearch) || ( yfirst > y0 && lowSearch) ) {
639 xmin = varmin;
640 }
641 // do upper extrapolation
642 if ( (ymax < y0 && lowSearch) || ( ylast > y0 && !lowSearch) ) {
643 xmax = varmax;
644 }
645 }
646
647 auto func = [&](double x) {
648 return (fInterpolOption == kSpline) ? graph.Eval(x, nullptr, "S") - y0 : graph.Eval(x) - y0;
649 };
650 ROOT::Math::Functor1D f1d(func);
651
653 brf.SetFunction(f1d,xmin,xmax);
654 brf.SetNpx(TMath::Max(graph.GetN()*2,100) );
655#ifdef DO_DEBUG
656 std::cout << "findind root for " << xmin << " , "<< xmax << "f(x) : " << graph.Eval(xmin) << " , " << graph.Eval(0.5*(xmax+xmin))
657 << " , " << graph.Eval(xmax) << " target " << y0 << std::endl;
658#endif
659
660 bool ret = brf.Solve(100, 1.E-16, 1.E-6);
661 if (!ret) {
662 ooccoutE(this,Eval) << "HypoTestInverterResult - interpolation failed for interval [" << xmin << "," << xmax
663 << " ] g(xmin,xmax) =" << graph.Eval(xmin) << "," << graph.Eval(xmax)
664 << " target=" << y0 << " return inf" << std::endl
665 << "One may try to clean up invalid points using HypoTestInverterResult::ExclusionCleanup()." << std::endl;
666 return TMath::Infinity();
667 }
668 double limit = brf.Root();
669
670#ifdef DO_DEBUG
671 if (lowSearch) std::cout << "lower limit search : ";
672 else std::cout << "Upper limit search : ";
673 std::cout << "interpolation done between " << xmin << " and " << xmax
674 << "\n Found limit using RootFinder is " << limit << std::endl;
675
676 TString fname = "graph_upper.root";
677 if (lowSearch) fname = "graph_lower.root";
678 {
679 std::unique_ptr<TFile> file{TFile::Open(fname,"RECREATE")};
680 graph.Write("graph");
681 }
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 coutW(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)
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)
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 now 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
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 = std::abs(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 (std::abs(GetYValue(indx[index1])-target) <= std::abs(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 HypoTestInverterResult::CalculateEstimatedError(double target, bool lower, double xmin, double xmax)
975{
976
977 if (ArraySize()==0) {
978 coutW(Eval) << "HypoTestInverterResult::CalculateEstimateError"
979 << "Empty result \n";
980 return 0;
981 }
982
983 if (ArraySize()<2) {
984 coutW(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;
1005 int np = 0;
1006 for (int i = 0; i < ArraySize(); ++i) {
1007 if ( (xmin < xmax) && ( GetXValue(indx[i]) >= xmin && GetXValue(indx[i]) <= xmax) ) {
1008 np++;
1009 // exclude points with zero or very small errors
1010 if (GetYError(indx[i] ) > 1.E-6) {
1011 graph.SetPoint(ip, GetXValue(indx[i]), GetYValue(indx[i] ) );
1012 graph.SetPointError(ip, 0., GetYError(indx[i]) );
1013 ip++;
1014 }
1015 }
1016 }
1017 if (graph.GetN() < 2) {
1018 if (np >= 2) coutW(Eval) << "HypoTestInverterResult::CalculateEstimatedError - no valid points - cannot estimate the " << type << " limit error " << std::endl;
1019 return 0;
1020 }
1021
1022 double minX = xmin;
1023 double maxX = xmax;
1024 if (xmin >= xmax) {
1025 minX = fXValues[ indx.front() ];
1026 maxX = fXValues[ indx.back() ];
1027 }
1028
1029 TF1 fct("fct", "exp([0] * (x - [2] ) + [1] * (x-[2])**2)", minX, maxX);
1030 double scale = maxX-minX;
1031 if (lower) {
1032 fct.SetParameters( 2./scale, 0.1/scale, graph.GetX()[0] );
1033 fct.SetParLimits(0,0,100./scale);
1034 fct.SetParLimits(1,0, 10./scale); }
1035 else {
1036 fct.SetParameters( -2./scale, -0.1/scale );
1037 fct.SetParLimits(0,-100./scale, 0);
1038 fct.SetParLimits(1,-100./scale, 0); }
1039
1040 if (graph.GetN() < 3) fct.FixParameter(1,0.);
1041
1042 // find the point closest to the limit
1043 double limit = (!lower) ? fUpperLimit : fLowerLimit;
1044 if (TMath::IsNaN(limit)) return 0; // cannot do if limit not computed
1045
1046
1047#ifdef DO_DEBUG
1048 TCanvas * c1 = new TCanvas();
1049 std::cout << "fitting for limit " << type << "between " << minX << " , " << maxX << " points considered " << graph.GetN() << std::endl;
1050 int fitstat = graph.Fit(&fct," EX0");
1051 graph.SetMarkerStyle(20);
1052 graph.Draw("AP");
1053 graph.Print();
1054 c1->SaveAs(TString::Format("graphFit_%s.pdf",type.Data()) );
1055 delete c1;
1056#else
1057 int fitstat = graph.Fit(&fct,"Q EX0");
1058#endif
1059
1060 int index = FindClosestPointIndex(target, 1, limit);
1061 double theError = 0;
1062 if (fitstat == 0) {
1063 double errY = GetYError(index);
1064 if (errY > 0) {
1065 double m = fct.Derivative( GetXValue(index) );
1066 theError = std::min(std::abs( GetYError(index) / m), maxX-minX);
1067 }
1068 }
1069 else {
1070 coutW(Eval) << "HypoTestInverterResult::CalculateEstimatedError - cannot estimate the " << type << " limit error " << std::endl;
1071 theError = 0;
1072 }
1073 if (lower) {
1074 fLowerLimitError = theError;
1075 } else {
1076 fUpperLimitError = theError;
1077 }
1078
1079#ifdef DO_DEBUG
1080 std::cout << "closes point to the limit is " << index << " " << GetXValue(index) << " and has error " << GetYError(index) << std::endl;
1081#endif
1082
1083 return theError;
1084}
1085
1086////////////////////////////////////////////////////////////////////////////////
1087/// need to have compute first lower limit
1088
1090{
1092 if (fLowerLimitError >= 0) return fLowerLimitError;
1093 // try to recompute the error
1094 return CalculateEstimatedError(1-ConfidenceLevel(), true);
1095}
1096
1097////////////////////////////////////////////////////////////////////////////////
1098
1100{
1102 if (fUpperLimitError >= 0) return fUpperLimitError;
1103 // try to recompute the error
1104 return CalculateEstimatedError(1-ConfidenceLevel(), false);
1105}
1106
1107////////////////////////////////////////////////////////////////////////////////
1108/// get the background test statistic distribution
1109
1111
1112 HypoTestResult * firstResult = static_cast<HypoTestResult*>(fYObjects.At(index));
1113 if (!firstResult) return nullptr;
1114 return firstResult->GetBackGroundIsAlt() ? firstResult->GetAltDistribution() : firstResult->GetNullDistribution();
1115}
1116
1117////////////////////////////////////////////////////////////////////////////////
1118/// get the signal and background test statistic distribution
1119
1122 if (!result) return nullptr;
1123 return !result->GetBackGroundIsAlt() ? result->GetAltDistribution() : result->GetNullDistribution();
1124}
1125
1126////////////////////////////////////////////////////////////////////////////////
1127/// get the expected p-value distribution at the scanned point index
1128
1130
1131 if (index < 0 || index >= ArraySize() ) return nullptr;
1132
1133 if (fExpPValues.GetSize() == ArraySize() ) {
1134 return static_cast<SamplingDistribution*>( fExpPValues.At(index)->Clone());
1135 }
1136
1137 static bool useFirstB = false;
1138 // get the S+B distribution
1139 int bIndex = (useFirstB) ? 0 : index;
1140
1141 SamplingDistribution * bDistribution = GetBackgroundTestStatDist(bIndex);
1143
1145
1146 if (bDistribution && sbDistribution) {
1147
1148 // create a new HypoTestResult
1149 HypoTestResult tempResult;
1150 tempResult.SetPValueIsRightTail( result->GetPValueIsRightTail() );
1151 tempResult.SetBackgroundAsAlt( true);
1152 // ownership of SamplingDistribution is in HypoTestResult class now
1153 tempResult.SetNullDistribution( new SamplingDistribution(*sbDistribution) );
1154 tempResult.SetAltDistribution( new SamplingDistribution(*bDistribution ) );
1155
1156 std::vector<double> values(bDistribution->GetSize());
1157 for (int i = 0; i < bDistribution->GetSize(); ++i) {
1158 tempResult.SetTestStatisticData( bDistribution->GetSamplingDistribution()[i] );
1159 values[i] = (fUseCLs) ? tempResult.CLs() : tempResult.CLsplusb();
1160 }
1161 return new SamplingDistribution("expected values","expected values",values);
1162 }
1163 // in case b abs sbDistribution are null assume is coming from the asymptotic calculator
1164 // hard -coded this value (no really needed to be used by user)
1167 const double smax = fgAsymptoticMaxSigma;
1168 const int npoints = fgAsymptoticNumPoints;
1169 const double dsig = 2* smax/ (npoints-1) ;
1170 std::vector<double> values(npoints);
1171 for (int i = 0; i < npoints; ++i) {
1172 double nsig = -smax + dsig*i;
1173 double pval = AsymptoticCalculator::GetExpectedPValues( result->NullPValue(), result->AlternatePValue(), nsig, fUseCLs, !fIsTwoSided);
1174 if (pval < 0) { return nullptr;}
1175 values[i] = pval;
1176 }
1177 return new SamplingDistribution("Asymptotic expected values","Asymptotic expected values",values);
1178
1179}
1180
1181////////////////////////////////////////////////////////////////////////////////
1182/// get the limit distribution (lower/upper depending on the flag)
1183/// by interpolating the expected p values for each point
1184
1186 if (ArraySize()<2) {
1187 coutE(Eval) << "HypoTestInverterResult::GetLimitDistribution"
1188 << " not enough points - return 0 " << std::endl;
1189 return nullptr;
1190 }
1191
1192 ooccoutD(this,Eval) << "HypoTestInverterResult - computing limit distribution...." << std::endl;
1193
1194
1195 // find optimal size by looking at the PValue distribution obtained
1196 int npoints = ArraySize();
1197 std::vector<SamplingDistribution*> distVec( npoints );
1198 double sum = 0;
1199 for (unsigned int i = 0; i < distVec.size(); ++i) {
1200 distVec[i] = GetExpectedPValueDist(i);
1201 // sort the distributions
1202 // hack (by calling InverseCDF(0) will sort the sampling distribution
1203 if (distVec[i] ) {
1204 distVec[i]->InverseCDF(0);
1205 sum += distVec[i]->GetSize();
1206 }
1207 }
1208 int size = int( sum/ npoints);
1209
1210 if (size < 10) {
1211 ooccoutW(this,InputArguments) << "HypoTestInverterResult - set a minimum size of 10 for limit distribution" << std::endl;
1212 size = 10;
1213 }
1214
1215
1216 double target = 1-fConfidenceLevel;
1217
1218 // vector with the quantiles of the p-values for each scanned poi point
1219 std::vector< std::vector<double> > quantVec(npoints );
1220 for (int i = 0; i < npoints; ++i) {
1221
1222 if (!distVec[i]) continue;
1223
1224 // make quantiles from the sampling distributions of the expected p values
1225 std::vector<double> pvalues = distVec[i]->GetSamplingDistribution();
1226 delete distVec[i]; distVec[i] = nullptr;
1227 std::sort(pvalues.begin(), pvalues.end());
1228 // find the quantiles of the distribution
1229 double p[1] = {0};
1230 double q[1] = {0};
1231
1232 quantVec[i] = std::vector<double>(size);
1233 for (int ibin = 0; ibin < size; ++ibin) {
1234 // exclude for a bug in TMath::Quantiles last item
1235 p[0] = std::min( (ibin+1) * 1./double(size), 1.0);
1236 // use the type 1 which give the point value
1237 TMath::Quantiles(pvalues.size(), 1, &pvalues[0], q, p, true, nullptr, 1 );
1238 (quantVec[i])[ibin] = q[0];
1239 }
1240
1241 }
1242
1243 // sort the values in x
1244 std::vector<unsigned int> index( npoints );
1245 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1246
1247 // SamplingDistribution * dist0 = distVec[index[0]];
1248 // SamplingDistribution * dist1 = distVec[index[1]];
1249
1250
1251 std::vector<double> limits(size);
1252 // loop on the p values and find the limit for each expected point in the quantiles vector
1253 for (int j = 0; j < size; ++j ) {
1254
1255 TGraph g;
1256 for (int k = 0; k < npoints ; ++k) {
1257 if (!quantVec[index[k]].empty() )
1258 g.SetPoint(g.GetN(), GetXValue(index[k]), (quantVec[index[k]])[j] );
1259 }
1260
1261 limits[j] = GetGraphX(g, target, lower);
1262
1263 }
1264
1265 if (lower) {
1266 return new SamplingDistribution("Expected lower Limit","Expected lower limits",limits);
1267 } else {
1268 return new SamplingDistribution("Expected upper Limit", "Expected upper limits", limits);
1269 }
1270}
1271
1272////////////////////////////////////////////////////////////////////////////////
1273/// Get the expected lower limit
1274/// nsig is used to specify which expected value of the UpperLimitDistribution
1275/// For example
1276/// - nsig = 0 (default value) returns the expected value
1277/// - nsig = -1 returns the lower band value at -1 sigma
1278/// - nsig + 1 return the upper value
1279/// - opt = "" (default) : compute limit by interpolating all the p values, find the corresponding limit distribution
1280/// and then find the quantiles in the limit distribution
1281/// ioption = "P" is the method used for plotting. One Finds the corresponding nsig quantile in the p values and then
1282/// interpolates them
1283
1284double HypoTestInverterResult::GetExpectedLowerLimit(double nsig, const char * opt ) const {
1285
1286 return GetExpectedLimit(nsig, true, opt);
1287}
1288
1289////////////////////////////////////////////////////////////////////////////////
1290/// Get the expected upper limit
1291/// nsig is used to specify which expected value of the UpperLimitDistribution
1292/// For example
1293/// - nsig = 0 (default value) returns the expected value
1294/// - nsig = -1 returns the lower band value at -1 sigma
1295/// - nsig + 1 return the upper value
1296/// - opt is an option specifying the type of method used for computing the upper limit
1297/// - opt = "" (default) : compute limit by interpolating all the p values, find the corresponding limit distribution
1298/// and then find the quantiles in the limit distribution
1299/// ioption = "P" is the method used for plotting. One Finds the corresponding nsig quantile in the p values and then
1300/// interpolates them
1301
1302double HypoTestInverterResult::GetExpectedUpperLimit(double nsig, const char * opt ) const {
1303
1304 return GetExpectedLimit(nsig, false, opt);
1305}
1306
1307////////////////////////////////////////////////////////////////////////////////
1308/// get expected limit (lower/upper) depending on the flag
1309/// for asymptotic is a special case (the distribution is generated an step in sigma values)
1310/// distinguish asymptotic looking at the hypotest results
1311/// if option = "P" get expected limit using directly quantiles of p value distribution
1312/// else (default) find expected limit by obtaining first a full limit distributions
1313/// The last one is in general more correct
1314
1315double HypoTestInverterResult::GetExpectedLimit(double nsig, bool lower, const char * opt ) const {
1316
1317 const int nEntries = ArraySize();
1318 if (nEntries <= 0) return (lower) ? 1 : 0; // return 1 for lower, 0 for upper
1319
1320 HypoTestResult * r = dynamic_cast<HypoTestResult *> (fYObjects.First() );
1321 assert(r != nullptr);
1322 if (!r->GetNullDistribution() && !r->GetAltDistribution() ) {
1323 // we are in the asymptotic case
1324 // get the limits obtained at the different sigma values
1325 SamplingDistribution * limitDist = GetLimitDistribution(lower);
1326 if (!limitDist) return 0;
1327 const std::vector<double> & values = limitDist->GetSamplingDistribution();
1328 if (values.size() <= 1) return 0;
1329 double dsig = 2* fgAsymptoticMaxSigma/ (values.size() -1) ;
1330 int i = (int) TMath::Floor ( (nsig + fgAsymptoticMaxSigma)/dsig + 0.5);
1331 return values[i];
1332 }
1333
1334 double p[1] = {0};
1335 double q[1] = {0};
1336 p[0] = ROOT::Math::normal_cdf(nsig,1);
1337
1338 // for CLs+b can get the quantiles of p-value distribution and
1339 // interpolate them
1340 // In the case of CLs (since it is not a real p-value anymore but a ratio)
1341 // then it is needed to get first all limit distribution values and then the quantiles
1342
1343 TString option(opt);
1344 option.ToUpper();
1345 if (option.Contains("P")) {
1346
1347 TGraph g;
1348
1349 // sort the arrays based on the x values
1350 std::vector<unsigned int> index(nEntries);
1351 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1352
1353 for (int j=0; j<nEntries; ++j) {
1354 int i = index[j]; // i is the order index
1356 if (!s) {
1357 ooccoutI(this,Eval) << "HypoTestInverterResult - cannot compute expected p value distribution for point, x = "
1358 << GetXValue(i) << " skip it " << std::endl;
1359 continue;
1360 }
1361 const std::vector<double> & values = s->GetSamplingDistribution();
1362 double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
1363 TMath::Quantiles(values.size(), 1, x,q,p,false);
1364 g.SetPoint(g.GetN(), fXValues[i], q[0] );
1365 delete s;
1366 }
1367 if (g.GetN() < 2) {
1368 ooccoutE(this,Eval) << "HypoTestInverterResult - cannot compute limits , not enough points, n = " << g.GetN() << std::endl;
1369 return 0;
1370 }
1371
1372 // interpolate the graph to obtain the limit
1373 double target = 1-fConfidenceLevel;
1374 return GetGraphX(g, target, lower);
1375
1376 }
1377 // here need to use the limit distribution
1378 SamplingDistribution * limitDist = GetLimitDistribution(lower);
1379 if (!limitDist) return 0;
1380 const std::vector<double> & values = limitDist->GetSamplingDistribution();
1381 double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
1382 TMath::Quantiles(values.size(), 1, x,q,p,false);
1383 return q[0];
1384
1385}
#define g(i)
Definition RSha256.hxx:105
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
#define coutW(a)
#define coutE(a)
#define ooccoutW(o, a)
#define ooccoutI(o, a)
#define ooccoutD(o, a)
#define ooccoutE(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t target
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char mode
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
char name[80]
Definition TGX11.cxx:110
float xmin
float * q
float ymin
float xmax
float ymax
Class for finding the root of a one dimensional function using the Brent algorithm.
bool SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup) override
Sets the function for the rest of the algorithms.
bool Solve(int maxIter=100, double absTol=1E-8, double relTol=1E-10) override
Returns the X value corresponding to the function value fy for (xmin<x<xmax).
double Root() const override
Returns root value.
void SetNpx(int npx)
Set the number of point used to bracket root using a grid.
Functor1D class for one-dimensional functions.
Definition Functor.h:95
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
RooAbsArg * first() const
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
TObject * clone(const char *newname) const override
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 GetYValue(int index) const
function to return the value of the confidence level for the i^th entry in the results
double LowerLimitEstimatedError()
rough estimation of the error on the computed bound of the confidence interval Estimate of lower limi...
double FindInterpolatedLimit(double target, bool lowSearch=false, double xmin=1, double xmax=0.0)
interpolate to find a limit value Use a linear or a spline interpolation depending on the interpolati...
double UpperLimitEstimatedError()
Estimate of lower limit error function evaluates only a rough error on the lower limit.
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
get the signal and background test statistic distribution
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
HypoTestInverterResult & operator=(const HypoTestInverterResult &other)
operator =
InterpolOption_t fInterpolOption
interpolation option (linear or spline)
int ArraySize() const
number of entries in the results array
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
HypoTestInverterResult(const char *name=nullptr)
default constructor
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results
double UpperLimit() override
return the interval upper limit
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 expected sampling distribution 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
bool fIsTwoSided
two sided scan (look for lower/upper limit)
double CalculateEstimatedError(double target, bool lower=true, double xmin=1, double xmax=0.0)
Return an error estimate on the upper(lower) limit.
static int fgAsymptoticNumPoints
number of points used to build expected p-values
static double fgAsymptoticMaxSigma
max sigma value used to scan asymptotic expected p values
int ExclusionCleanup()
remove points that appear to have failed.
double LowerLimit() override
lower and upper bound of the confidence interval (to get upper/lower limits, multiply the size( = 1-c...
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
TList fYObjects
list of HypoTestResult for each point
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.
virtual double CLsplusb() const
Convert AlternatePValue into a "confidence level".
virtual void Append(const HypoTestResult *other)
add values from another HypoTestResult
void SetBackgroundAsAlt(bool l=true)
void SetNullDistribution(SamplingDistribution *null)
bool GetBackGroundIsAlt(void) const
void SetTestStatisticData(const double tsd)
void SetAltDistribution(SamplingDistribution *alt)
SamplingDistribution * GetNullDistribution(void) const
virtual double CLs() const
is simply (not a method, but a quantity)
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 > & GetSamplingDistribution() const
Get test statistics values.
SimpleInterval is a concrete implementation of the ConfInterval interface.
double fUpperLimit
upper interval limit
RooArgSet fParameters
set containing the parameter of interest
double ConfidenceLevel() const override
return the confidence interval
double fLowerLimit
lower interval limit
double fConfidenceLevel
confidence level
SimpleInterval & operator=(const SimpleInterval &other)
default constructor
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:233
virtual Double_t Derivative(Double_t x, Double_t *params=nullptr, 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:1115
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set lower and upper limits for parameter ipar.
Definition TF1.cxx:3507
virtual void SetParameters(const Double_t *params)
Definition TF1.h:669
virtual void FixParameter(Int_t ipar, Double_t value)
Fix the value of a parameter for a fit operation The specified value will be used in the fit and the ...
Definition TF1.cxx:1559
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition TFile.cxx:4070
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:2315
void Add(TObject *obj) override
Definition TList.h:81
TObject * First() const override
Return the first object in the list. Returns 0 when list is empty.
Definition TList.cxx:659
void Delete(Option_t *option="") override
Remove all objects from the list AND delete all heap based objects.
Definition TList.cxx:470
TObject * At(Int_t idx) const override
Returns the object at position idx. Returns 0 if idx is out of range.
Definition TList.cxx:357
TObject * Clone(const char *newname="") const override
Make a clone of an object using the Streamer facility.
Definition TNamed.cxx:74
const char * GetName() const override
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:223
virtual TObject * RemoveAt(Int_t idx)
Basic string class.
Definition TString.h:139
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2356
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
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 JSONIO.h:26
@ InputArguments
Namespace for the RooStats classes.
Definition Asimov.h:19
TMath.
Definition TMathBase.h:35
Bool_t IsNaN(Double_t x)
Definition TMath.h:892
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
Sort the n1 elements of the Short_t array defined by its iterators.
Definition TMathBase.h:406
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754.
Definition TMath.h:902
Double_t Floor(Double_t x)
Rounds x downward, returning the largest integral value that is not greater than x.
Definition TMath.h:680
T MinElement(Long64_t n, const T *a)
Returns minimum of array a of length n.
Definition TMath.h:960
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=nullptr, Int_t type=7)
Computes sample quantiles, corresponding to the given probabilities.
Definition TMath.cxx:1207
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Comparing floating points.
Definition TMath.h:426
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Comparing floating points.
Definition TMath.h:418
T MaxElement(Long64_t n, const T *a)
Returns maximum of array a of length n.
Definition TMath.h:968
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Binary search in an array of n values to locate value.
Definition TMathBase.h:347
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition TMath.h:917
Definition graph.py:1
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345