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