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