Logo ROOT   6.12/07
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 
59  SimpleInterval(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),
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),
160  fFittedLowerLimit(false),
161  fFittedUpperLimit(false),
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  double xmin = axmin;
614  double xmax = axmax;
615 
616  // case no range is given check if need to extrapolate to lower/upper values
617  if (axmin >= axmax ) {
618 
619 #ifdef DO_DEBUG
620  std::cout << "No rage given - check if extrapolation is needed " << std::endl;
621 #endif
622 
623  xmin = graph.GetX()[0];
624  xmax = graph.GetX()[n-1];
625 
626  double yfirst = graph.GetY()[0];
627  double ylast = graph.GetY()[n-1];
628 
629  // case we have lower/upper limits
630 
631  // find ymin and ymax and corresponding values
632  //double ymin = TMath::MinElement(n,y);
633  double ymax = TMath::MaxElement(n,y);
634  // do lower extrapolation
635 
636  if ( (ymax < y0 && !lowSearch) || ( yfirst > y0 && lowSearch) ) {
637  xmin = varmin;
638  }
639  // do upper extrapolation
640  if ( (ymax < y0 && lowSearch) || ( ylast > y0 && !lowSearch) ) {
641  xmax = varmax;
642  }
643  }
644 
645  auto func = [&](double x) {
646  return (fInterpolOption == kSpline) ? graph.Eval(x, nullptr, "S") - y0 : graph.Eval(x) - y0;
647  };
648  ROOT::Math::Functor1D f1d(func);
649 
651  brf.SetFunction(f1d,xmin,xmax);
652  brf.SetNpx(TMath::Max(graph.GetN()*2,100) );
653  bool ret = brf.Solve(100, 1.E-16, 1.E-6);
654  if (!ret) {
655  ooccoutE(this,Eval) << "HypoTestInverterResult - interpolation failed - return inf" << std::endl;
656  return TMath::Infinity();
657  }
658  double limit = brf.Root();
659 
660  // auto grfunc = [&](double * x, double *) {
661  // return (fInterpolOption == kSpline) ? graph.Eval(*x, nullptr, "S") - y0 : graph.Eval(*x) - y0;
662  // };
663  // TF1 tgrfunc("tgrfunc",grfunc,xmin,xmax,0);
664  // double limit = tgrfunc.GetX(0,xmin,xmax);
665 
666 #ifdef DO_DEBUG
667  if (lowSearch) std::cout << "lower limit search : ";
668  else std::cout << "Upper limit search : ";
669  std::cout << "interpolation done between " << xmin << " and " << xmax
670  << "\n Found limit using RootFinder is " << limit << std::endl;
671 
672  TString fname = "graph_upper.root";
673  if (lowSearch) fname = "graph_lower.root";
674  auto file = TFile::Open(fname,"RECREATE");
675  graph.Write("graph");
676  file->Close();
677 
678 #endif
679 
680  // look in case if a new intersection exists
681  // only when boundaries are not given
682  if (axmin >= axmax) {
683  int index = TMath::BinarySearch(n, graph.GetX(), limit);
684 #ifdef DO_DEBUG
685  std::cout << "do new interpolation dividing from " << index << " and " << y[index] << std::endl;
686 #endif
687 
688  if (lowSearch && index >= 1 && (y[0] - y0) * ( y[index]- y0) < 0) {
689  //search if another intersection exists at a lower value
690  limit = GetGraphX(graph, y0, lowSearch, graph.GetX()[0], graph.GetX()[index] );
691  }
692  else if (!lowSearch && index < n-2 && (y[n-1] - y0) * ( y[index+1]- y0) < 0) {
693  // another intersection exists at an higher value
694  limit = GetGraphX(graph, y0, lowSearch, graph.GetX()[index+1], graph.GetX()[n-1] );
695  }
696  }
697  // return also xmin, xmax values
698  axmin = xmin;
699  axmax = xmax;
700 
701  return limit;
702 }
703 
704 ////////////////////////////////////////////////////////////////////////////////
705 /// interpolate to find a limit value
706 /// Use a linear or a spline interpolation depending on the interpolation option
707 
708 double HypoTestInverterResult::FindInterpolatedLimit(double target, bool lowSearch, double xmin, double xmax )
709 {
710 
711  // variable minimum and maximum
712  double varmin = - TMath::Infinity();
713  double varmax = TMath::Infinity();
714  const RooRealVar* var = dynamic_cast<RooRealVar*>( fParameters.first() );
715  if (var) {
716  varmin = var->getMin();
717  varmax = var->getMax();
718  }
719 
720  if (ArraySize()<2) {
721  double val = (lowSearch) ? xmin : xmax;
722  oocoutW(this,Eval) << "HypoTestInverterResult::FindInterpolatedLimit"
723  << " - not enough points to get the inverted interval - return "
724  << val << std::endl;
725  fLowerLimit = varmin;
726  fUpperLimit = varmax;
727  return (lowSearch) ? fLowerLimit : fUpperLimit;
728  }
729 
730  // sort the values in x
731  int n = ArraySize();
732  std::vector<unsigned int> index(n );
733  TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
734  // make a graph with the sorted point
735  TGraph graph(n);
736  for (int i = 0; i < n; ++i)
737  graph.SetPoint(i, GetXValue(index[i]), GetYValue(index[i] ) );
738 
739 
740  //std::cout << " search for " << lowSearch << " xmin = " << xmin << " xmax " << xmax << std::endl;
741 
742 
743  // search first for min/max in the given range
744  if (xmin >= xmax) {
745 
746 
747  // search for maximum between the point
748  double * itrmax = std::max_element(graph.GetY() , graph.GetY() +n);
749  double ymax = *itrmax;
750  int iymax = itrmax - graph.GetY();
751  double xwithymax = graph.GetX()[iymax];
752 
753 #ifdef DO_DEBUG
754  std::cout << " max of y " << iymax << " " << xwithymax << " " << ymax << " target is " << target << std::endl;
755 #endif
756  // look if maximum is above/below target
757  if (ymax > target) {
758  if (lowSearch) {
759  if ( iymax > 0) {
760  // low search (minimum is first point or minimum range)
761  xmin = ( graph.GetY()[0] <= target ) ? graph.GetX()[0] : varmin;
762  xmax = xwithymax;
763  }
764  else {
765  // no room for lower limit
766  fLowerLimit = varmin;
767  return fLowerLimit;
768  }
769  }
770  if (!lowSearch ) {
771  // up search
772  if ( iymax < n-1 ) {
773  xmin = xwithymax;
774  xmax = ( graph.GetY()[n-1] <= target ) ? graph.GetX()[n-1] : varmax;
775  }
776  else {
777  // no room for upper limit
778  fUpperLimit = varmax;
779  return fUpperLimit;
780  }
781  }
782  }
783  else {
784  // in case is below the target
785  // find out if is a lower or upper search
786  if (iymax <= (n-1)/2 ) {
787  lowSearch = false;
788  fLowerLimit = varmin;
789  }
790  else {
791  lowSearch = true;
792  fUpperLimit = varmax;
793  }
794  }
795 
796 #ifdef DO_DEBUG
797  std::cout << " found xmin, xmax = " << xmin << " " << xmax << " for search " << lowSearch << std::endl;
798 #endif
799 
800  // now come here if I have already found a lower/upper limit
801  // i.e. I am calling routine for the second time
802 #ifdef ISNEEDED
803  // should not really come here
804  if (lowSearch && fUpperLimit < varmax) {
805  xmin = fXValues[ index.front() ];
806  // find xmax (is first point before upper limit)
807  int upI = FindClosestPointIndex(target, 2, fUpperLimit);
808  if (upI < 1) return xmin;
809  xmax = GetXValue(upI);
810  }
811  else if (!lowSearch && fLowerLimit > varmin ) {
812  // find xmin (is first point after lower limit)
813  int lowI = FindClosestPointIndex(target, 3, fLowerLimit);
814  if (lowI >= n-1) return xmax;
815  xmin = GetXValue(lowI);
816  xmax = fXValues[ index.back() ];
817  }
818 #endif
819  }
820 
821 #ifdef DO_DEBUG
822  std::cout << "finding " << lowSearch << " limit between " << xmin << " " << xmax << endl;
823 #endif
824 
825  // compute noe the limit using the TGraph interpolations routine
826  double limit = GetGraphX(graph, target, lowSearch, xmin, xmax);
827  if (lowSearch) fLowerLimit = limit;
828  else fUpperLimit = limit;
829  // estimate the error
830  double error = CalculateEstimatedError( target, lowSearch, xmin, xmax);
831 
832  TString limitType = (lowSearch) ? "lower" : "upper";
833  ooccoutD(this,Eval) << "HypoTestInverterResult::FindInterpolateLimit "
834  << "the computed " << limitType << " limit is " << limit << " +/- " << error << std::endl;
835 
836 #ifdef DO_DEBUG
837  std::cout << "Found limit is " << limit << " +/- " << error << astd::endl;
838 #endif
839 
840 
841  if (lowSearch) return fLowerLimit;
842  return fUpperLimit;
843 
844 
845 // if (lowSearch && !TMath::IsNaN(fUpperLimit)) return fLowerLimit;
846 // if (!lowSearch && !TMath::IsNaN(fLowerLimit)) return fUpperLimit;
847 // // is this needed ?
848 // // we call again the function for the upper limits
849 
850 // // now perform the opposite search on the complement interval
851 // if (lowSearch) {
852 // xmin = xmax;
853 // xmax = varmax;
854 // } else {
855 // xmax = xmin;
856 // xmin = varmin;
857 // }
858 // double limit2 = GetGraphX(graph, target, !lowSearch, xmin, xmax);
859 // if (!lowSearch) fLowerLimit = limit2;
860 // else fUpperLimit = limit2;
861 
862 // CalculateEstimatedError( target, !lowSearch, xmin, xmax);
863 
864 // #ifdef DO_DEBUG
865 // std::cout << "other limit is " << limit2 << std::endl;
866 // #endif
867 
868 // return (lowSearch) ? fLowerLimit : fUpperLimit;
869 
870 }
871 
872 ////////////////////////////////////////////////////////////////////////////////
873 /// - if mode = 0
874 /// find closest point to target in Y, the object closest to the target which is 3 sigma from the target
875 /// and has smaller error
876 /// - if mode = 1
877 /// find 2 closest point to target in X and between these two take the one closer to the target
878 /// - if mode = 2 as in mode = 1 but return the lower point not the closest one
879 /// - if mode = 3 as in mode = 1 but return the upper point not the closest one
880 
881 int HypoTestInverterResult::FindClosestPointIndex(double target, int mode, double xtarget)
882 {
883 
884  int bestIndex = -1;
885  int closestIndex = -1;
886  if (mode == 0) {
887  double smallestError = 2; // error must be < 1
888  double bestValue = 2;
889  for (int i=0; i<ArraySize(); i++) {
890  double dist = fabs(GetYValue(i)-target);
891  if ( dist <3 *GetYError(i) ) { // less than 1 sigma from target CL
892  if (GetYError(i) < smallestError ) {
893  smallestError = GetYError(i);
894  bestIndex = i;
895  }
896  }
897  if ( dist < bestValue) {
898  bestValue = dist;
899  closestIndex = i;
900  }
901  }
902  if (bestIndex >=0) return bestIndex;
903  // if no points found just return the closest one to the target
904  return closestIndex;
905  }
906  // else mode = 1,2,3
907  // find the two closest points to limit value
908  // sort the array first
909  int n = fXValues.size();
910  std::vector<unsigned int> indx(n);
911  TMath::SortItr(fXValues.begin(), fXValues.end(), indx.begin(), false);
912  std::vector<double> xsorted( n);
913  for (int i = 0; i < n; ++i) xsorted[i] = fXValues[indx[i] ];
914  int index1 = TMath::BinarySearch( n, &xsorted[0], xtarget);
915 
916 #ifdef DO_DEBUG
917  std::cout << "finding closest point to " << xtarget << " is " << index1 << " " << indx[index1] << std::endl;
918 #endif
919 
920  // case xtarget is outside the range (before or afterwards)
921  if (index1 < 0) return indx[0];
922  if (index1 >= n-1) return indx[n-1];
923  int index2 = index1 +1;
924 
925  if (mode == 2) return (GetXValue(indx[index1]) < GetXValue(indx[index2])) ? indx[index1] : indx[index2];
926  if (mode == 3) return (GetXValue(indx[index1]) > GetXValue(indx[index2])) ? indx[index1] : indx[index2];
927  // get smaller point of the two (mode == 1)
928  if (fabs(GetYValue(indx[index1])-target) <= fabs(GetYValue(indx[index2])-target) )
929  return indx[index1];
930  return indx[index2];
931 
932 }
933 
934 ////////////////////////////////////////////////////////////////////////////////
935 
937 {
938  if (fFittedLowerLimit) return fLowerLimit;
939  //std::cout << "finding point with cl = " << 1-(1-ConfidenceLevel())/2 << endl;
940  if ( fInterpolateLowerLimit ) {
941  // find both lower/upper limit
943  } else {
944  //LM: I think this is never called
946  }
947  return fLowerLimit;
948 }
949 
950 ////////////////////////////////////////////////////////////////////////////////
951 
953 {
954  //std::cout << "finding point with cl = " << (1-ConfidenceLevel())/2 << endl;
955  if (fFittedUpperLimit) return fUpperLimit;
956  if ( fInterpolateUpperLimit ) {
958  } else {
959  //LM: I think this is never called
961  }
962  return fUpperLimit;
963 }
964 
965 ////////////////////////////////////////////////////////////////////////////////
966 /// Return an error estimate on the upper(lower) limit. This is the error on
967 /// either CLs or CLsplusb divided by an estimate of the slope at this
968 /// point.
969 
970 Double_t HypoTestInverterResult::CalculateEstimatedError(double target, bool lower, double xmin, double xmax)
971 {
972 
973  if (ArraySize()==0) {
974  oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimateError"
975  << "Empty result \n";
976  return 0;
977  }
978 
979  if (ArraySize()<2) {
980  oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimateError"
981  << " only points - return its error\n";
982  return GetYError(0);
983  }
984 
985  // it does not make sense in case of asymptotic which do not have point errors
986  if (!GetNullTestStatDist(0) ) return 0;
987 
988  TString type = (!lower) ? "upper" : "lower";
989 
990 #ifdef DO_DEBUG
991  std::cout << "calculate estimate error " << type << " between " << xmin << " and " << xmax << std::endl;
992  std::cout << "computed limit is " << ( (lower) ? fLowerLimit : fUpperLimit ) << std::endl;
993 #endif
994 
995  // make a TGraph Errors with the sorted points
996  std::vector<unsigned int> indx(fXValues.size());
997  TMath::SortItr(fXValues.begin(), fXValues.end(), indx.begin(), false);
998  // make a graph with the sorted point
1000  int ip = 0, np = 0;
1001  for (int i = 0; i < ArraySize(); ++i) {
1002  if ( (xmin < xmax) && ( GetXValue(indx[i]) >= xmin && GetXValue(indx[i]) <= xmax) ) {
1003  np++;
1004  // exclude points with zero or very small errors
1005  if (GetYError(indx[i] ) > 1.E-6) {
1006  graph.SetPoint(ip, GetXValue(indx[i]), GetYValue(indx[i] ) );
1007  graph.SetPointError(ip, 0., GetYError(indx[i]) );
1008  ip++;
1009  }
1010  }
1011  }
1012  if (graph.GetN() < 2) {
1013  if (np >= 2) oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimatedError - no valid points - cannot estimate the " << type << " limit error " << std::endl;
1014  return 0;
1015  }
1016 
1017  double minX = xmin;
1018  double maxX = xmax;
1019  if (xmin >= xmax) {
1020  minX = fXValues[ indx.front() ];
1021  maxX = fXValues[ indx.back() ];
1022  }
1023 
1024  TF1 fct("fct", "exp([0] * (x - [2] ) + [1] * (x-[2])**2)", minX, maxX);
1025  double scale = maxX-minX;
1026  if (lower) {
1027  fct.SetParameters( 2./scale, 0.1/scale, graph.GetX()[0] );
1028  fct.SetParLimits(0,0,100./scale);
1029  fct.SetParLimits(1,0, 10./scale); }
1030  else {
1031  fct.SetParameters( -2./scale, -0.1/scale );
1032  fct.SetParLimits(0,-100./scale, 0);
1033  fct.SetParLimits(1,-100./scale, 0); }
1034 
1035  if (graph.GetN() < 3) fct.FixParameter(1,0.);
1036 
1037  // find the point closest to the limit
1038  double limit = (!lower) ? fUpperLimit : fLowerLimit;
1039  if (TMath::IsNaN(limit)) return 0; // cannot do if limit not computed
1040 
1041 
1042 #ifdef DO_DEBUG
1043  TCanvas * c1 = new TCanvas();
1044  std::cout << "fitting for limit " << type << "between " << minX << " , " << maxX << " points considered " << graph.GetN() << std::endl;
1045  int fitstat = graph.Fit(&fct," EX0");
1046  graph.SetMarkerStyle(20);
1047  graph.Draw("AP");
1048  graph.Print();
1049  c1->SaveAs(TString::Format("graphFit_%s.pdf",type.Data()) );
1050  delete c1;
1051 #else
1052  int fitstat = graph.Fit(&fct,"Q EX0");
1053 #endif
1054 
1055  int index = FindClosestPointIndex(target, 1, limit);
1056  double theError = 0;
1057  if (fitstat == 0) {
1058  double errY = GetYError(index);
1059  if (errY > 0) {
1060  double m = fct.Derivative( GetXValue(index) );
1061  theError = std::min(fabs( GetYError(index) / m), maxX-minX);
1062  }
1063  }
1064  else {
1065  oocoutW(this,Eval) << "HypoTestInverterResult::CalculateEstimatedError - cannot estimate the " << type << " limit error " << std::endl;
1066  theError = 0;
1067  }
1068  if (lower)
1069  fLowerLimitError = theError;
1070  else
1071  fUpperLimitError = theError;
1072 
1073 #ifdef DO_DEBUG
1074  std::cout << "closes point to the limit is " << index << " " << GetXValue(index) << " and has error " << GetYError(index) << std::endl;
1075 #endif
1076 
1077  return theError;
1078 }
1079 
1080 ////////////////////////////////////////////////////////////////////////////////
1081 /// need to have compute first lower limit
1082 
1084 {
1086  if (fLowerLimitError >= 0) return fLowerLimitError;
1087  // try to recompute the error
1088  return CalculateEstimatedError(1-ConfidenceLevel(), true);
1089 }
1090 
1091 ////////////////////////////////////////////////////////////////////////////////
1092 
1094 {
1096  if (fUpperLimitError >= 0) return fUpperLimitError;
1097  // try to recompute the error
1098  return CalculateEstimatedError(1-ConfidenceLevel(), false);
1099 }
1100 
1101 ////////////////////////////////////////////////////////////////////////////////
1102 /// get the background test statistic distribution
1103 
1105 
1106  HypoTestResult * firstResult = (HypoTestResult*) fYObjects.At(index);
1107  if (!firstResult) return 0;
1108  return firstResult->GetBackGroundIsAlt() ? firstResult->GetAltDistribution() : firstResult->GetNullDistribution();
1109 }
1110 
1111 ////////////////////////////////////////////////////////////////////////////////
1112 /// get the signal and background test statistic distribution
1113 
1115  HypoTestResult * result = (HypoTestResult*) fYObjects.At(index);
1116  if (!result) return 0;
1117  return !result->GetBackGroundIsAlt() ? result->GetAltDistribution() : result->GetNullDistribution();
1118 }
1119 
1120 ////////////////////////////////////////////////////////////////////////////////
1121 /// get the expected p-value distribution at the scanned point index
1122 
1124 
1125  if (index < 0 || index >= ArraySize() ) return 0;
1126 
1127  if (fExpPValues.GetSize() == ArraySize() ) {
1128  return (SamplingDistribution*) fExpPValues.At(index)->Clone();
1129  }
1130 
1131  static bool useFirstB = false;
1132  // get the S+B distribution
1133  int bIndex = (useFirstB) ? 0 : index;
1134 
1135  SamplingDistribution * bDistribution = GetBackgroundTestStatDist(bIndex);
1136  SamplingDistribution * sbDistribution = GetSignalAndBackgroundTestStatDist(index);
1137 
1138  HypoTestResult * result = (HypoTestResult*) fYObjects.At(index);
1139 
1140  if (bDistribution && sbDistribution) {
1141 
1142  // create a new HypoTestResult
1143  HypoTestResult tempResult;
1144  tempResult.SetPValueIsRightTail( result->GetPValueIsRightTail() );
1145  tempResult.SetBackgroundAsAlt( true);
1146  // ownership of SamplingDistribution is in HypoTestResult class now
1147  tempResult.SetNullDistribution( new SamplingDistribution(*sbDistribution) );
1148  tempResult.SetAltDistribution( new SamplingDistribution(*bDistribution ) );
1149 
1150  std::vector<double> values(bDistribution->GetSize());
1151  for (int i = 0; i < bDistribution->GetSize(); ++i) {
1152  tempResult.SetTestStatisticData( bDistribution->GetSamplingDistribution()[i] );
1153  values[i] = (fUseCLs) ? tempResult.CLs() : tempResult.CLsplusb();
1154  }
1155  return new SamplingDistribution("expected values","expected values",values);
1156  }
1157  // in case b abs sbDistribution are null assume is coming from the asymptotic calculator
1158  // hard -coded this value (no really needed to be used by user)
1161  const double smax = fgAsymptoticMaxSigma;
1162  const int npoints = fgAsymptoticNumPoints;
1163  const double dsig = 2* smax/ (npoints-1) ;
1164  std::vector<double> values(npoints);
1165  for (int i = 0; i < npoints; ++i) {
1166  double nsig = -smax + dsig*i;
1167  double pval = AsymptoticCalculator::GetExpectedPValues( result->NullPValue(), result->AlternatePValue(), nsig, fUseCLs, !fIsTwoSided);
1168  if (pval < 0) { return 0;}
1169  values[i] = pval;
1170  }
1171  return new SamplingDistribution("Asymptotic expected values","Asymptotic expected values",values);
1172 
1173 }
1174 
1175 ////////////////////////////////////////////////////////////////////////////////
1176 /// get the limit distribution (lower/upper depending on the flag)
1177 /// by interpolating the expected p values for each point
1178 
1180  if (ArraySize()<2) {
1181  oocoutE(this,Eval) << "HypoTestInverterResult::GetLimitDistribution"
1182  << " not enough points - return 0 " << std::endl;
1183  return 0;
1184  }
1185 
1186  ooccoutD(this,Eval) << "HypoTestInverterResult - computing limit distribution...." << std::endl;
1187 
1188 
1189  // find optimal size by looking at the PValue distribution obtained
1190  int npoints = ArraySize();
1191  std::vector<SamplingDistribution*> distVec( npoints );
1192  double sum = 0;
1193  for (unsigned int i = 0; i < distVec.size(); ++i) {
1194  distVec[i] = GetExpectedPValueDist(i);
1195  // sort the distributions
1196  // hack (by calling InverseCDF(0) will sort the sampling distribution
1197  if (distVec[i] ) {
1198  distVec[i]->InverseCDF(0);
1199  sum += distVec[i]->GetSize();
1200  }
1201  }
1202  int size = int( sum/ npoints);
1203 
1204  if (size < 10) {
1205  ooccoutW(this,InputArguments) << "HypoTestInverterResult - set a minimum size of 10 for limit distribution" << std::endl;
1206  size = 10;
1207  }
1208 
1209 
1210  double target = 1-fConfidenceLevel;
1211 
1212  // vector with the quantiles of the p-values for each scanned poi point
1213  std::vector< std::vector<double> > quantVec(npoints );
1214  for (int i = 0; i < npoints; ++i) {
1215 
1216  if (!distVec[i]) continue;
1217 
1218  // make quantiles from the sampling distributions of the expected p values
1219  std::vector<double> pvalues = distVec[i]->GetSamplingDistribution();
1220  delete distVec[i]; distVec[i] = 0;
1221  std::sort(pvalues.begin(), pvalues.end());
1222  // find the quantiles of the distribution
1223  double p[1] = {0};
1224  double q[1] = {0};
1225 
1226  quantVec[i] = std::vector<double>(size);
1227  for (int ibin = 0; ibin < size; ++ibin) {
1228  // exclude for a bug in TMath::Quantiles last item
1229  p[0] = std::min( (ibin+1) * 1./double(size), 1.0);
1230  // use the type 1 which give the point value
1231  TMath::Quantiles(pvalues.size(), 1, &pvalues[0], q, p, true, 0, 1 );
1232  (quantVec[i])[ibin] = q[0];
1233  }
1234 
1235  }
1236 
1237  // sort the values in x
1238  std::vector<unsigned int> index( npoints );
1239  TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1240 
1241  // SamplingDistribution * dist0 = distVec[index[0]];
1242  // SamplingDistribution * dist1 = distVec[index[1]];
1243 
1244 
1245  std::vector<double> limits(size);
1246  // loop on the p values and find the limit for each expected point in the quantiles vector
1247  for (int j = 0; j < size; ++j ) {
1248 
1249  TGraph g;
1250  for (int k = 0; k < npoints ; ++k) {
1251  if (quantVec[index[k]].size() > 0 )
1252  g.SetPoint(g.GetN(), GetXValue(index[k]), (quantVec[index[k]])[j] );
1253  }
1254 
1255  limits[j] = GetGraphX(g, target, lower);
1256 
1257  }
1258 
1259 
1260  if (lower)
1261  return new SamplingDistribution("Expected lower Limit","Expected lower limits",limits);
1262  else
1263  return new SamplingDistribution("Expected upper Limit","Expected upper limits",limits);
1264 
1265 }
1266 
1267 ////////////////////////////////////////////////////////////////////////////////
1268 /// Get the expected lower limit
1269 /// nsig is used to specify which expected value of the UpperLimitDistribution
1270 /// For example
1271 /// - nsig = 0 (default value) returns the expected value
1272 /// - nsig = -1 returns the lower band value at -1 sigma
1273 /// - nsig + 1 return the upper value
1274 /// - opt = "" (default) : compute limit by interpolating all the p values, find the corresponding limit distribution
1275 /// and then find the quantiles in the limit distribution
1276 /// ioption = "P" is the method used for plotting. One Finds the corresponding nsig quantile in the p values and then
1277 /// interpolates them
1278 
1279 double HypoTestInverterResult::GetExpectedLowerLimit(double nsig, const char * opt ) const {
1280 
1281  return GetExpectedLimit(nsig, true, opt);
1282 }
1283 
1284 ////////////////////////////////////////////////////////////////////////////////
1285 /// Get the expected upper limit
1286 /// nsig is used to specify which expected value of the UpperLimitDistribution
1287 /// For example
1288 /// - nsig = 0 (default value) returns the expected value
1289 /// - nsig = -1 returns the lower band value at -1 sigma
1290 /// - nsig + 1 return the upper value
1291 /// - opt is an option specifying the type of method used for computing the upper limit
1292 /// - opt = "" (default) : compute limit by interpolating all the p values, find the corresponding limit distribution
1293 /// and then find the quantiles in the limit distribution
1294 /// ioption = "P" is the method used for plotting. One Finds the corresponding nsig quantile in the p values and then
1295 /// interpolates them
1296 
1297 double HypoTestInverterResult::GetExpectedUpperLimit(double nsig, const char * opt ) const {
1298 
1299  return GetExpectedLimit(nsig, false, opt);
1300 }
1301 
1302 ////////////////////////////////////////////////////////////////////////////////
1303 /// get expected limit (lower/upper) depending on the flag
1304 /// for asymptotic is a special case (the distribution is generated an step in sigma values)
1305 /// distinguish asymptotic looking at the hypotest results
1306 /// if option = "P" get expected limit using directly quantiles of p value distribution
1307 /// else (default) find expected limit by obtaining first a full limit distributions
1308 /// The last one is in general more correct
1309 
1310 double HypoTestInverterResult::GetExpectedLimit(double nsig, bool lower, const char * opt ) const {
1311 
1312  const int nEntries = ArraySize();
1313  if (nEntries <= 0) return (lower) ? 1 : 0; // return 1 for lower, 0 for upper
1314 
1315  HypoTestResult * r = dynamic_cast<HypoTestResult *> (fYObjects.First() );
1316  assert(r != 0);
1317  if (!r->GetNullDistribution() && !r->GetAltDistribution() ) {
1318  // we are in the asymptotic case
1319  // get the limits obtained at the different sigma values
1320  SamplingDistribution * limitDist = GetLimitDistribution(lower);
1321  if (!limitDist) return 0;
1322  const std::vector<double> & values = limitDist->GetSamplingDistribution();
1323  if (values.size() <= 1) return 0;
1324  double dsig = 2* fgAsymptoticMaxSigma/ (values.size() -1) ;
1325  int i = (int) TMath::Floor ( (nsig + fgAsymptoticMaxSigma)/dsig + 0.5);
1326  return values[i];
1327  }
1328 
1329  double p[1] = {0};
1330  double q[1] = {0};
1331  p[0] = ROOT::Math::normal_cdf(nsig,1);
1332 
1333  // for CLs+b can get the quantiles of p-value distribution and
1334  // interpolate them
1335  // In the case of CLs (since it is not a real p-value anymore but a ratio)
1336  // then it is needed to get first all limit distribution values and then the quantiles
1337 
1338  TString option(opt);
1339  option.ToUpper();
1340  if (option.Contains("P")) {
1341 
1342  TGraph g;
1343 
1344  // sort the arrays based on the x values
1345  std::vector<unsigned int> index(nEntries);
1346  TMath::SortItr(fXValues.begin(), fXValues.end(), index.begin(), false);
1347 
1348  for (int j=0; j<nEntries; ++j) {
1349  int i = index[j]; // i is the order index
1351  if (!s) {
1352  ooccoutI(this,Eval) << "HypoTestInverterResult - cannot compute expected p value distribution for point, x = "
1353  << GetXValue(i) << " skip it " << std::endl;
1354  continue;
1355  }
1356  const std::vector<double> & values = s->GetSamplingDistribution();
1357  double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
1358  TMath::Quantiles(values.size(), 1, x,q,p,false);
1359  g.SetPoint(g.GetN(), fXValues[i], q[0] );
1360  delete s;
1361  }
1362  if (g.GetN() < 2) {
1363  ooccoutE(this,Eval) << "HypoTestInverterResult - cannot compute limits , not enough points, n = " << g.GetN() << std::endl;
1364  return 0;
1365  }
1366 
1367  // interpolate the graph to obtain the limit
1368  double target = 1-fConfidenceLevel;
1369  return GetGraphX(g, target, lower);
1370 
1371  }
1372  // here need to use the limit distribution
1373  SamplingDistribution * limitDist = GetLimitDistribution(lower);
1374  if (!limitDist) return 0;
1375  const std::vector<double> & values = limitDist->GetSamplingDistribution();
1376  double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
1377  TMath::Quantiles(values.size(), 1, x,q,p,false);
1378  return q[0];
1379 
1380 }
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
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:785
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:2173
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:628
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:467
#define ooccoutI(o, a)
Definition: RooMsgService.h:51
Double_t Floor(Double_t x)
Definition: TMath.h:599
virtual Double_t getMax(const char *name=0) const
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).
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)
Double_t QuietNaN()
Definition: TMath.h:783
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:44
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.
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
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:777
STL namespace.
#define ooccoutW(o, a)
Definition: RooMsgService.h:53
virtual TObject * RemoveAt(Int_t idx)
#define ooccoutE(o, a)
Definition: RooMsgService.h:54
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=1, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3950
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:1002
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:2365
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:47
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 Parameters: x -the data sample n ...
Definition: TMath.cxx:1178
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
virtual Double_t Eval(Double_t x, TSpline *spline=0, Option_t *option="") const
Interpolate points in this graph at x using a TSpline.
Definition: TGraph.cxx:863
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
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()
Definition: TMath.h:796
int ExclusionCleanup()
remove points that appear to have failed.
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition: TMath.h:318
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:322
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:3385
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:655
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:1447
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:90
SamplingDistribution * GetNullDistribution(void) const
Int_t GetN() const
Definition: TGraph.h:122
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:354
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()
Definition: TMath.h:74
virtual Double_t CLs() const
is simply (not a method, but a quantity)
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
Double_t * GetX() const
Definition: TGraph.h:129
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:31
virtual Double_t CLsplusb() const
Convert AlternatePValue into a "confidence level".
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:359
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
double Double_t
Definition: RtypesCore.h:55
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)
Definition: TMath.h:836
Double_t y[n]
Definition: legend1.C:17
static constexpr double s
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:50
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:46
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TObject.cxx:144
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.
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2184
SimpleInterval is a concrete implementation of the ConfInterval interface.
Double_t * GetY() const
Definition: TGraph.h:130
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:200
1-Dim function class
Definition: TF1.h:211
SamplingDistribution * GetNullTestStatDist(int index) const
same in terms of alt and null
A Graph is a graphics 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:487
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save Pad contents in a file in one of various formats.
Definition: TPad.cxx:5486
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
Definition: TCollection.h:180
float * q
Definition: THbookFile.cxx:87
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: TMath.h:1092
char name[80]
Definition: TGX11.cxx:109
void SetNullDistribution(SamplingDistribution *null)
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
Definition: TMath.h:1125
std::vector< double > fXValues
number of points used to build expected p-values
static constexpr double g