#include "TRolke.h"
#include "TMath.h"
#include "Riostream.h"
ClassImp(TRolke)
TRolke::TRolke(Double_t CL, Option_t * )
{
fUpperLimit = 0.0;
fLowerLimit = 0.0;
fCL = CL;
fSwitch = 0;
}
TRolke::~TRolke()
{
}
Double_t TRolke::CalculateInterval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em,Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
{
Int_t done = 0;
Double_t limit[2];
limit[1] = Interval(x,y,z,bm,em,e,mid, sde,sdb,tau,b,m);
if (limit[1] > 0) {
done = 1;
}
if (fSwitch == 0) {
Int_t trial_x = x;
while (done == 0) {
trial_x++;
limit[1] = Interval(trial_x,y,z,bm,em,e,mid, sde,sdb,tau,b,m);
if (limit[1] > 0) done = 1;
}
}
return limit[1];
}
Double_t TRolke::Interval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em,Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
{
Double_t dchi2 = TMath::ChisquareQuantile(fCL, 1);
Double_t tempxy[2],limits[2] = {0,0};
Double_t slope,fmid,low,flow,high,fhigh,test,ftest,mu0,maximum,target,l,f0;
Double_t med = 0;
Double_t maxiter=1000, acc = 0.00001;
Int_t i;
Int_t bp = 0;
if ((mid != 3) && (mid != 5)) bm = (Double_t)y;
if ((mid == 3) || (mid == 5)) {
if (bm == 0) bm = 0.00001;
}
if ((mid <= 2) || (mid == 4)) bp = 1;
if (bp == 1 && x == 0 && bm > 0 ){
for(Int_t i = 0; i < 2; i++) {
x++;
tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
}
slope = tempxy[1] - tempxy[0];
limits[1] = tempxy[0] - slope;
limits[0] = 0.0;
if (limits[1] < 0) limits[1] = 0.0;
goto done;
}
if (bp != 1 && x == 0){
for(Int_t i = 0; i < 2; i++) {
x++;
tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
}
slope = tempxy[1] - tempxy[0];
limits[1] = tempxy[0] - slope;
limits[0] = 0.0;
if (limits[1] < 0) limits[1] = 0.0;
goto done;
}
if (bp != 1 && bm == 0){
for(Int_t i = 0; i < 2; i++) {
bm++;
limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
tempxy[i] = limits[1];
}
slope = tempxy[1] - tempxy[0];
limits[1] = tempxy[0] - slope;
if (limits[1] < 0) limits[1] = 0;
goto done;
}
if (x == 0 && bm == 0){
x++;
bm++;
limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
tempxy[0] = limits[1];
x = 1;
bm = 2;
limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
tempxy[1] = limits[1];
x = 2;
bm = 1;
limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
limits[1] = 3*tempxy[0] -tempxy[1] - limits[1];
if (limits[1] < 0) limits[1] = 0;
goto done;
}
mu0 = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,1);
maximum = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,2);
test = 0;
f0 = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
if ( fSwitch == 1 ) {
if ( mu0 < 0 ) maximum = f0;
}
target = maximum - dchi2;
if (f0 > target) {
limits[0] = 0;
} else {
if (mu0 < 0){
limits[0] = 0;
limits[1] = 0;
}
low = 0;
flow = f0;
high = mu0;
fhigh = maximum;
for(Int_t i = 0; i < maxiter; i++) {
l = (target-fhigh)/(flow-fhigh);
if (l < 0.2) l = 0.2;
if (l > 0.8) l = 0.8;
med = l*low + (1-l)*high;
if(med < 0.01){
limits[1]=0.0;
goto done;
}
fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
if (fmid > target) {
high = med;
fhigh = fmid;
} else {
low = med;
flow = fmid;
}
if ((high-low) < acc*high) break;
}
limits[0] = med;
}
if(mu0 > 0) {
low = mu0;
flow = maximum;
} else {
low = 0;
flow = f0;
}
test = low +1 ;
ftest = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
if (ftest < target) {
high = test;
fhigh = ftest;
} else {
slope = (ftest - flow)/(test - low);
high = test + (target -ftest)/slope;
fhigh = Likelihood(high,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
}
for(i = 0; i < maxiter; i++) {
l = (target-fhigh)/(flow-fhigh);
if (l < 0.2) l = 0.2;
if (l > 0.8) l = 0.8;
med = l * low + (1.-l)*high;
fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
if (fmid < target) {
high = med;
fhigh = fmid;
} else {
low = med;
flow = fmid;
}
if (high-low < acc*high) break;
}
limits[1] = med;
done:
if ( (mid == 4) || (mid==5) ) {
limits[0] /= e;
limits[1] /= e;
}
fUpperLimit = limits[1];
fLowerLimit = TMath::Max(limits[0],0.0);
return limits[1];
}
Double_t TRolke::Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm,Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m, Int_t what)
{
switch (mid) {
case 1: return EvalLikeMod1(mu,x,y,z,e,tau,b,m,what);
case 2: return EvalLikeMod2(mu,x,y,em,e,sde,tau,b,what);
case 3: return EvalLikeMod3(mu,x,bm,em,e,sde,sdb,b,what);
case 4: return EvalLikeMod4(mu,x,y,tau,b,what);
case 5: return EvalLikeMod5(mu,x,bm,sdb,b,what);
case 6: return EvalLikeMod6(mu,x,z,e,b,m,what);
case 7: return EvalLikeMod7(mu,x,em,e,sde,b,what);
}
return 0;
}
Double_t TRolke::EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t e, Double_t tau, Double_t b, Int_t m, Int_t what)
{
Double_t f = 0;
Double_t zm = Double_t(z)/m;
if (what == 1) {
f = (x-y/tau)/zm;
}
if (what == 2) {
mu = (x-y/tau)/zm;
b = y/tau;
Double_t e = zm;
f = LikeMod1(mu,b,e,x,y,z,tau,m);
}
if (what == 3) {
if (mu == 0){
b = (x+y)/(1.0+tau);
e = zm;
f = LikeMod1(mu,b,e,x,y,z,tau,m);
} else {
TRolke g;
g.ProfLikeMod1(mu,b,e,x,y,z,tau,m);
f = LikeMod1(mu,b,e,x,y,z,tau,m);
}
}
return f;
}
Double_t TRolke::LikeMod1(Double_t mu,Double_t b, Double_t e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
{
return 2*(x*TMath::Log(e*mu+b)-(e*mu +b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) + TMath::Log(TMath::Factorial(m)) - TMath::Log(TMath::Factorial(m-z)) - TMath::Log(TMath::Factorial(z))+ z * TMath::Log(e) + (m-z)*TMath::Log(1-e));
}
void TRolke::ProfLikeMod1(Double_t mu,Double_t &b,Double_t &e,Int_t x,Int_t y, Int_t z,Double_t tau,Int_t m)
{
Double_t med = 0.0,fmid;
Int_t maxiter =1000;
Double_t acc = 0.00001;
Double_t emin = ((m+mu*tau)-TMath::Sqrt((m+mu*tau)*(m+mu*tau)-4 * mu* tau * z))/2/mu/tau;
Double_t low = TMath::Max(1e-10,emin+1e-10);
Double_t high = 1 - 1e-10;
for(Int_t i = 0; i < maxiter; i++) {
med = (low+high)/2.;
fmid = LikeGradMod1(med,mu,x,y,z,tau,m);
if(high < 0.5) acc = 0.00001*high;
else acc = 0.00001*(1-high);
if ((high - low) < acc*high) break;
if(fmid > 0) low = med;
else high = med;
}
e = med;
Double_t eta = Double_t(z)/e -Double_t(m-z)/(1-e);
b = Double_t(y)/(tau -eta/mu);
}
Double_t TRolke::LikeGradMod1(Double_t e, Double_t mu, Int_t x,Int_t y,Int_t z,Double_t tau,Int_t m)
{
Double_t eta, etaprime, bprime,f;
eta = static_cast<double>(z)/e - static_cast<double>(m-z)/(1.0 - e);
etaprime = (-1) * (static_cast<double>(m-z)/((1.0 - e)*(1.0 - e)) + static_cast<double>(z)/(e*e));
Double_t b = y/(tau - eta/mu);
bprime = (b*b * etaprime)/mu/y;
f = (mu + bprime) * (x/(e * mu + b) - 1)+(y/b - tau) * bprime + eta;
return f;
}
Double_t TRolke::EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t e,Double_t sde, Double_t tau, Double_t b, Int_t what)
{
Double_t v = sde*sde;
Double_t coef[4],roots[3];
Double_t f = 0;
if (what == 1) {
f = (x-y/tau)/em;
}
if (what == 2) {
mu = (x-y/tau)/em;
b = y/tau;
e = em;
f = LikeMod2(mu,b,e,x,y,em,tau,v);
}
if (what == 3) {
if (mu == 0 ) {
b = (x+y)/(1+tau);
f = LikeMod2(mu,b,e,x,y,em,tau,v);
} else {
coef[3] = mu;
coef[2] = mu*mu*v-2*em*mu-mu*mu*v*tau;
coef[1] = ( - x)*mu*v - mu*mu*mu*v*v*tau - mu*mu*v*em + em*mu*mu*v*tau + em*em*mu - y*mu*v;
coef[0] = x*mu*mu*v*v*tau + x*em*mu*v - y*mu*mu*v*v + y*em*mu*v;
TMath::RootsCubic(coef,roots[0],roots[1],roots[2]);
e = roots[1];
b = y/(tau + (em - e)/mu/v);
f = LikeMod2(mu,b,e,x,y,em,tau,v);
}
}
return f;
}
Double_t TRolke::LikeMod2(Double_t mu, Double_t b, Double_t e,Int_t x,Int_t y,Double_t em,Double_t tau, Double_t v)
{
return 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
}
Double_t TRolke::EvalLikeMod3(Double_t mu, Int_t x, Double_t bm, Double_t em, Double_t e, Double_t sde, Double_t sdb, Double_t b, Int_t what)
{
Double_t f = 0.;
Double_t v = sde*sde;
Double_t u = sdb*sdb;
if (what == 1) {
f = (x-bm)/em;
}
if (what == 2) {
mu = (x-bm)/em;
b = bm;
e = em;
f = LikeMod3(mu,b,e,x,bm,em,u,v);
}
if(what == 3) {
if(mu == 0.0){
b = ((bm-u)+TMath::Sqrt((bm-u)*(bm-u)+4*x*u))/2.;
e = em;
f = LikeMod3(mu,b,e,x,bm,em,u,v);
} else {
Double_t temp[3];
temp[0] = mu*mu*v+u;
temp[1] = mu*mu*mu*v*v+mu*v*u-mu*mu*v*em+mu*v*bm-2*u*em;
temp[2] = mu*mu*v*v*bm-mu*v*u*em-mu*v*bm*em+u*em*em-mu*mu*v*v*x;
e = (-temp[1]+TMath::Sqrt(temp[1]*temp[1]-4*temp[0]*temp[2]))/2/temp[0];
b = bm-(u*(em-e))/v/mu;
f = LikeMod3(mu,b,e,x,bm,em,u,v);
}
}
return f;
}
Double_t TRolke::LikeMod3(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t bm,Double_t em,Double_t u,Double_t v)
{
return 2*(x * TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)-1.837877-TMath::Log(u)/2-(bm-b)*(bm-b)/u/2-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
}
Double_t TRolke::EvalLikeMod4(Double_t mu, Int_t x, Int_t y, Double_t tau, Double_t b, Int_t what)
{
Double_t f = 0.0;
if (what == 1) f = x-y/tau;
if (what == 2) {
mu = x-y/tau;
b = Double_t(y)/tau;
f = LikeMod4(mu,b,x,y,tau);
}
if (what == 3) {
if (mu == 0.0) {
b = Double_t(x+y)/(1+tau);
f = LikeMod4(mu,b,x,y,tau);
} else {
b = (x+y-(1+tau)*mu+sqrt((x+y-(1+tau)*mu)*(x+y-(1+tau)*mu)+4*(1+tau)*y*mu))/2/(1+tau);
f = LikeMod4(mu,b,x,y,tau);
}
}
return f;
}
Double_t TRolke::LikeMod4(Double_t mu,Double_t b,Int_t x,Int_t y,Double_t tau)
{
return 2*(x*TMath::Log(mu+b)-(mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) );
}
Double_t TRolke::EvalLikeMod5(Double_t mu, Int_t x, Double_t bm, Double_t sdb, Double_t b, Int_t what)
{
Double_t u=sdb*sdb;
Double_t f = 0;
if(what == 1) {
f = x - bm;
}
if(what == 2) {
mu = x-bm;
b = bm;
f = LikeMod5(mu,b,x,bm,u);
}
if (what == 3) {
b = ((bm-u-mu)+TMath::Sqrt((bm-u-mu)*(bm-u-mu)-4*(mu*u-mu*bm-u*x)))/2;
f = LikeMod5(mu,b,x,bm,u);
}
return f;
}
Double_t TRolke::LikeMod5(Double_t mu,Double_t b,Int_t x,Double_t bm,Double_t u)
{
return 2*(x*TMath::Log(mu+b)-(mu + b)-LogFactorial(x)-0.9189385-TMath::Log(u)/2-((bm-b)*(bm-b))/u/2);
}
Double_t TRolke::EvalLikeMod6(Double_t mu, Int_t x, Int_t z, Double_t e, Double_t b, Int_t m, Int_t what)
{
Double_t coef[4],roots[3];
Double_t f = 0.;
Double_t zm = Double_t(z)/m;
if(what==1){
f = (x-b)/zm;
}
if(what==2){
mu = (x-b)/zm;
e = zm;
f = LikeMod6(mu,b,e,x,z,m);
}
if(what == 3){
if(mu==0){
e = zm;
} else {
coef[3] = mu*mu;
coef[2] = mu * b - mu * x - mu*mu - mu * m;
coef[1] = mu * x - mu * b + mu * z - m * b;
coef[0] = b * z;
TMath::RootsCubic(coef,roots[0],roots[1],roots[2]);
e = roots[1];
}
f =LikeMod6(mu,b,e,x,z,m);
}
return f;
}
Double_t TRolke::LikeMod6(Double_t mu,Double_t b,Double_t e,Int_t x,Int_t z,Int_t m)
{
Double_t f = 0.0;
if (z > 100 || m > 100) {
f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+(m*TMath::Log(m) - m)-(z*TMath::Log(z) - z) - ((m-z)*TMath::Log(m-z) - m + z)+z*TMath::Log(e)+(m-z)*TMath::Log(1-e));
} else {
f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+TMath::Log(TMath::Factorial(m))-TMath::Log(TMath::Factorial(z))-TMath::Log(TMath::Factorial(m-z))+z*TMath::Log(e)+(m-z)*TMath::Log(1-e));
}
return f;
}
Double_t TRolke::EvalLikeMod7(Double_t mu, Int_t x, Double_t em, Double_t e, Double_t sde, Double_t b, Int_t what)
{
Double_t v=sde*sde;
Double_t f = 0.;
if(what == 1) {
f = (x-b)/em;
}
if(what == 2) {
mu = (x-b)/em;
e = em;
f = LikeMod7(mu, b, e, x, em, v);
}
if(what == 3) {
if(mu==0) {
e = em;
} else {
e = ( -(mu*em-b-mu*mu*v)-TMath::Sqrt((mu*em-b-mu*mu*v)*(mu*em-b-mu*mu*v)+4*mu*(x*mu*v-mu*b*v + b * em)))/( - mu)/2;
}
f = LikeMod7(mu, b, e, x, em, v);
}
return f;
}
Double_t TRolke::LikeMod7(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t em,Double_t v)
{
return 2*(x*TMath::Log(e*mu+b)-(e*mu + b)-LogFactorial(x)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
}
Double_t TRolke::EvalPolynomial(Double_t x, const Int_t coef[], Int_t N)
{
const Int_t *p;
p = coef;
Double_t ans = *p++;
Int_t i = N;
do
ans = ans * x + *p++;
while( --i );
return ans;
}
Double_t TRolke::EvalMonomial(Double_t x, const Int_t coef[], Int_t N)
{
Double_t ans;
const Int_t *p;
p = coef;
ans = x + *p++;
Int_t i = N-1;
do
ans = ans * x + *p++;
while( --i );
return ans;
}
Double_t TRolke::LogFactorial(Int_t n)
{
if (n>69)
return n*TMath::Log(n)-n + 0.5*TMath::Log(TMath::TwoPi()*n);
else
return TMath::Log(TMath::Factorial(n));
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.