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

/* ======================================================================== *\
!
! *
! * 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));
}