source: trunk/MagicSoft/Mars/mcalib/MHGausEvent.cc@ 3152

/* ======================================================================== *\
!
! *
! * 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): Markus Gaug 11/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MHGausEvent
//
//  A base class for all kind of event which follow a Gaussian distribution 
//  with time, i.e. observables containing white noise.
//
//  The class provides the basic tools for fitting, 
//  spectrum analysis, etc. 
//
//////////////////////////////////////////////////////////////////////////////
#include "MHGausEvent.h"

#include <TH1.h>
#include <TF1.h>

#include "MFFT.h"
#include "MArray.h"

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MHGausEvent);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHGausEvent::MHGausEvent(const char *name, const char *title)
    : fHGausHist(NULL), fHPowerProbability(NULL), 
      fFGausFit(NULL), fFExpFit(NULL),
      fEvents(NULL), fPowerSpectrum(NULL)
{ 

    fName  = name  ? name  : "MHGausEvent";
    fTitle = title ? title : "Events which follow a Gaussian distribution";
    
    Clear();
}


MHGausEvent::~MHGausEvent()
{
  Clear();

  if (fHGausHist)
    delete fHGausHist;
}
      


void MHGausEvent::Clear(Option_t *o)
{
  
  fGausHistBins      = 50;
  fGausHistAxisFirst = 0.;
  fGausHistAxisLast  = 100.;

  fPowerProbabilityBins = 30;

  fProbLimit         = 0.01;
  fGausFitOK         = kFALSE;
  fExpFitOK          = kFALSE;
  fOscillating       = kFALSE;

  fMean              = 0.;
  fSigma             = 0.;
  fMeanErr           = 0.;
  fSigmaErr          = 0.;

  fProb              = 0.;

  if (fHPowerProbability)
    delete fHPowerProbability;
  if (fFGausFit)
    delete fFGausFit; 
  if (fEvents)
    delete fEvents;
  if (fPowerSpectrum)  
    delete fPowerSpectrum;     
}


void MHGausEvent::Reset()
{
  
  Clear();
  fHGausHist->Reset();

}

const Double_t MHGausEvent::GetChiSquare()  const 
{
  return ( fFGausFit ? fFGausFit->GetChisquare() : 0.);
}

const Int_t MHGausEvent::GetNdf() const 
{
  return ( fFGausFit ? fFGausFit->GetNDF() : 0);
}


const Double_t MHGausEvent::GetExpChiSquare()  const 
{
  return ( fFExpFit ? fFExpFit->GetChisquare() : 0.);
}


const Int_t MHGausEvent::GetExpNdf()  const 
{
  return ( fFExpFit ? fFExpFit->GetNDF() : 0);
}

const Double_t MHGausEvent::GetExpProb()  const 
{
  return ( fFExpFit ? fFExpFit->GetProb() : 0.);
}

const Double_t MHGausEvent::GetOffset()  const 
{
  return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.);
}

const Double_t MHGausEvent::GetSlope()  const 
{
  return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.);
}



Bool_t MHGausEvent::CheckOscillations()
{

  if (fFExpFit)
    return IsOscillating();

  if (!fEvents)
    return kFALSE;

  //
  // The number of entries HAS to be a potence of 2, 
  // so we can only cut out from the last potence of 2 to the rest. 
  // Another possibility would be to fill everything with 
  // zeros, but that gives a low frequency peak, which we would 
  // have to cut out later again. 
  //
  // So, we have to live with the possibility that at the end 
  // of the calibration run, something has happened without noticing 
  // it...
  //
  
  // This cuts only the non-used zero's, but MFFT will later cut the rest
  MArray::StripZeros(*fEvents);

  MFFT fourier;

  fPowerSpectrum     = fourier.PowerSpectrumDensity(fEvents);
  fHPowerProbability = ProjectArray(*fPowerSpectrum, fPowerProbabilityBins,
                                    "PowerProbability",
                                    "Probability of Power occurrance");
  //
  // First guesses for the fit (should be as close to reality as possible, 
  //
  const Double_t xmax = fHPowerProbability->GetXaxis()->GetXmax();

  fFExpFit = new TF1("FExpFit","exp([0]-[1]*x)",0.,xmax);

  const Double_t slope_guess  = (TMath::Log(fHPowerProbability->GetEntries())+1.)/xmax;
  const Double_t offset_guess = slope_guess*xmax;

  fFExpFit->SetParameters(offset_guess, slope_guess);
  fFExpFit->SetParNames("Offset","Slope");
  fFExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess);
  fFExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess);
  fFExpFit->SetRange(0.,xmax);

  fHPowerProbability->Fit(fFExpFit,"RQL0");
  
  if (GetExpProb() < fProbLimit)
    fExpFitOK = kFALSE;
  
  // For the moment, this is the only check, later we can add more...
  fOscillating = fExpFitOK;

  return fOscillating;
}


Bool_t MHGausEvent::IsEmpty() const
{
    return !(fHGausHist->GetEntries());
}

Bool_t MHGausEvent::IsOscillating()
{

  if (fFExpFit)
    return fOscillating;

  return CheckOscillations();

}


Bool_t MHGausEvent::FitGaus(Option_t *option)
{

  if (IsGausFitOK())
    return kTRUE;

  //
  // First, cut the edges which contain only zeros and rebin 
  // to about 20 bins. 
  //
  // (ATTENTION: The Chisquare method is more sensitive, 
  // the _less_ bins, you have!)
  //
  Int_t newbins = 20;
  MH::StripZeros(fHGausHist,newbins);
  
  //
  // Get the fitting ranges
  //
  Axis_t rmin = fHGausHist->GetXaxis()->GetFirst();
  Axis_t rmax = fHGausHist->GetXaxis()->GetLast(); 

  //
  // First guesses for the fit (should be as close to reality as possible, 
  //
  const Stat_t   entries     = fHGausHist->Integral("width");
  const Double_t mu_guess    = fHGausHist->GetBinCenter(fHGausHist->GetMaximumBin());
  const Double_t sigma_guess = (rmax-rmin)/2.;
  const Double_t area_guess  = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess;

  fFGausFit = new TF1("GausFit","gaus",rmin,rmax);

  if (!fFGausFit) 
    {
    *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl;
    return kFALSE;
    }
  
  fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
  fFGausFit->SetParNames("Area","#mu","#sigma");
  fFGausFit->SetParLimits(0,0.,entries);
  fFGausFit->SetParLimits(1,rmin,rmax);
  fFGausFit->SetParLimits(2,0.,rmax-rmin);
  fFGausFit->SetRange(rmin,rmax);

  fHGausHist->Fit(fFGausFit,option);
  
  fMean     = fFGausFit->GetParameter(1);
  fSigma    = fFGausFit->GetParameter(2);
  fMeanErr  = fFGausFit->GetParError(1);
  fSigmaErr = fFGausFit->GetParError(2);

  fProb     = fFGausFit->GetProb();
  //
  // The fit result is accepted under condition:
  // 1) The results are not nan's
  // 2) The NDF is not smaller than fNDFLimit (5)
  // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
  //
  if (    TMath::IsNaN(fMean) 
      || TMath::IsNaN(fMeanErr)
      || TMath::IsNaN(fProb)    
      || TMath::IsNaN(fSigma)
      || TMath::IsNaN(fSigmaErr) )
    {
      fGausFitOK = kFALSE;
      return kFALSE;
    }
  
  fGausFitOK = kTRUE;
  return kTRUE;
}

void MHGausEvent::Print(const Option_t *o) const 
{
  
  *fLog << all                                                        << endl;
  *fLog << all << "Results of the Gauss Fit: "                        << endl;
  *fLog << all << "Mean: "             << GetMean()                   << endl;
  *fLog << all << "Sigma: "            << GetSigma()                  << endl;
  *fLog << all << "Chisquare: "        << GetChiSquare()              << endl;
  *fLog << all << "DoF: "              << GetNdf()                    << endl;
  *fLog << all << "Probability: "      << GetProb()                   << endl;
  *fLog << all                                                        << endl;
  
}