Double_t bi=1;
Double_t result=0;
Int_t m=2*((n+int(sqrt(float(iacc*n))))); // Downward recurrence from even m
- for (Int_t j=m; j<=1; j--)
+ for (Int_t j=m; j>=1; j--)
{
bim=bip+double(j)*tox*bi;
bip=bi;
//
// Note : <T>=0 Var(T)=ndf/(ndf-2)
- Double_t pi=acos(-1.);
-
TF1 tdist("tdist","(TMath::Gamma(([0]+1.)/2.)/(sqrt(pi*[0])*TMath::Gamma([0]/2.)))*pow(1.+x*x/[0],-([0]+1.)/2.)");
TString title="Student's t distribution for ndf = ";
title+=ndf;
return val;
}
///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::Nfac(Int_t n,Int_t mode) const
+{
+// Compute n!.
+// The algorithm can be selected by the "mode" input argument.
+//
+// mode : 0 ==> Calculation by means of straightforward multiplication
+// : 1 ==> Calculation by means of Stirling's approximation
+// : 2 ==> Calculation by means of n!=Gamma(n+1)
+//
+// For large n the calculation modes 1 and 2 will in general be faster.
+// By default mode=0 is used.
+// For n<0 the value 0 will be returned.
+//
+// Note : Because of Double_t value overflow the maximum value is n=170.
+//
+//--- NvE 20-jan-2007 Utrecht University
+
+ if (n<0) return 0;
+ if (n==0) return 1;
+
+ Double_t twopi=2.*acos(-1.);
+ Double_t z=0;
+ Double_t nfac=1;
+ Int_t i=n;
+
+ switch (mode)
+ {
+ case 0: // Straightforward multiplication
+ while (i>1)
+ {
+ nfac*=Double_t(i);
+ i--;
+ }
+ break;
+
+ case 1: // Stirling's approximation
+ z=n;
+ nfac=sqrt(twopi)*pow(z,z+0.5)*exp(-z)*(1.+1./(12.*z));
+ break;
+
+ case 2: // Use of Gamma(n+1)
+ z=n+1;
+ nfac=Gamma(z);
+ break;
+
+ default:
+ nfac=0;
+ break;
+ }
+
+ return nfac;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::LnNfac(Int_t n,Int_t mode) const
+{
+// Compute ln(n!).
+// The algorithm can be selected by the "mode" input argument.
+//
+// mode : 0 ==> Calculation via evaluation of n! followed by taking ln(n!)
+// : 1 ==> Calculation via Stirling's approximation ln(n!)=0.5*ln(2*pi)+(n+0.5)*ln(n)-n+1/(12*n)
+// : 2 ==> Calculation by means of ln(n!)=LnGamma(n+1)
+//
+// Note : Because of Double_t value overflow the maximum value is n=170 for mode=0.
+//
+// For mode=2 rather accurate results are obtained for both small and large n.
+// By default mode=2 is used.
+// For n<1 the value 0 will be returned.
+//
+//--- NvE 20-jan-2007 Utrecht University
+
+ if (n<=1) return 0;
+
+ Double_t twopi=2.*acos(-1.);
+ Double_t z=0;
+ Double_t lognfac=0;
+
+ switch (mode)
+ {
+ case 0: // Straightforward ln(n!)
+ z=Nfac(n);
+ lognfac=log(z);
+ break;
+
+ case 1: // Stirling's approximation
+ z=n;
+ lognfac=0.5*log(twopi)+(z+0.5)*log(z)-z+1./(12.*z);
+ break;
+
+ case 2: // Use of LnGamma(n+1)
+ z=n+1;
+ lognfac=LnGamma(z);
+ break;
+
+ default:
+ lognfac=0;
+ break;
+ }
+
+ return lognfac;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::LogNfac(Int_t n,Int_t mode) const
+{
+// Compute log_10(n!).
+// First ln(n!) is evaluated via invokation of LnNfac(n,mode).
+// Then the algorithm log_10(z)=ln(z)*log_10(e) is used.
+//
+//--- NvE 20-jan-2007 Utrecht University
+
+ Double_t e=exp(1.);
+
+ Double_t val=LnNfac(n,mode);
+ val*=log10(e);
+
+ return val;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::PsiValue(Int_t m,Int_t* n,Double_t* p,Int_t f) const
+{
+// Provide the Bayesian Psi value of observations of a counting experiment
+// w.r.t. a Bernoulli class hypothesis B_m.
+// The hypothesis B_m represents a counting experiment with m different
+// possible outcomes and is completely defined by the probabilities
+// of the various outcomes (and the requirement that the sum of all these
+// probabilities equals 1).
+//
+// The Psi value provides (in dB scale) the amount of support that the
+// data can maximally give to any Bernoulli class hypothesis different
+// from the currently specified B_m.
+//
+// To be specific : Psi=-10*log[p(D|B_m I)]
+//
+// where p(D|B_m I) represents the likelihood of the data D under the condition
+// that B_m and some prior I are true.
+//
+// A Psi value of zero indicates a perfect match between the observations
+// and the specified hypothesis.
+// Further mathematical details can be found in astro-ph/0702029.
+//
+// m : The number of different possible outcomes of the counting experiment
+// n : The observed number of different outcomes
+// p : The probabilities of the different outcomes according to the hypothesis
+// f : Flag to indicate the use of a frequentist (Stirling) approximation (f=1)
+// or the exact Bayesian expression (f=0).
+//
+// In case no probabilities are given (i.e. pk=0), a uniform distribution
+// is assumed.
+//
+// The default values are pk=0 and f=0.
+//
+// In the case of inconsistent input, a Psi value of -1 is returned.
+//
+//--- NvE 03-oct-2007 Utrecht University
+
+ Double_t psi=-1;
+
+ if (m<=0 || !n) return psi;
+
+ Int_t ntot=0;
+ for (Int_t j=0; j<m; j++)
+ {
+ if (n[j]>0) ntot+=n[j];
+ }
+
+ psi=0;
+ Double_t pk=1./float(m); // Prob. of getting outcome k for a uniform distr.
+ for (Int_t i=0; i<m; i++)
+ {
+ if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
+ if (n[i]>0 && pk>0)
+ {
+ if (!f) // Exact Bayesian expression
+ {
+ psi+=double(n[i])*log10(pk)-LogNfac(n[i]);
+ }
+ else // Frequentist (Stirling) approximation
+ {
+ if (ntot>0) psi+=double(n[i])*log10(n[i]/(ntot*pk));
+ }
+ }
+ }
+
+ if (!f) // Exact Bayesian expression
+ {
+ psi+=LogNfac(ntot);
+ psi*=-10.;
+ }
+ else // Frequentist (Stirling) approximation
+ {
+ psi*=10.;
+ }
+
+ return psi;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::PsiValue(TH1F* his,TH1F* hyp,TF1* pdf,Int_t f) const
+{
+// Provide the Bayesian Psi value of observations of a counting experiment
+// (in histogram format) w.r.t. a Bernoulli class hypothesis B_m.
+// The hypothesis B_m represents a counting experiment with m different
+// possible outcomes and is completely defined by the probabilities
+// of the various outcomes (and the requirement that the sum of all these
+// probabilities equals 1).
+// The specification of a hypothesis B_m can be provided either in
+// histogram format (hyp) or via a probability distribution function (pdf),
+// as outlined below.
+// Note : The pdf does not need to be normalised.
+//
+// The Psi value provides (in dB scale) the amount of support that the
+// data can maximally give to any Bernoulli class hypothesis different
+// from the currently specified B_m.
+//
+// To be specific : Psi=-10*log[p(D|B_m I)]
+//
+// where p(D|B_m I) represents the likelihood of the data D under the condition
+// that B_m and some prior I are true.
+//
+// A Psi value of zero indicates a perfect match between the observations
+// and the specified hypothesis.
+// Further mathematical details can be found in astro-ph/0702029.
+//
+// his : The experimental observations in histogram format
+// hyp : Hypothetical observations according to some hypothesis
+// pdf : Probability distribution function for the hypothesis
+// f : Flag to indicate the use of a frequentist (Stirling) approximation (f=1)
+// or the exact Bayesian expression (f=0).
+//
+// In case no hypothesis is specified (i.e. hyp=0 and pdf=0), a uniform
+// background distribution is assumed.
+//
+// Default values are : hyp=0, pdf=0 and f=0.
+//
+// In the case of inconsistent input, a Psi value of -1 is returned.
+//
+//--- NvE 03-oct-2007 Utrecht University
+
+ Double_t psi=-1;
+
+ if (!his) return psi;
+
+ TAxis* xaxis=his->GetXaxis();
+ Double_t xmin=xaxis->GetXmin();
+ Double_t xmax=xaxis->GetXmax();
+ Double_t range=xmax-xmin;
+ Int_t nbins=his->GetNbinsX();
+ Double_t nensig=his->GetEntries();
+
+ if (nbins<=0 || nensig<=0 || range<=0) return psi;
+
+ Int_t* n=new Int_t[nbins];
+ Double_t* p=new Double_t[nbins];
+ Int_t nk;
+ Double_t pk;
+
+ // Uniform hypothesis distribution
+ if (!hyp && !pdf)
+ {
+ for (Int_t i=1; i<=nbins; i++)
+ {
+ nk=int(his->GetBinContent(i));
+ pk=his->GetBinWidth(i)/range;
+ n[i-1]=0;
+ p[i-1]=0;
+ if (nk>0) n[i-1]=nk;
+ if (pk>0) p[i-1]=pk;
+ }
+ psi=PsiValue(nbins,n,p,f);
+ }
+
+ // Hypothesis specified via a pdf
+ if (pdf && !hyp)
+ {
+ pdf->SetRange(xmin,xmax);
+ Double_t ftot=pdf->Integral(xmin,xmax);
+ if (ftot>0)
+ {
+ Double_t x1,x2;
+ for (Int_t ipdf=1; ipdf<=nbins; ipdf++)
+ {
+ nk=int(his->GetBinContent(ipdf));
+ x1=his->GetBinLowEdge(ipdf);
+ x2=x1+his->GetBinWidth(ipdf);
+ pk=pdf->Integral(x1,x2)/ftot;
+ n[ipdf-1]=0;
+ p[ipdf-1]=0;
+ if (nk>0) n[ipdf-1]=nk;
+ if (pk>0) p[ipdf-1]=pk;
+ }
+ psi=PsiValue(nbins,n,p,f);
+ }
+ }
+
+ // Hypothesis specified via a histogram
+ if (hyp && !pdf)
+ {
+ TH1F* href=new TH1F(*his);
+ href->Reset();
+ Double_t nenhyp=hyp->GetEntries();
+ Float_t x,y;
+ for (Int_t ihyp=1; ihyp<=hyp->GetNbinsX(); ihyp++)
+ {
+ x=hyp->GetBinCenter(ihyp);
+ y=hyp->GetBinContent(ihyp);
+ href->Fill(x,y);
+ }
+ for (Int_t j=1; j<=nbins; j++)
+ {
+ nk=int(his->GetBinContent(j));
+ pk=href->GetBinContent(j)/nenhyp;
+ n[j-1]=0;
+ p[j-1]=0;
+ if (nk>0) n[j-1]=nk;
+ if (pk>0) p[j-1]=pk;
+ }
+ psi=PsiValue(nbins,n,p,f);
+ delete href;
+ }
+
+ delete [] n;
+ delete [] p;
+
+ return psi;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::Chi2Value(Int_t m,Int_t* n,Double_t* p,Int_t* ndf) const
+{
+// Provide the frequentist chi-squared value of observations of a counting
+// experiment w.r.t. a Bernoulli class hypothesis B_m.
+// The hypothesis B_m represents a counting experiment with m different
+// possible outcomes and is completely defined by the probabilities
+// of the various outcomes (and the requirement that the sum of all these
+// probabilities equals 1).
+//
+// Further mathematical details can be found in astro-ph/0702029.
+//
+// m : The number of different possible outcomes of the counting experiment
+// n : The observed number of different outcomes
+// p : The probabilities of the different outcomes according to the hypothesis
+// ndf : The returned number of degrees of freedom
+//
+// Note : Both the arrays "n" and (when provided) "p" should be of dimension "m".
+//
+// In case no probabilities are given (i.e. pk=0), a uniform distribution
+// is assumed.
+//
+// The default values are pk=0 and ndf=0.
+//
+// In the case of inconsistent input, a chi-squared and ndf value of -1 is returned.
+//
+//--- NvE 03-oct-2007 Utrecht University
+
+ Double_t chi=-1;
+ if (ndf) *ndf=-1;
+
+ if (m<=0 || !n) return chi;
+
+ Int_t ntot=0;
+ for (Int_t j=0; j<m; j++)
+ {
+ if (n[j]>0) ntot+=n[j];
+ }
+
+ chi=0;
+ Double_t pk=1./float(m); // Prob. of getting outcome k for a uniform distr.
+ for (Int_t i=0; i<m; i++)
+ {
+ if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
+ if (n[i]>0 && pk>0 && ntot>0)
+ {
+ chi+=pow(double(n[i])-double(ntot)*pk,2)/(ntot*pk);
+ if (ndf) (*ndf)=(*ndf)+1;
+ }
+ }
+
+ return chi;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::Chi2Value(TH1F* his,TH1F* hyp,TF1* pdf,Int_t* ndf) const
+{
+// Provide the frequentist chi-squared value of observations of a counting
+// experiment (in histogram format) w.r.t. a Bernoulli class hypothesis B_m.
+// The hypothesis B_m represents a counting experiment with m different
+// possible outcomes and is completely defined by the probabilities
+// of the various outcomes (and the requirement that the sum of all these
+// probabilities equals 1).
+// The specification of a hypothesis B_m can be provided either in
+// histogram format (hyp) or via a probability distribution function (pdf),
+// as outlined below.
+// Note : The pdf does not need to be normalised.
+//
+// Further mathematical details can be found in astro-ph/0702029.
+//
+// his : The experimental observations in histogram format
+// hyp : Hypothetical observations according to some hypothesis
+// pdf : Probability distribution function for the hypothesis
+// ndf : The returned number of degrees of freedom
+//
+// In case no hypothesis is specified (i.e. hyp=0 and pdf=0), a uniform
+// background distribution is assumed.
+//
+// Default values are : hyp=0, pdf=0 and ndf=0.
+//
+// In the case of inconsistent input, a chi-squared and ndf value of -1 is returned.
+//
+//--- NvE 03-oct-2007 Utrecht University
+
+ Double_t chi=-1;
+
+ if (!his) return chi;
+
+ TAxis* xaxis=his->GetXaxis();
+ Double_t xmin=xaxis->GetXmin();
+ Double_t xmax=xaxis->GetXmax();
+ Double_t range=xmax-xmin;
+ Int_t nbins=his->GetNbinsX();
+ Double_t nensig=his->GetEntries();
+
+ if (nbins<=0 || nensig<=0 || range<=0) return chi;
+
+ Int_t* n=new Int_t[nbins];
+ Double_t* p=new Double_t[nbins];
+ Int_t nk;
+ Double_t pk;
+
+ // Uniform hypothesis distribution
+ if (!hyp && !pdf)
+ {
+ for (Int_t i=1; i<=nbins; i++)
+ {
+ nk=int(his->GetBinContent(i));
+ pk=his->GetBinWidth(i)/range;
+ n[i-1]=0;
+ p[i-1]=0;
+ if (nk>0) n[i-1]=nk;
+ if (pk>0) p[i-1]=pk;
+ }
+ chi=Chi2Value(nbins,n,p,ndf);
+ }
+
+ // Hypothesis specified via a pdf
+ if (pdf && !hyp)
+ {
+ pdf->SetRange(xmin,xmax);
+ Double_t ftot=pdf->Integral(xmin,xmax);
+ if (ftot>0)
+ {
+ Double_t x1,x2;
+ for (Int_t ipdf=1; ipdf<=nbins; ipdf++)
+ {
+ nk=int(his->GetBinContent(ipdf));
+ x1=his->GetBinLowEdge(ipdf);
+ x2=x1+his->GetBinWidth(ipdf);
+ pk=pdf->Integral(x1,x2)/ftot;
+ n[ipdf-1]=0;
+ p[ipdf-1]=0;
+ if (nk>0) n[ipdf-1]=nk;
+ if (pk>0) p[ipdf-1]=pk;
+ }
+ chi=Chi2Value(nbins,n,p,ndf);
+ }
+ }
+
+ // Hypothesis specified via a histogram
+ if (hyp && !pdf)
+ {
+ TH1F* href=new TH1F(*his);
+ href->Reset();
+ Double_t nenhyp=hyp->GetEntries();
+ Float_t x,y;
+ for (Int_t ihyp=1; ihyp<=hyp->GetNbinsX(); ihyp++)
+ {
+ x=hyp->GetBinCenter(ihyp);
+ y=hyp->GetBinContent(ihyp);
+ href->Fill(x,y);
+ }
+ for (Int_t j=1; j<=nbins; j++)
+ {
+ nk=int(his->GetBinContent(j));
+ pk=href->GetBinContent(j)/nenhyp;
+ n[j-1]=0;
+ p[j-1]=0;
+ if (nk>0) n[j-1]=nk;
+ if (pk>0) p[j-1]=pk;
+ }
+ chi=Chi2Value(nbins,n,p,ndf);
+ delete href;
+ }
+
+ delete [] n;
+ delete [] p;
+
+ return chi;
+}
+///////////////////////////////////////////////////////////////////////////
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff