62 fInterpolateLowerLimit(true),
63 fInterpolateUpperLimit(true),
64 fFittedLowerLimit(false),
65 fFittedUpperLimit(false),
66 fInterpolOption(kLinear),
69 fCLsCleanupThreshold(0.005)
99 for (
int i = 0; i < nOther; ++i)
135 for (
int i=0; i < nOther; ++i) {
202 std::vector<double> qv;
211 bool resultIsAsymptotic(
false);
216 resultIsAsymptotic =
true;
220 int nPointsRemoved(0);
222 double CLsobsprev(1.0);
223 std::vector<double>::iterator itr =
fXValues.begin();
238 oocoutE(
this,
Eval) <<
"HypoTestInverterResult::ExclusionCleanup - invalid size of sampling distribution" << std::endl;
245 if (resultIsAsymptotic) {
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 );
260 double *
z =
const_cast<double *
>( &values[0] );
268 for (
int j=0; j<5; ++j) { qv[j]=q[j]; }
275 double CLsobs = qv[5];
279 bool removeThisPoint(
false);
282 if (!removeThisPoint && resultIsAsymptotic && i>=1 && CLsobs>CLsobsprev) {
284 removeThisPoint =
true;
285 }
else { CLsobsprev = CLsobs; }
288 if (!removeThisPoint && i>=1 && CLsobs>=0.9999) {
290 removeThisPoint =
true;
296 if (removeThisPoint) {
313 return nPointsRemoved;
331 if (nOther == 0)
return true;
339 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - merging result from " << otherResult.
GetName()
340 <<
" in " <<
GetName() << std::endl;
345 if (addExpPValues || mergeExpPValues)
346 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - merging also the expected p-values from pseudo-data" << std::endl;
353 for (
int i = 0; i < nOther; ++i)
360 for (
int i = 0; i < nOther; ++i) {
361 double otherVal = otherResult.
fXValues[i];
363 if (otherHTR == 0)
continue;
364 bool sameXFound =
false;
365 for (
int j = 0; j < nThis; ++j) {
372 thisHTR->
Append(otherHTR);
374 if (mergeExpPValues) {
379 if (thisNToys != otherNToys )
380 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::Add expected p values have been generated with different toys " << thisNToys <<
" , " << otherNToys << std::endl;
399 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - new number of points is " <<
fXValues.size()
402 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - new toys/point is " 424 if (!r)
return false;
461 return result->CLsplusb();
476 return result->CLsError();
478 return result->CLsplusbError();
503 return result->CLsplusb();
527 return result->CLbError();
539 return result->CLsplusbError();
551 return result->CLsError();
574 const double tol = 1.E-12;
593 std::cout <<
"using graph for search " << lowSearch <<
" min " << axmin <<
" max " << axmax << std::endl;
598 const double *
y = graph.
GetY();
599 int n = graph.
GetN();
601 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult::GetGraphX - need at least 2 points for interpolation (n=" << n <<
")\n";
602 return (n>0) ? y[0] : 0;
617 if (axmin >= axmax ) {
620 std::cout <<
"No rage given - check if extrapolation is needed " << std::endl;
623 xmin = graph.
GetX()[0];
624 xmax = graph.
GetX()[n-1];
626 double yfirst = graph.
GetY()[0];
627 double ylast = graph.
GetY()[n-1];
636 if ( (ymax < y0 && !lowSearch) || ( yfirst > y0 && lowSearch) ) {
640 if ( (ymax < y0 && lowSearch) || ( ylast > y0 && !lowSearch) ) {
645 auto func = [&](
double x) {
653 bool ret = brf.
Solve(100, 1.
E-16, 1.
E-6);
655 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult - interpolation failed - return inf" << std::endl;
658 double limit = brf.
Root();
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;
672 TString fname =
"graph_upper.root";
673 if (lowSearch) fname =
"graph_lower.root";
675 graph.
Write(
"graph");
682 if (axmin >= axmax) {
685 std::cout <<
"do new interpolation dividing from " << index <<
" and " << y[index] << std::endl;
688 if (lowSearch && index >= 1 && (y[0] - y0) * ( y[index]- y0) < 0) {
692 else if (!lowSearch && index < n-2 && (y[n-1] - y0) * ( y[index+1]- y0) < 0) {
694 limit =
GetGraphX(graph, y0, lowSearch, graph.
GetX()[index+1], graph.
GetX()[n-1] );
721 double val = (lowSearch) ? xmin : xmax;
722 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::FindInterpolatedLimit" 723 <<
" - not enough points to get the inverted interval - return " 732 std::vector<unsigned int> index(n );
736 for (
int i = 0; i <
n; ++i)
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];
754 std::cout <<
" max of y " << iymax <<
" " << xwithymax <<
" " << ymax <<
" target is " << target << std::endl;
761 xmin = ( graph.GetY()[0] <= target ) ? graph.GetX()[0] : varmin;
774 xmax = ( graph.GetY()[n-1] <= target ) ? graph.GetX()[n-1] : varmax;
786 if (iymax <= (n-1)/2 ) {
797 std::cout <<
" found xmin, xmax = " << xmin <<
" " << xmax <<
" for search " << lowSearch << std::endl;
808 if (upI < 1)
return xmin;
814 if (lowI >= n-1)
return xmax;
822 std::cout <<
"finding " << lowSearch <<
" limit between " << xmin <<
" " << xmax << endl;
826 double limit =
GetGraphX(graph, target, lowSearch, xmin, xmax);
832 TString limitType = (lowSearch) ?
"lower" :
"upper";
833 ooccoutD(
this,
Eval) <<
"HypoTestInverterResult::FindInterpolateLimit " 834 <<
"the computed " << limitType <<
" limit is " << limit <<
" +/- " << error << std::endl;
837 std::cout <<
"Found limit is " << limit <<
" +/- " << error << astd::endl;
885 int closestIndex = -1;
887 double smallestError = 2;
888 double bestValue = 2;
897 if ( dist < bestValue) {
902 if (bestIndex >=0)
return bestIndex;
910 std::vector<unsigned int> indx(n);
912 std::vector<double> xsorted( n);
913 for (
int i = 0; i <
n; ++i) xsorted[i] =
fXValues[indx[i] ];
917 std::cout <<
"finding closest point to " << xtarget <<
" is " << index1 <<
" " << indx[index1] << std::endl;
921 if (index1 < 0)
return indx[0];
922 if (index1 >= n-1)
return indx[n-1];
923 int index2 = index1 +1;
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];
974 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimateError" 975 <<
"Empty result \n";
980 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimateError" 981 <<
" only points - return its error\n";
988 TString
type = (!lower) ?
"upper" :
"lower";
991 std::cout <<
"calculate estimate error " << type <<
" between " << xmin <<
" and " << xmax << std::endl;
996 std::vector<unsigned int> indx(
fXValues.size());
1007 graph.SetPointError(ip, 0.,
GetYError(indx[i]) );
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;
1024 TF1 fct(
"fct",
"exp([0] * (x - [2] ) + [1] * (x-[2])**2)", minX, maxX);
1025 double scale = maxX-minX;
1038 double limit = (!lower) ? fUpperLimit :
fLowerLimit;
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);
1052 int fitstat = graph.Fit(&fct,
"Q EX0");
1056 double theError = 0;
1061 theError = std::min(
fabs(
GetYError(index) / m), maxX-minX);
1065 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimatedError - cannot estimate the " << type <<
" limit error " << std::endl;
1074 std::cout <<
"closes point to the limit is " << index <<
" " <<
GetXValue(index) <<
" and has error " <<
GetYError(index) << std::endl;
1107 if (!firstResult)
return 0;
1116 if (!result)
return 0;
1125 if (index < 0 || index >=
ArraySize() )
return 0;
1131 static bool useFirstB =
false;
1133 int bIndex = (useFirstB) ? 0 : index;
1140 if (bDistribution && sbDistribution) {
1150 std::vector<double> values(bDistribution->
GetSize());
1151 for (
int i = 0; i < bDistribution->
GetSize(); ++i) {
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;
1168 if (pval < 0) {
return 0;}
1171 return new SamplingDistribution(
"Asymptotic expected values",
"Asymptotic expected values",values);
1181 oocoutE(
this,
Eval) <<
"HypoTestInverterResult::GetLimitDistribution" 1182 <<
" not enough points - return 0 " << std::endl;
1186 ooccoutD(
this,
Eval) <<
"HypoTestInverterResult - computing limit distribution...." << std::endl;
1191 std::vector<SamplingDistribution*> distVec( npoints );
1193 for (
unsigned int i = 0; i < distVec.size(); ++i) {
1198 distVec[i]->InverseCDF(0);
1199 sum += distVec[i]->GetSize();
1202 int size = int( sum/ npoints);
1205 ooccoutW(
this,
InputArguments) <<
"HypoTestInverterResult - set a minimum size of 10 for limit distribution" << std::endl;
1213 std::vector< std::vector<double> > quantVec(npoints );
1214 for (
int i = 0; i < npoints; ++i) {
1216 if (!distVec[i])
continue;
1219 std::vector<double> pvalues = distVec[i]->GetSamplingDistribution();
1220 delete distVec[i]; distVec[i] = 0;
1221 std::sort(pvalues.begin(), pvalues.end());
1226 quantVec[i] = std::vector<double>(size);
1227 for (
int ibin = 0; ibin < size; ++ibin) {
1229 p[0] = std::min( (ibin+1) * 1./
double(size), 1.0);
1232 (quantVec[i])[ibin] = q[0];
1238 std::vector<unsigned int> index( npoints );
1245 std::vector<double> limits(size);
1247 for (
int j = 0; j < size; ++j ) {
1250 for (
int k = 0; k < npoints ; ++k) {
1251 if (quantVec[index[k]].size() > 0 )
1255 limits[j] =
GetGraphX(g, target, lower);
1313 if (nEntries <= 0)
return (lower) ? 1 : 0;
1321 if (!limitDist)
return 0;
1323 if (values.size() <= 1)
return 0;
1338 TString option(opt);
1340 if (option.Contains(
"P")) {
1345 std::vector<unsigned int> index(nEntries);
1348 for (
int j=0; j<nEntries; ++j) {
1352 ooccoutI(
this,
Eval) <<
"HypoTestInverterResult - cannot compute expected p value distribution for point, x = " 1353 <<
GetXValue(i) <<
" skip it " << std::endl;
1357 double *
x =
const_cast<double *
>(&values[0]);
1363 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult - cannot compute limits , not enough points, n = " << g.
GetN() << std::endl;
1374 if (!limitDist)
return 0;
1376 double *
x =
const_cast<double *
>(&values[0]);
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.
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
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)
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
virtual void SetParameters(const Double_t *params)
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
get the signal and background test statistic distribution
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 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).
SamplingDistribution * GetExpectedPValueDist(int index) const
return expected distribution of p-values (Cls or Clsplusb)
virtual TObject * clone(const char *newname) const
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 ...
int FindClosestPointIndex(double target, int mode=0, double xtarget=0)
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
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...
virtual TObject * RemoveAt(Int_t idx)
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=1, Int_t netopt=0)
Create / open a file.
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...
SimpleInterval & operator=(const SimpleInterval &other)
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)
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
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 ...
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 ~HypoTestInverterResult()
destructor
virtual Double_t Eval(Double_t x, TSpline *spline=0, Option_t *option="") const
Interpolate points in this graph at x using a TSpline.
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)
int ExclusionCleanup()
remove points that appear to have failed.
RooRealVar represents a fundamental (non-derived) real valued object.
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
double CLs(int index) const
return the observed CLb value for the i-th entry
double Root() const
Returns root value.
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)
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
SamplingDistribution * GetBackgroundTestStatDist(int index) const
get the background test statistic distribution
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
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
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...
SamplingDistribution * GetNullDistribution(void) const
static int fgAsymptoticNumPoints
max sigma value used to scan asymptotic expected p values
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
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
virtual Double_t CLs() const
is simply (not a method, but a quantity)
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
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 CLsplusb() const
Convert AlternatePValue into a "confidence level".
Namespace for the RooStats classes.
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
Double_t LowerLimitEstimatedError()
rough estimation of the error on the computed bound of the confidence interval Estimate of lower limi...
T MaxElement(Long64_t n, const T *a)
InterpolOption_t fInterpolOption
double func(double *x, double *p)
double GetYValue(int index) const
function to return the value of the confidence level for the i^th entry in the results ...
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
you should not use this method at all Int_t Int_t z
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.
SimpleInterval is a concrete implementation of the ConfInterval interface.
HypoTestInverterResult & operator=(const HypoTestInverterResult &other)
operator =
virtual void Add(TObject *obj)
double fLowerLimitError
interpolation option (linear or spline)
virtual void RemoveAll(TCollection *col)
Remove all objects in collection col from this collection.
Short_t Max(Short_t a, Short_t b)
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.
A TGraphErrors is a TGraph with error bars.
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.
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save Pad contents in a file in one of various formats.
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
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
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)
std::vector< double > fXValues
number of points used to build expected p-values