49 fInterpolateLowerLimit(true),
50 fInterpolateUpperLimit(true),
51 fFittedLowerLimit(false),
52 fFittedUpperLimit(false),
53 fInterpolOption(kLinear),
56 fCLsCleanupThreshold(0.005)
63 fExpPValues.SetOwner();
69 fUseCLs(other.fUseCLs),
70 fIsTwoSided(other.fIsTwoSided),
71 fInterpolateLowerLimit(other.fInterpolateLowerLimit),
72 fInterpolateUpperLimit(other.fInterpolateUpperLimit),
73 fFittedLowerLimit(other.fFittedLowerLimit),
74 fFittedUpperLimit(other.fFittedUpperLimit),
75 fInterpolOption(other.fInterpolOption),
76 fLowerLimitError(other.fLowerLimitError),
77 fUpperLimitError(other.fUpperLimitError),
78 fCLsCleanupThreshold(other.fCLsCleanupThreshold)
86 for (
int i = 0; i < nOther; ++i)
121 for (
int i=0; i < nOther; ++i) {
183 std::vector<double> qv;
192 bool resultIsAsymptotic(
false);
197 resultIsAsymptotic =
true;
201 int nPointsRemoved(0);
203 double CLsobsprev(1.0);
204 std::vector<double>::iterator itr =
fXValues.begin();
219 oocoutE(
this,
Eval) <<
"HypoTestInverterResult::ExclusionCleanup - invalid size of sampling distribution" << std::endl;
226 if (resultIsAsymptotic) {
228 double dsig = 2.*maxSigma / (values.size() -1) ;
229 int i0 = (int)
TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5 );
230 int i1 = (int)
TMath::Floor ( ( -nsig1 + maxSigma )/dsig + 0.5 );
231 int i2 = (int)
TMath::Floor ( ( maxSigma )/dsig + 0.5 );
232 int i3 = (int)
TMath::Floor ( ( nsig1 + maxSigma )/dsig + 0.5 );
233 int i4 = (int)
TMath::Floor ( ( nsig2 + maxSigma )/dsig + 0.5 );
241 double *
z =
const_cast<double *
>( &values[0] );
249 for (
int j=0; j<5; ++j) { qv[j]=q[j]; }
256 double CLsobs = qv[5];
260 bool removeThisPoint(
false);
263 if (!removeThisPoint && resultIsAsymptotic && i>=1 && CLsobs>CLsobsprev) {
265 removeThisPoint =
true;
266 }
else { CLsobsprev = CLsobs; }
269 if (!removeThisPoint && i>=1 && CLsobs>=0.9999) {
271 removeThisPoint =
true;
277 if (removeThisPoint) {
294 return nPointsRemoved;
312 if (nOther == 0)
return true;
320 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - merging result from " << otherResult.
GetName()
321 <<
" in " <<
GetName() << std::endl;
326 if (addExpPValues || mergeExpPValues)
327 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - merging also the expected p-values from pseudo-data" << std::endl;
334 for (
int i = 0; i < nOther; ++i)
341 for (
int i = 0; i < nOther; ++i) {
342 double otherVal = otherResult.
fXValues[i];
344 if (otherHTR == 0)
continue;
345 bool sameXFound =
false;
346 for (
int j = 0; j < nThis; ++j) {
353 thisHTR->
Append(otherHTR);
355 if (mergeExpPValues) {
360 if (thisNToys != otherNToys )
361 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::Add expexted p values have been generated with different toys " << thisNToys <<
" , " << otherNToys << std::endl;
380 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - new number of points is " <<
fXValues.size()
383 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - new toys/point is " 403 if (!r)
return false;
438 return result->CLsplusb();
451 return result->CLsError();
453 return result->CLsplusbError();
474 return result->CLsplusb();
493 return result->CLbError();
503 return result->CLsplusbError();
513 return result->CLsError();
532 const double tol = 1.E-12;
543 struct InterpolatedGraph {
544 InterpolatedGraph(
const TGraph & g,
double target,
const char * interpOpt) :
545 fGraph(g), fTarget(target), fInterpOpt(interpOpt) {}
548 double operator() (
double x)
const {
549 return fGraph.Eval(x, (
TSpline*) 0, fInterpOpt) - fTarget;
566 InterpolatedGraph
f(graph,y0,opt);
571 const double *
y = graph.
GetY();
572 int n = graph.
GetN();
574 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult::GetGraphX - need at least 2 points for interpolation (n=" << n <<
")\n";
575 return (n>0) ? y[0] : 0;
589 if (axmin >= axmax ) {
590 xmin = graph.
GetX()[0];
591 xmax = graph.
GetX()[n-1];
599 if ( (ymax < y0 && !lowSearch) || ( ymin > y0 && lowSearch) ) {
603 if ( (ymax < y0 && lowSearch) || ( ymin > y0 && !lowSearch) ) {
609 #ifdef ISREALLYNEEDED //?? 611 bool isCross =
false;
614 for (
int i = 0; i<
n; ++i) {
617 if (xp >= xmin && xp <= xmax) {
619 if (prod < 0 && !first) {
627 return (lowSearch) ? varmin : varmax;
636 bool ret = brf.
Solve(100, 1.
E-16, 1.
E-6);
638 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult - interpolation failed - return inf" << std::endl;
641 double limit = brf.
Root();
645 if (lowSearch) std::cout <<
"lower limit search : ";
646 else std::cout <<
"Upper limit search : ";
647 std::cout <<
"do interpolation between " << xmin <<
" and " << xmax <<
" limit is " << limit << std::endl;
652 if (axmin >= axmax) {
655 std::cout <<
"do new interpolation dividing from " << index <<
" and " << y[index] << std::endl;
658 if (lowSearch && index >= 1 && (y[0] - y0) * ( y[index]- y0) < 0) {
662 else if (!lowSearch && index < n-2 && (y[n-1] - y0) * ( y[index+1]- y0) < 0) {
664 limit =
GetGraphX(graph, y0, lowSearch, graph.
GetX()[index+1], graph.
GetX()[n-1] );
680 <<
"Interpolate the upper limit between the 2 results closest to the target confidence level" 693 double val = (lowSearch) ? xmin : xmax;
694 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::FindInterpolatedLimit" 695 <<
" - not enough points to get the inverted interval - return " 704 std::vector<unsigned int> index(n );
708 for (
int i = 0; i <
n; ++i)
720 double * itrmax = std::max_element(graph.GetY() , graph.GetY() +
n);
721 double ymax = *itrmax;
722 int iymax = itrmax - graph.GetY();
723 double xwithymax = graph.GetX()[iymax];
732 xmin = ( graph.GetY()[0] <= target ) ? graph.GetX()[0] : varmin;
745 xmax = ( graph.GetY()[n-1] <= target ) ? graph.GetX()[n-1] : varmax;
759 if (iymax <= (n-1)/2 ) {
770 std::cout <<
" found xmin, xmax = " << xmin <<
" " << xmax <<
" for search " << lowSearch << std::endl;
781 if (upI < 1)
return xmin;
787 if (lowI >= n-1)
return xmax;
795 std::cout <<
"finding " << lowSearch <<
" limit betweem " << xmin <<
" " << xmax << endl;
798 double limit =
GetGraphX(graph, target, lowSearch, xmin, xmax);
804 std::cout <<
"limit is " << limit << std::endl;
818 double limit2 =
GetGraphX(graph, target, !lowSearch, xmin, xmax);
825 std::cout <<
"other limit is " << limit2 << std::endl;
844 int closestIndex = -1;
846 double smallestError = 2;
847 double bestValue = 2;
856 if ( dist < bestValue) {
861 if (bestIndex >=0)
return bestIndex;
869 std::vector<unsigned int> indx(n);
871 std::vector<double> xsorted( n);
872 for (
int i = 0; i <
n; ++i) xsorted[i] =
fXValues[indx[i] ];
876 std::cout <<
"finding closest point to " << xtarget <<
" is " << index1 <<
" " << indx[index1] << std::endl;
880 if (index1 < 0)
return indx[0];
881 if (index1 >= n-1)
return indx[n-1];
882 int index2 = index1 +1;
884 if (mode == 2)
return (
GetXValue(indx[index1]) <
GetXValue(indx[index2])) ? indx[index1] : indx[index2];
885 if (mode == 3)
return (
GetXValue(indx[index1]) >
GetXValue(indx[index2])) ? indx[index1] : indx[index2];
928 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimateError" 929 <<
"Empty result \n";
934 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimateError" 935 <<
" only points - return its error\n";
945 std::cout <<
"calculate estimate error " << type <<
" between " << xmin <<
" and " << xmax << std::endl;
950 std::vector<unsigned int> indx(
fXValues.size());
961 graph.SetPointError(ip, 0.,
GetYError(indx[i]) );
966 if (graph.GetN() < 2) {
967 if (np >= 2)
oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimatedError - no valid points - cannot estimate the " << type <<
" limit error " << std::endl;
980 TF1 fct(
"fct",
"exp([0] * (x - [2] ) + [1] * (x-[2])**2)", minX, maxX);
981 double scale = maxX-minX;
994 double limit = (!lower) ? fUpperLimit :
fLowerLimit;
1000 std::cout <<
"fitting for limit " << type <<
"between " << minX <<
" , " << maxX <<
" points considered " << graph.GetN() << std::endl;
1001 int fitstat = graph.Fit(&fct,
" EX0");
1002 graph.SetMarkerStyle(20);
1008 int fitstat = graph.Fit(&fct,
"Q EX0");
1012 double theError = 0;
1017 theError = std::min(
fabs(
GetYError(index) / m), maxX-minX);
1021 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimatedError - cannot estimate the " << type <<
" limit error " << std::endl;
1030 std::cout <<
"closes point to the limit is " << index <<
" " <<
GetXValue(index) <<
" and has error " <<
GetYError(index) << std::endl;
1059 if (!firstResult)
return 0;
1066 if (!result)
return 0;
1073 if (index < 0 || index >=
ArraySize() )
return 0;
1079 static bool useFirstB =
false;
1081 int bIndex = (useFirstB) ? 0 : index;
1088 if (bDistribution && sbDistribution) {
1098 std::vector<double> values(bDistribution->
GetSize());
1099 for (
int i = 0; i < bDistribution->
GetSize(); ++i) {
1111 const double dsig = 2* smax/ (npoints-1) ;
1112 std::vector<double> values(npoints);
1113 for (
int i = 0; i < npoints; ++i) {
1114 double nsig = -smax + dsig*i;
1116 if (pval < 0) {
return 0;}
1119 return new SamplingDistribution(
"Asymptotic expected values",
"Asymptotic expected values",values);
1131 oocoutE(
this,
Eval) <<
"HypoTestInverterResult::GetLimitDistribution" 1132 <<
" not enought points - return 0 " << std::endl;
1136 ooccoutD(
this,
Eval) <<
"HypoTestInverterResult - computing limit distribution...." << std::endl;
1141 std::vector<SamplingDistribution*> distVec( npoints );
1143 for (
unsigned int i = 0; i < distVec.size(); ++i) {
1148 distVec[i]->InverseCDF(0);
1149 sum += distVec[i]->GetSize();
1152 int size = int( sum/ npoints);
1155 ooccoutW(
this,
InputArguments) <<
"HypoTestInverterResult - set a minimum size of 10 for limit distribution" << std::endl;
1163 std::vector< std::vector<double> > quantVec(npoints );
1164 for (
int i = 0; i < npoints; ++i) {
1166 if (!distVec[i])
continue;
1169 std::vector<double> pvalues = distVec[i]->GetSamplingDistribution();
1170 delete distVec[i]; distVec[i] = 0;
1171 std::sort(pvalues.begin(), pvalues.end());
1176 quantVec[i] = std::vector<double>(size);
1177 for (
int ibin = 0; ibin < size; ++ibin) {
1179 p[0] = std::min( (ibin+1) * 1./
double(size), 1.0);
1182 (quantVec[i])[ibin] = q[0];
1188 std::vector<unsigned int> index( npoints );
1195 std::vector<double> limits(size);
1197 for (
int j = 0; j < size; ++j ) {
1200 for (
int k = 0; k < npoints ; ++k) {
1201 if (quantVec[index[k]].size() > 0 )
1205 limits[j] =
GetGraphX(g, target, lower);
1261 if (nEntries <= 0)
return (lower) ? 1 : 0;
1269 if (!limitDist)
return 0;
1271 if (values.size() <= 1)
return 0;
1293 std::vector<unsigned int> index(nEntries);
1296 for (
int j=0; j<nEntries; ++j) {
1300 ooccoutI(
this,
Eval) <<
"HypoTestInverterResult - cannot compute expected p value distribution for point, x = " 1301 <<
GetXValue(i) <<
" skip it " << std::endl;
1305 double *
x =
const_cast<double *
>(&values[0]);
1311 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult - cannot compute limits , not enough points, n = " << g.
GetN() << std::endl;
1322 if (!limitDist)
return 0;
1324 double *
x =
const_cast<double *
>(&values[0]);
SamplingDistribution * GetLimitDistribution(bool lower) const
void SetBackgroundAsAlt(Bool_t l=kTRUE)
double dist(Rotation3D const &r1, Rotation3D const &r2)
void SetAltDistribution(SamplingDistribution *alt)
static long int sum(long int i)
virtual void SetParameters(const Double_t *params)
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Double_t Floor(Double_t x)
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).
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
int ArraySize() const
number of entries in the results array
virtual Double_t CLs() const
is simply (not a method, but a quantity)
int FindClosestPointIndex(double target, int mode=0, double xtarget=0)
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 //.
HypoTestInverterResult(const char *name=0)
default constructor
HypoTestResult is a base class for results from hypothesis tests.
void ToUpper()
Change string to upper case.
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
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".
SamplingDistribution * GetBackgroundTestStatDist(int index) const
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.
virtual TObject * RemoveAt(Int_t idx)
Template class to wrap any C++ callable object which takes one argument i.e.
const char * Data() const
RooAbsArg * first() const
SimpleInterval & operator=(const SimpleInterval &other)
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
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
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Add element to an owning set.
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 ...
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
Double_t fConfidenceLevel
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save Pad contents in a file in one of various formats.
virtual ~HypoTestInverterResult()
destructor
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.
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)
int ExclusionCleanup()
remove points that appear to have failed.
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
RooRealVar represents a fundamental (non-derived) real valued object.
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)
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
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 ...
double fCLsCleanupThreshold
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. ...
TList fExpPValues
list of HypoTestResult for each point
virtual const char * GetName() const
Returns name of object.
static int fgAsymptoticNumPoints
max sigma value used to scan asymptotic expected p values
Bool_t GetPValueIsRightTail(void) const
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
static double fgAsymptoticMaxSigma
Class for finding the root of a one dimensional function using the Brent algorithm.
void SetPValueIsRightTail(Bool_t pr)
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's extrapolation metho...
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results ...
Namespace for the RooStats classes.
virtual Int_t GetSize() const
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_t LowerLimitEstimatedError()
rough estimation of the error on the computed bound of the confidence interval Estimate of lower limi...
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)
InterpolOption_t fInterpolOption
SamplingDistribution * GetNullDistribution(void) const
virtual Int_t GetPoint(Int_t i, Double_t &x, Double_t &y) const
Get x and y values for point number i.
double CLsplusbError(int index) const
return the observed CLsplusb value for the i-th entry
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
you should not use this method at all Int_t Int_t z
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.
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
double GetGraphX(const TGraph &g, double y0, bool lowSearch, double &xmin, double &xmax) const
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
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)
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
A Graph is a graphics object made of two arrays X and Y with npoints each.
A TGraphErrors is a TGraph with error bars.
double GetExpectedLimit(double nsig, bool lower, const char *opt="") const
bool Add(const HypoTestInverterResult &otherResult)
merge with the content of another HypoTestInverterResult object
int FindIndex(double xvalue) const
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)
double CLb(int index) const
return the observed CLb value for the i-th entry
bool fInterpolateUpperLimit
Long64_t BinarySearch(Long64_t n, const T *array, T value)
void SetNullDistribution(SamplingDistribution *null)
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
SamplingDistribution * GetNullTestStatDist(int index) const
same in terms of alt and null
std::vector< double > fXValues
number of points used to build expected p-values
T MinElement(Long64_t n, const T *a)
Bool_t GetBackGroundIsAlt(void) const