ChisqConsistency.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//  $Id: ChisqConsistency.cxx,v 1.1.1.1 2005/04/21 01:17:17 zhangy Exp $
//
// Description:
//
// Environment:
//  Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
//  Bob Jacobsen, Ed Iskander
//
// Copyright Information:
//  Copyright (C) 1996  Lawrence Berkeley Laboratory
//
//------------------------------------------------------------------------
 
//-----------------------
// This Class's Header --
//-----------------------
// #include "BaBar/BaBar.hh"
#include "ProbTools/ChisqConsistency.h"
 
#include <float.h>
#include <iostream>
#include <math.h>
 
#include "ProbTools/NumRecipes.h"
 
// #include "ErrLogger/ErrLog.hh"
 
// prototype the cernlib function
extern "C" {
float  chisin_( const float&, const int& );
double dgausn_( double& arg );
}
 
ChisqConsistency::ChisqConsistency() : _chisq( -1.0 ), _nDof( 0 ) {}
 
ChisqConsistency::ChisqConsistency( double chisq, double nDof )
    : _chisq( chisq ), _nDof( nDof ) {
  double z2 = 0.5 * _chisq;
  double n2 = 0.5 * _nDof;
 
  if ( n2 <= 0 || z2 < 0 )
  {
    std::cout << "ErrMsg(warning)"
              << " Got unphysical values: chisq = " << chisq << " #dof = " << nDof
              << std::endl;
    _value      = 0;
    _likelihood = 0;
    setStatus( Consistency::unPhysical );
    return;
  }
  setStatus( Consistency::OK );
 
  // given that n2>0 && z2>=0, gammq will NOT abort
  _value = NumRecipes::gammq( n2, z2 );
 
  if ( _chisq == 0 ) { _likelihood = 1; }
  else
  {
    double loglike = ( n2 - 1 ) * log( z2 ) - z2 - NumRecipes::gammln( n2 );
    if ( loglike < DBL_MIN_EXP )
    {
      _likelihood = 0;
      setStatus( Consistency::underFlow );
    }
    else { _likelihood = 0.5 * exp( loglike ); }
  }
}
 
ChisqConsistency::ChisqConsistency( unsigned nDof, double prob ) : _nDof( nDof ) {
  if ( prob >= 0.0 || prob <= 1.0 || nDof < 0 ) _value = prob;
  else
  {
    std::cout << "ErrMsg(warning)"
              << " Got unphysical values: prob = " << prob << " #dof = " << nDof << std::endl;
    _value      = 0;
    _likelihood = 0;
    setStatus( Consistency::unPhysical );
    return;
  }
  setStatus( Consistency::OK );
  if ( prob != 1.0 )
  {
    // use the cernlib function to get chisq.  Note the funny convention on prob!!
    float value = 1.0 - float( _value );
    int   ndof  = nDof;
    if ( value < 1.0 ) _chisq = chisin_( value, ndof );
    else _chisq = log( double( FLT_MAX ) );
    // use the same algorithm as above to get loglikelihood
    double z2 = 0.5 * _chisq;
    double n2 = 0.5 * _nDof;
    if ( _chisq == 0 ) { _likelihood = 1; }
    else
    {
      double loglike = ( n2 - 1 ) * log( z2 ) - z2 - NumRecipes::gammln( n2 );
      if ( loglike < DBL_MIN_EXP )
      {
        _likelihood = 0;
        setStatus( Consistency::underFlow );
      }
      else { _likelihood = 0.5 * exp( loglike ); }
    }
  }
}
 
ChisqConsistency::ChisqConsistency( const ChisqConsistency& other )
    : Consistency( other ), _chisq( other._chisq ), _nDof( other._nDof ) {}
 
ChisqConsistency& ChisqConsistency::operator=( const ChisqConsistency& other ) {
  if ( this != &other )
  {
    Consistency::operator=( other );
    _chisq = other._chisq;
    _nDof  = other._nDof;
  }
  return *this;
}