source: trunk/MagicSoft/Mars/mbase/MMath.cc@ 4966

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz  3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MMath
//
/////////////////////////////////////////////////////////////////////////////
#include "MMath.h"

ClassImp(MMath);

using namespace std;

// --------------------------------------------------------------------------
//
// Calculate Significance as
// significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
//
// s: total number of events in signal region
// b: number of background events in signal region
// 
Double_t MMath::Significance(Double_t s, Double_t b)
{
    const Double_t k = b==0 ? 0 : s/b;
    const Double_t f = s+k*k*b;

    return f==0 ? 0 : (s-b)/Sqrt(f);
}

// --------------------------------------------------------------------------
//
// Symmetrized significance - this is somehow analog to
// SignificanceLiMaSigned
//
// Returns Significance(s,b) if s>b otherwise -Significance(b, s);
// 
Double_t MMath::SignificanceSym(Double_t s, Double_t b)
{
    return s>b ? Significance(s, b) : -Significance(b, s);
}

// --------------------------------------------------------------------------
//
//  calculates the significance according to Li & Ma
//  ApJ 272 (1983) 317, Formula 17
//
//  s                    // s: number of on events
//  b                    // b: number of off events
//  alpha = t_on/t_off;  // t: observation time
//
//  The significance has the same (positive!) value for s>b and b>s.
//
//  Returns -1 if sum<0 or alpha<0 or the argument of sqrt<0
//  Returns  0 if s+b==0
//
Double_t MMath::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha)
{
    const Double_t sum = s+b;

    if (sum==0)
        return 0;

    if (sum<0 || alpha<=0)
        return -1;

    const Double_t l = s*Log(s/sum*(alpha+1)/alpha);
    const Double_t m = b*Log(b/sum*(alpha+1)      );

    return l+m<0 ? -1 : Sqrt((l+m)*2);
}

// --------------------------------------------------------------------------
//
// Calculates MMath::SignificanceLiMa(s, b, alpha). Returns 0 if the
// calculation has failed. Otherwise the Li/Ma significance which was
// calculated. If s<b a negative value is returned.
//
Double_t MMath::SignificanceLiMaSigned(Double_t s, Double_t b, Double_t alpha)
{
    const Double_t sig = SignificanceLiMa(s, b, alpha);
    if (sig<=0)
        return 0;

    return Sign(sig, s-alpha*b);
}

// --------------------------------------------------------------------------
//
// Returns: 2/(sigma*sqrt(2))*integral[0,x](exp(-(x-mu)^2/(2*sigma^2)))
//
Double_t MMath::GaussProb(Double_t x, Double_t sigma, Double_t mean)
{
    static const Double_t sqrt2 = Sqrt(2.);
    return Erf((x-mean)/(sigma*sqrt2));
}

// -------------------------------------------------------------------------
//
// Quadratic interpolation
//
// calculate the parameters of a parabula such that
//    y(i) = a + b*x(i) + c*x(i)^2
//
// If the determinant==0 an empty TVector3 is returned.
//
TVector3 MMath::GetParab(const TVector3 &x, const TVector3 &y)
{
    const TVector3 sq(x(0)*x(0), x(1)*x(1), x(2)*x(2));

    const TVector3 ai2 = sq.Cross(sq);

    const Double_t det = x.Dot(ai2);
    if (det==0)
        return TVector3();

    const TVector3 ai1 = x.Cross(sq);
    const TVector3 ai3 = x.Cross(x);

    TVector3 res(y.Dot(ai1), y.Dot(ai2), y.Dot(ai3));
    res *= 1./det;

    return res;
}

Double_t MMath::InterpolParabLin(const TVector3 &vx, const TVector3 &vy, Double_t x)
{
    const TVector3 c = GetParab(vx, vy);
    return c(0) + c(1)*x + c(2)*x*x;
}

Double_t MMath::InterpolParabLog(const TVector3 &vx, const TVector3 &vy, Double_t x)
{
    const Double_t l0 = TMath::Log10(vx(0));
    const Double_t l1 = TMath::Log10(vx(1));
    const Double_t l2 = TMath::Log10(vx(2));

    const TVector3 vx0(l0, l1, l2);
    return pow(10, InterpolParabLin(vx0, vy, TMath::Log10(x)));
}

Double_t MMath::InterpolParabCos(const TVector3 &vx, const TVector3 &vy, Double_t x)
{
    const Double_t l0 = TMath::Cos(vx(0));
    const Double_t l1 = TMath::Cos(vx(1));
    const Double_t l2 = TMath::Cos(vx(2));

    const TVector3 vx0(l0, l1, l2);
    return TMath::ACos(InterpolParabLin(vx0, vy, TMath::Cos(x)));
}