Logo ROOT  
Reference Guide
HypoTestInverterResult.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
3/*************************************************************************
4 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11
12/** \class RooStats::HypoTestInverterResult
13 \ingroup Roostats
14
15HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence interval.
16Based on the RatioFinder code available in the RooStatsCms package developed by Gregory Schott and Danilo Piparo
17Ported and adapted to RooStats by Gregory Schott
18Some contributions to this class have been written by Matthias Wolf (error estimation)
19
20*/
21
22// include header file of this class
24
28#include "RooMsgService.h"
29#include "RooGlobalFunc.h"
30
31#include "TF1.h"
32#include "TGraphErrors.h"
33#include <cmath>
36#include "Math/Functor.h"
37
38#include "TCanvas.h"
39#include "TFile.h"
41
42#include <algorithm>
43
45
46using namespace RooStats;
47using namespace RooFit;
48using namespace std;
49
50
51// initialize static value
52double HypoTestInverterResult::fgAsymptoticMaxSigma = 5;
53int HypoTestInverterResult::fgAsymptoticNumPoints = 11;
54
55////////////////////////////////////////////////////////////////////////////////
56/// default constructor
57
58HypoTestInverterResult::HypoTestInverterResult(const char * name ) :
60 fUseCLs(false),
61 fIsTwoSided(false),
62 fInterpolateLowerLimit(true),
63 fInterpolateUpperLimit(true),
64 fFittedLowerLimit(false),
65 fFittedUpperLimit(false),
66 fInterpolOption(kLinear),
67 fLowerLimitError(-1),
68 fUpperLimitError(-1),
69 fCLsCleanupThreshold(0.005)
70{
73
76}
77
78////////////////////////////////////////////////////////////////////////////////
79/// copy constructor
80
82 SimpleInterval(other,name),
83 fUseCLs(other.fUseCLs),
84 fIsTwoSided(other.fIsTwoSided),
85 fInterpolateLowerLimit(other.fInterpolateLowerLimit),
86 fInterpolateUpperLimit(other.fInterpolateUpperLimit),
87 fFittedLowerLimit(other.fFittedLowerLimit),
88 fFittedUpperLimit(other.fFittedUpperLimit),
89 fInterpolOption(other.fInterpolOption),
90 fLowerLimitError(other.fLowerLimitError),
91 fUpperLimitError(other.fUpperLimitError),
92 fCLsCleanupThreshold(other.fCLsCleanupThreshold)
93{
96 int nOther = other.ArraySize();
97
98 fXValues = other.fXValues;
99 for (int i = 0; i < nOther; ++i)
100 fYObjects.Add( other.fYObjects.At(i)->Clone() );
101 for (int i = 0; i < fExpPValues.GetSize() ; ++i)
102 fExpPValues.Add( other.fExpPValues.At(i)->Clone() );
103
106}
107
108////////////////////////////////////////////////////////////////////////////////
109
112{
113 if (&other==this) {
114 return *this ;
115 }
116
118 fLowerLimit = other.fLowerLimit;
119 fUpperLimit = other.fUpperLimit;
120 fUseCLs = other.fUseCLs;
121 fIsTwoSided = other.fIsTwoSided;
130
131 int nOther = other.ArraySize();
132 fXValues = other.fXValues;
133
135 for (int i=0; i < nOther; ++i) {
136 fYObjects.Add( other.fYObjects.At(i)->Clone() );
137 }
139 for (int i=0; i < fExpPValues.GetSize() ; ++i) {
140 fExpPValues.Add( other.fExpPValues.At(i)->Clone() );
141 }
142
145
146 return *this;
147}
148
149////////////////////////////////////////////////////////////////////////////////
150/// constructor
151
153 const RooRealVar& scannedVariable,
154 double cl ) :
155 SimpleInterval(name,scannedVariable,TMath::QuietNaN(),TMath::QuietNaN(),cl),
156 fUseCLs(false),
157 fIsTwoSided(false),
158 fInterpolateLowerLimit(true),
159 fInterpolateUpperLimit(true),
160 fFittedLowerLimit(false),
161 fFittedUpperLimit(false),
162 fInterpolOption(kLinear),
163 fLowerLimitError(-1),
164 fUpperLimitError(-1),
165 fCLsCleanupThreshold(0.005)
166{
169
170 // put a cloned copy of scanned variable to set in the interval
171 // to avoid I/O problem of the Result class -
172 // make the set owning the cloned copy (use clone instead of Clone to not copying all links)
175 fParameters.addOwned(*((RooRealVar *) scannedVariable.clone(scannedVariable.GetName()) ));
176}
177
178////////////////////////////////////////////////////////////////////////////////
179/// destructor
180
182{
183 // explicitly empty the TLists - these contain pointers, not objects
186
189}
190
191////////////////////////////////////////////////////////////////////////////////
192/// 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 coutE(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 coutI(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 coutI(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 coutW(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 coutI(Eval) << "HypoTestInverterResult::Add - new number of points is " << fXValues.size()
394 << std::endl;
395 else
396 coutI(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
433{
434 if ( index >= ArraySize() || index<0 ) {
435 coutE(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
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
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
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
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
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
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
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 coutE(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 {
680 std::unique_ptr<TFile> file{TFile::Open(fname,"RECREATE")};
681 graph.Write("graph");
682 }
683#endif
684
685 // look in case if a new intersection exists
686 // only when boundaries are not given
687 if (axmin >= axmax) {
688 int index = TMath::BinarySearch(n, graph.GetX(), limit);
689#ifdef DO_DEBUG
690 std::cout << "do new interpolation dividing from " << index << " and " << y[index] << std::endl;
691#endif
692
693 if (lowSearch && index >= 1 && (y[0] - y0) * ( y[index]- y0) < 0) {
694 //search if another intersection exists at a lower value
695 limit = GetGraphX(graph, y0, lowSearch, graph.GetX()[0], graph.GetX()[index] );
696 }
697 else if (!lowSearch && index < n-2 && (y[n-1] - y0) * ( y[index+1]- y0) < 0) {
698 // another intersection exists at an higher value
699 limit = GetGraphX(graph, y0, lowSearch, graph.GetX()[index+1], graph.GetX()[n-1] );
700 }
701 }
702 // return also xmin, xmax values
703 axmin = xmin;
704 axmax = xmax;
705
706 return limit;
707}
708
709////////////////////////////////////////////////////////////////////////////////
710/// interpolate to find a limit value
711/// Use a linear or a spline interpolation depending on the interpolation option
712
713double HypoTestInverterResult::FindInterpolatedLimit(double target, bool lowSearch, double xmin, double xmax )
714{
715
716 // variable minimum and maximum
717 double varmin = - TMath::Infinity();
718 double varmax = TMath::Infinity();
719 const RooRealVar* var = dynamic_cast<RooRealVar*>( fParameters.first() );
720 if (var) {
721 varmin = var->getMin();
722 varmax = var->getMax();
723 }
724
725 if (ArraySize()<2) {
726 double val = (lowSearch) ? xmin : xmax;
727 coutW(Eval) << "HypoTestInverterResult::FindInterpolatedLimit"
728 << " - not enough points to get the inverted interval - return "
729 << val << std::endl;
730 fLowerLimit = varmin;
731 fUpperLimit = varmax;
732 return (lowSearch) ? fLowerLimit : fUpperLimit;
733 }
734
735 // sort the values in x
736 int n = ArraySize();
737 std::vector<unsigned int> index(n );
738 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
739 // make a graph with the sorted point
740 TGraph graph(n);
741 for (int i = 0; i < n; ++i)
742 graph.SetPoint(i, GetXValue(index[i]), GetYValue(index[i] ) );
743
744
745 //std::cout << " search for " << lowSearch << " xmin = " << xmin << " xmax " << xmax << std::endl;
746
747
748 // search first for min/max in the given range
749 if (xmin >= xmax) {
750
751
752 // search for maximum between the point
753 double * itrmax = std::max_element(graph.GetY() , graph.GetY() +n);
754 double ymax = *itrmax;
755 int iymax = itrmax - graph.GetY();
756 double xwithymax = graph.GetX()[iymax];
757
758#ifdef DO_DEBUG
759 std::cout << " max of y " << iymax << " " << xwithymax << " " << ymax << " target is " << target << std::endl;
760#endif
761 // look if maximum is above/below target
762 if (ymax > target) {
763 if (lowSearch) {
764 if ( iymax > 0) {
765 // low search (minimum is first point or minimum range)
766 xmin = ( graph.GetY()[0] <= target ) ? graph.GetX()[0] : varmin;
767 xmax = xwithymax;
768 }
769 else {
770 // no room for lower limit
771 fLowerLimit = varmin;
772 return fLowerLimit;
773 }
774 }
775 if (!lowSearch ) {
776 // up search
777 if ( iymax < n-1 ) {
778 xmin = xwithymax;
779 xmax = ( graph.GetY()[n-1] <= target ) ? graph.GetX()[n-1] : varmax;
780 }
781 else {
782 // no room for upper limit
783 fUpperLimit = varmax;
784 return fUpperLimit;
785 }
786 }
787 }
788 else {
789 // in case is below the target
790 // find out if is a lower or upper search
791 if (iymax <= (n-1)/2 ) {
792 lowSearch = false;
793 fLowerLimit = varmin;
794 }
795 else {
796 lowSearch = true;
797 fUpperLimit = varmax;
798 }
799 }
800
801#ifdef DO_DEBUG
802 std::cout << " found xmin, xmax = " << xmin << " " << xmax << " for search " << lowSearch << std::endl;
803#endif
804
805 // now come here if I have already found a lower/upper limit
806 // i.e. I am calling routine for the second time
807#ifdef ISNEEDED
808 // should not really come here
809 if (lowSearch && fUpperLimit < varmax) {
810 xmin = fXValues[ index.front() ];
811 // find xmax (is first point before upper limit)
813 if (upI < 1) return xmin;
814 xmax = GetXValue(upI);
815 }
816 else if (!lowSearch && fLowerLimit > varmin ) {
817 // find xmin (is first point after lower limit)
819 if (lowI >= n-1) return xmax;
820 xmin = GetXValue(lowI);
821 xmax = fXValues[ index.back() ];
822 }
823#endif
824 }
825
826#ifdef DO_DEBUG
827 std::cout << "finding " << lowSearch << " limit between " << xmin << " " << xmax << endl;
828#endif
829
830 // compute noe the limit using the TGraph interpolations routine
831 double limit = GetGraphX(graph, target, lowSearch, xmin, xmax);
832 if (lowSearch) fLowerLimit = limit;
833 else fUpperLimit = limit;
834 // estimate the error
835 double error = CalculateEstimatedError( target, lowSearch, xmin, xmax);
836
837 TString limitType = (lowSearch) ? "lower" : "upper";
838 ooccoutD(this,Eval) << "HypoTestInverterResult::FindInterpolateLimit "
839 << "the computed " << limitType << " limit is " << limit << " +/- " << error << std::endl;
840
841#ifdef DO_DEBUG
842 std::cout << "Found limit is " << limit << " +/- " << error << std::endl;
843#endif
844
845
846 if (lowSearch) return fLowerLimit;
847 return fUpperLimit;
848
849
850// if (lowSearch && !TMath::IsNaN(fUpperLimit)) return fLowerLimit;
851// if (!lowSearch && !TMath::IsNaN(fLowerLimit)) return fUpperLimit;
852// // is this needed ?
853// // we call again the function for the upper limits
854
855// // now perform the opposite search on the complement interval
856// if (lowSearch) {
857// xmin = xmax;
858// xmax = varmax;
859// } else {
860// xmax = xmin;
861// xmin = varmin;
862// }
863// double limit2 = GetGraphX(graph, target, !lowSearch, xmin, xmax);
864// if (!lowSearch) fLowerLimit = limit2;
865// else fUpperLimit = limit2;
866
867// CalculateEstimatedError( target, !lowSearch, xmin, xmax);
868
869// #ifdef DO_DEBUG
870// std::cout << "other limit is " << limit2 << std::endl;
871// #endif
872
873// return (lowSearch) ? fLowerLimit : fUpperLimit;
874
875}
876
877////////////////////////////////////////////////////////////////////////////////
878/// - if mode = 0
879/// find closest point to target in Y, the object closest to the target which is 3 sigma from the target
880/// and has smaller error
881/// - if mode = 1
882/// find 2 closest point to target in X and between these two take the one closer to the target
883/// - if mode = 2 as in mode = 1 but return the lower point not the closest one
884/// - if mode = 3 as in mode = 1 but return the upper point not the closest one
885
887{
888
889 int bestIndex = -1;
890 int closestIndex = -1;
891 if (mode == 0) {
892 double smallestError = 2; // error must be < 1
893 double bestValue = 2;
894 for (int i=0; i<ArraySize(); i++) {
895 double dist = fabs(GetYValue(i)-target);
896 if ( dist <3 *GetYError(i) ) { // less than 1 sigma from target CL
897 if (GetYError(i) < smallestError ) {
898 smallestError = GetYError(i);
899 bestIndex = i;
900 }
901 }
902 if ( dist < bestValue) {
903 bestValue = dist;
904 closestIndex = i;
905 }
906 }
907 if (bestIndex >=0) return bestIndex;
908 // if no points found just return the closest one to the target
909 return closestIndex;
910 }
911 // else mode = 1,2,3
912 // find the two closest points to limit value
913 // sort the array first
914 int n = fXValues.size();
915 std::vector<unsigned int> indx(n);
916 TMath::SortItr(fXValues.begin(), fXValues.end(), indx.begin(), false);
917 std::vector<double> xsorted( n);
918 for (int i = 0; i < n; ++i) xsorted[i] = fXValues[indx[i] ];
919 int index1 = TMath::BinarySearch( n, &xsorted[0], xtarget);
920
921#ifdef DO_DEBUG
922 std::cout << "finding closest point to " << xtarget << " is " << index1 << " " << indx[index1] << std::endl;
923#endif
924
925 // case xtarget is outside the range (before or afterwards)
926 if (index1 < 0) return indx[0];
927 if (index1 >= n-1) return indx[n-1];
928 int index2 = index1 +1;
929
930 if (mode == 2) return (GetXValue(indx[index1]) < GetXValue(indx[index2])) ? indx[index1] : indx[index2];
931 if (mode == 3) return (GetXValue(indx[index1]) > GetXValue(indx[index2])) ? indx[index1] : indx[index2];
932 // get smaller point of the two (mode == 1)
933 if (fabs(GetYValue(indx[index1])-target) <= fabs(GetYValue(indx[index2])-target) )
934 return indx[index1];
935 return indx[index2];
936
937}
938
939////////////////////////////////////////////////////////////////////////////////
940
942{
943 if (fFittedLowerLimit) return fLowerLimit;
944 //std::cout << "finding point with cl = " << 1-(1-ConfidenceLevel())/2 << endl;
946 // find both lower/upper limit
948 } else {
949 //LM: I think this is never called
951 }
952 return fLowerLimit;
953}
954
955////////////////////////////////////////////////////////////////////////////////
956
958{
959 //std::cout << "finding point with cl = " << (1-ConfidenceLevel())/2 << endl;
960 if (fFittedUpperLimit) return fUpperLimit;
963 } else {
964 //LM: I think this is never called
966 }
967 return fUpperLimit;
968}
969
970////////////////////////////////////////////////////////////////////////////////
971/// Return an error estimate on the upper(lower) limit. This is the error on
972/// either CLs or CLsplusb divided by an estimate of the slope at this
973/// point.
974
975double HypoTestInverterResult::CalculateEstimatedError(double target, bool lower, double xmin, double xmax)
976{
977
978 if (ArraySize()==0) {
979 coutW(Eval) << "HypoTestInverterResult::CalculateEstimateError"
980 << "Empty result \n";
981 return 0;
982 }
983
984 if (ArraySize()<2) {
985 coutW(Eval) << "HypoTestInverterResult::CalculateEstimateError"
986 << " only points - return its error\n";
987 return GetYError(0);
988 }
989
990 // it does not make sense in case of asymptotic which do not have point errors
991 if (!GetNullTestStatDist(0) ) return 0;
992
993 TString type = (!lower) ? "upper" : "lower";
994
995#ifdef DO_DEBUG
996 std::cout << "calculate estimate error " << type << " between " << xmin << " and " << xmax << std::endl;
997 std::cout << "computed limit is " << ( (lower) ? fLowerLimit : fUpperLimit ) << std::endl;
998#endif
999
1000 // make a TGraph Errors with the sorted points
1001 std::vector<unsigned int> indx(fXValues.size());
1002 TMath::SortItr(fXValues.begin(), fXValues.end(), indx.begin(), false);
1003 // make a graph with the sorted point
1005 int ip = 0, 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(fabs( 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#ifdef DO_DEBUG
1079 std::cout << "closes point to the limit is " << index << " " << GetXValue(index) << " and has error " << GetYError(index) << std::endl;
1080#endif
1081
1082 return theError;
1083}
1084
1085////////////////////////////////////////////////////////////////////////////////
1086/// need to have compute first lower limit
1087
1089{
1091 if (fLowerLimitError >= 0) return fLowerLimitError;
1092 // try to recompute the error
1093 return CalculateEstimatedError(1-ConfidenceLevel(), true);
1094}
1095
1096////////////////////////////////////////////////////////////////////////////////
1097
1099{
1101 if (fUpperLimitError >= 0) return fUpperLimitError;
1102 // try to recompute the error
1103 return CalculateEstimatedError(1-ConfidenceLevel(), false);
1104}
1105
1106////////////////////////////////////////////////////////////////////////////////
1107/// get the background test statistic distribution
1108
1110
1111 HypoTestResult * firstResult = (HypoTestResult*) fYObjects.At(index);
1112 if (!firstResult) return 0;
1113 return firstResult->GetBackGroundIsAlt() ? firstResult->GetAltDistribution() : firstResult->GetNullDistribution();
1114}
1115
1116////////////////////////////////////////////////////////////////////////////////
1117/// get the signal and background test statistic distribution
1118
1121 if (!result) return 0;
1122 return !result->GetBackGroundIsAlt() ? result->GetAltDistribution() : result->GetNullDistribution();
1123}
1124
1125////////////////////////////////////////////////////////////////////////////////
1126/// get the expected p-value distribution at the scanned point index
1127
1129
1130 if (index < 0 || index >= ArraySize() ) return 0;
1131
1132 if (fExpPValues.GetSize() == ArraySize() ) {
1134 }
1135
1136 static bool useFirstB = false;
1137 // get the S+B distribution
1138 int bIndex = (useFirstB) ? 0 : index;
1139
1140 SamplingDistribution * bDistribution = GetBackgroundTestStatDist(bIndex);
1142
1144
1145 if (bDistribution && sbDistribution) {
1146
1147 // create a new HypoTestResult
1148 HypoTestResult tempResult;
1149 tempResult.SetPValueIsRightTail( result->GetPValueIsRightTail() );
1150 tempResult.SetBackgroundAsAlt( true);
1151 // ownership of SamplingDistribution is in HypoTestResult class now
1152 tempResult.SetNullDistribution( new SamplingDistribution(*sbDistribution) );
1153 tempResult.SetAltDistribution( new SamplingDistribution(*bDistribution ) );
1154
1155 std::vector<double> values(bDistribution->GetSize());
1156 for (int i = 0; i < bDistribution->GetSize(); ++i) {
1157 tempResult.SetTestStatisticData( bDistribution->GetSamplingDistribution()[i] );
1158 values[i] = (fUseCLs) ? tempResult.CLs() : tempResult.CLsplusb();
1159 }
1160 return new SamplingDistribution("expected values","expected values",values);
1161 }
1162 // in case b abs sbDistribution are null assume is coming from the asymptotic calculator
1163 // hard -coded this value (no really needed to be used by user)
1166 const double smax = fgAsymptoticMaxSigma;
1167 const int npoints = fgAsymptoticNumPoints;
1168 const double dsig = 2* smax/ (npoints-1) ;
1169 std::vector<double> values(npoints);
1170 for (int i = 0; i < npoints; ++i) {
1171 double nsig = -smax + dsig*i;
1172 double pval = AsymptoticCalculator::GetExpectedPValues( result->NullPValue(), result->AlternatePValue(), nsig, fUseCLs, !fIsTwoSided);
1173 if (pval < 0) { return 0;}
1174 values[i] = pval;
1175 }
1176 return new SamplingDistribution("Asymptotic expected values","Asymptotic expected values",values);
1177
1178}
1179
1180////////////////////////////////////////////////////////////////////////////////
1181/// get the limit distribution (lower/upper depending on the flag)
1182/// by interpolating the expected p values for each point
1183
1185 if (ArraySize()<2) {
1186 coutE(Eval) << "HypoTestInverterResult::GetLimitDistribution"
1187 << " not enough points - return 0 " << std::endl;
1188 return 0;
1189 }
1190
1191 ooccoutD(this,Eval) << "HypoTestInverterResult - computing limit distribution...." << std::endl;
1192
1193
1194 // find optimal size by looking at the PValue distribution obtained
1195 int npoints = ArraySize();
1196 std::vector<SamplingDistribution*> distVec( npoints );
1197 double sum = 0;
1198 for (unsigned int i = 0; i < distVec.size(); ++i) {
1199 distVec[i] = GetExpectedPValueDist(i);
1200 // sort the distributions
1201 // hack (by calling InverseCDF(0) will sort the sampling distribution
1202 if (distVec[i] ) {
1203 distVec[i]->InverseCDF(0);
1204 sum += distVec[i]->GetSize();
1205 }
1206 }
1207 int size = int( sum/ npoints);
1208
1209 if (size < 10) {
1210 ooccoutW(this,InputArguments) << "HypoTestInverterResult - set a minimum size of 10 for limit distribution" << std::endl;
1211 size = 10;
1212 }
1213
1214
1215 double target = 1-fConfidenceLevel;
1216
1217 // vector with the quantiles of the p-values for each scanned poi point
1218 std::vector< std::vector<double> > quantVec(npoints );
1219 for (int i = 0; i < npoints; ++i) {
1220
1221 if (!distVec[i]) continue;
1222
1223 // make quantiles from the sampling distributions of the expected p values
1224 std::vector<double> pvalues = distVec[i]->GetSamplingDistribution();
1225 delete distVec[i]; distVec[i] = 0;
1226 std::sort(pvalues.begin(), pvalues.end());
1227 // find the quantiles of the distribution
1228 double p[1] = {0};
1229 double q[1] = {0};
1230
1231 quantVec[i] = std::vector<double>(size);
1232 for (int ibin = 0; ibin < size; ++ibin) {
1233 // exclude for a bug in TMath::Quantiles last item
1234 p[0] = std::min( (ibin+1) * 1./double(size), 1.0);
1235 // use the type 1 which give the point value
1236 TMath::Quantiles(pvalues.size(), 1, &pvalues[0], q, p, true, 0, 1 );
1237 (quantVec[i])[ibin] = q[0];
1238 }
1239
1240 }
1241
1242 // sort the values in x
1243 std::vector<unsigned int> index( npoints );
1244 TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1245
1246 // SamplingDistribution * dist0 = distVec[index[0]];
1247 // SamplingDistribution * dist1 = distVec[index[1]];
1248
1249
1250 std::vector<double> limits(size);
1251 // loop on the p values and find the limit for each expected point in the quantiles vector
1252 for (int j = 0; j < size; ++j ) {
1253
1254 TGraph g;
1255 for (int k = 0; k < npoints ; ++k) {
1256 if (quantVec[index[k]].size() > 0 )
1257 g.SetPoint(g.GetN(), GetXValue(index[k]), (quantVec[index[k]])[j] );
1258 }
1259
1260 limits[j] = GetGraphX(g, target, lower);
1261
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 != 0);
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}
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
Definition: RooMsgService.h:34
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
#define ooccoutW(o, a)
Definition: RooMsgService.h:59
#define ooccoutI(o, a)
Definition: RooMsgService.h:57
#define ooccoutD(o, a)
Definition: RooMsgService.h:56
#define ooccoutE(o, a)
Definition: RooMsgService.h:60
#define ClassImp(name)
Definition: Rtypes.h:375
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 g
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
Definition: THbookFile.cxx:95
float * q
Definition: THbookFile.cxx:89
float ymin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
float ymax
Definition: THbookFile.cxx:95
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:506
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.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
TObject * clone(const char *newname) const override
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 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)
void SetPValueIsRightTail(bool pr)
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.
Definition: TCollection.h:184
1-Dim function class
Definition: TF1.h:213
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:1114
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set lower and upper limits for parameter ipar.
Definition: TF1.cxx:3536
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:649
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:1564
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:4053
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
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:136
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:2345
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1778
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
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
@ InputArguments
Definition: RooGlobalFunc.h:62
Namespace for the RooStats classes.
Definition: Asimov.h:19
static constexpr double s
TMath.
Definition: TMathBase.h:35
Bool_t IsNaN(Double_t x)
Definition: TMath.h:890
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:900
Double_t Floor(Double_t x)
Rounds x downward, returning the largest integral value that is not greater than x.
Definition: TMath.h:678
T MinElement(Long64_t n, const T *a)
Returns minimum of array a of length n.
Definition: TMath.h:958
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:425
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:966
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:915
Definition: file.py:1
Definition: graph.py:1
TMarker m
Definition: textangle.C:8