| 1 |
#include <cmath>
|
| 2 |
#include <iostream>
|
| 3 |
#include "solve.h"
|
| 4 |
#include "cls.h"
|
| 5 |
#include "significance.h"
|
| 6 |
|
| 7 |
|
| 8 |
template <class TClass>
|
| 9 |
double TSolve<TClass>::_bisection(int n)
|
| 10 |
{
|
| 11 |
int i=0;
|
| 12 |
double a=_left, b=_right;
|
| 13 |
do {
|
| 14 |
double xmid = a+(b-a)*0.5;
|
| 15 |
double mid = (*_pt2Object.*_function)(xmid,_par);
|
| 16 |
//std::cout<<"TSolve<TClass>::bisection(left="<<a<<", right="<<b<<"); f(mid)="<<mid<<std::endl;
|
| 17 |
if ((*_pt2Object.*_function)(a,_par)*mid<0.) b = xmid;
|
| 18 |
if ((*_pt2Object.*_function)(b,_par)*mid<0.) a = xmid;
|
| 19 |
if (b-a<epsilon) return xmid;
|
| 20 |
} while ( i<n );
|
| 21 |
return 0;
|
| 22 |
}
|
| 23 |
|
| 24 |
template <class TClass>
|
| 25 |
double TSolve<TClass>::operator()()
|
| 26 |
{
|
| 27 |
//std::cout<<"TSolve<TClass>::operator() " <<std::endl;
|
| 28 |
if (!_isBracketed(_left,_right)) BracketFunction(_left,_right);
|
| 29 |
return _bisection();
|
| 30 |
}
|
| 31 |
|
| 32 |
template <class TClass>
|
| 33 |
bool TSolve<TClass>::_isBracketed(double l, double r)
|
| 34 |
{
|
| 35 |
//std::cerr<<"TSolve<TClass>::isBracketed(double l="<<l<<", double r="<<r<<") ? "
|
| 36 |
// <<"CLs(l*sigma) ="<<(*_pt2Object.*_function)(l,_par)+0.05
|
| 37 |
// <<" <-> CLs(r*sigma)="<< (*_pt2Object.*_function)(r,_par)+0.05<<std::endl;
|
| 38 |
return (*_pt2Object.*_function)(l,_par) * (*_pt2Object.*_function)(r,_par) <= 0.0;
|
| 39 |
}
|
| 40 |
|
| 41 |
|
| 42 |
template <class TClass>
|
| 43 |
void TSolve<TClass>::BracketFunction(double l, double r)
|
| 44 |
{
|
| 45 |
//std::cout<<"TSolve<TClass>::BracketFunction(double l="<<l<<", double r="<<r<<")"<<std::endl;
|
| 46 |
//double ix = epsilon;
|
| 47 |
double ix = 1;
|
| 48 |
int i=0;
|
| 49 |
do {
|
| 50 |
if (i>5) throw TSolveException("Not converged! Cannot be bracketed!");
|
| 51 |
|
| 52 |
ix *= 2.;
|
| 53 |
++i;
|
| 54 |
// } while(!_isBracketed(l-ix,r+ix) && ix<100000);
|
| 55 |
// _left = l-ix;
|
| 56 |
// _right = r+ix;
|
| 57 |
} while(!_isBracketed(l/ix,r*ix) && ix<100000 );
|
| 58 |
_left = l/ix;
|
| 59 |
_right = r*ix;
|
| 60 |
}
|
| 61 |
|
| 62 |
template class TSolve<cls>;
|
| 63 |
template class TSolve<TSignificance>;
|
