source: trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc@ 3625

/* ======================================================================== *\
!
! *
! * 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 02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MHGausEvents
//
//  A base class for events which are believed to follow a Gaussian distribution 
//  with time, e.g. calibration events, observables containing white noise, ...
//
//  MHGausEvents derives from MH, thus all features of MH can be used by a class 
//  deriving from MHGausEvents, especially the filling functions. 
//
//  The central objects are: 
//
//  1) The TH1F fHGausHist: 
//     ====================
//   
//     It is created with a default name and title and resides in directory NULL.
//     - Any class deriving from MHGausEvents needs to apply a binning to fHGausHist
//       (e.g. with the function TH1F::SetBins(..) )
//     - The histogram is filled with the functions FillHist() or FillHistAndArray(). 
//     - The histogram can be fitted with the function FitGaus(). This involves stripping 
//       of all zeros at the lower and upper end of the histogram and re-binning to 
//       a new number of bins, specified in the variable fBinsAfterStripping.      
//     - The fit result's probability is compared to a reference probability fProbLimit
//       The NDF is compared to fNDFLimit and a check is made whether results are NaNs. 
//       Accordingly, a flag IsGausFitOK() is set.
// 
//  2) The TArrayF fEvents:
//     ==========================
// 
//     It is created with 0 entries and not expanded unless FillArray() or FillHistAndArray()
//     are called.
//     - A first call to FillArray() or FillHistAndArray() initializes fEvents by default 
//       to 512 entries. 
//     - Any later call to FillArray() or FillHistAndArray() fills up the array. 
//       Reaching the limit, the array is expanded by a factor 2.
//     - The array can be fourier-transformed into the array fPowerSpectrum. 
//       Note that any FFT accepts only number of events which are a power of 2. 
//       Thus, fEvents is cut to the next power of 2 smaller than its actual number of entries. 
//       Be aware that you might lose information at the end of your analysis. 
//     - Calling the function CreateFourierSpectrum() creates the array fPowerSpectrum 
//       and its projection fHPowerProbability which in turn is fit to an exponential. 
//     - The fit result's probability is compared to a referenc probability fProbLimit 
//       and accordingly the flag IsExpFitOK() is set.
//     - The flag IsFourierSpectrumOK() is set accordingly to IsExpFitOK(). 
//       Later, a closer check will be installed. 
//     - You can display all arrays by calls to: CreateGraphEvents() and 
//       CreateGraphPowerSpectrum() and successive calls to the corresponding Getters.
//
//////////////////////////////////////////////////////////////////////////////
#include "MHGausEvents.h"

#include <TH1.h>
#include <TF1.h>
#include <TGraph.h>
#include <TPad.h>
#include <TVirtualPad.h>
#include <TCanvas.h>
#include <TStyle.h>

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

#include "MH.h"

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

ClassImp(MHGausEvents);

using namespace std;

const Float_t  MHGausEvents::fgProbLimit            = 0.001;
const Int_t    MHGausEvents::fgNDFLimit             = 2;
const Int_t    MHGausEvents::fgPowerProbabilityBins = 20;
const Int_t    MHGausEvents::fgBinsAfterStripping   = 40;
// --------------------------------------------------------------------------
//
// Default Constructor. 
// Sets: 
// - the default Probability Bins for fPowerProbabilityBins (fgPowerProbabilityBins)
// - the default Probability Limit for fProbLimit  (fgProbLimit)
// - the default NDF Limit for fNDFLimit  (fgNDFLimit)
// - the default number of bins after stripping for fBinsAfterStipping (fgBinsAfterStipping)
// - the default name of the fHGausHist   ("HGausHist")
// - the default title of the fHGausHist  ("Histogram of Events with Gaussian Distribution")
// - the default directory of the fHGausHist (NULL)
//
// Initializes:
// - fEvents to 0 entries
// - fHGausHist()
// - all other pointers to NULL
// - all variables to 0.
// - all flags to kFALSE
//
MHGausEvents::MHGausEvents(const char *name, const char *title)
    : fEventFrequency(0), fHPowerProbability(NULL), 
      fPowerSpectrum(NULL),
      fGraphEvents(NULL), fGraphPowerSpectrum(NULL),
      fHGausHist(), fEvents(0),
      fFGausFit(NULL), fFExpFit(NULL)
{ 

  fName  = name  ? name  : "MHGausEvents";
  fTitle = title ? title : "Events with expected Gaussian distributions";

  Clear();
  
  SetPowerProbabilityBins();
  SetProbLimit();
  SetNDFLimit();
  SetBinsAfterStripping();

  fHGausHist.SetName("HGausHist");
  fHGausHist.SetTitle("Histogram of Events with Gaussian Distribution");
  // important, other ROOT will not draw the axes:
  fHGausHist.UseCurrentStyle();
  fHGausHist.SetDirectory(NULL);
}



// --------------------------------------------------------------------------
//
// Default Destructor. 
//
// Deletes (if Pointer is not NULL):
// 
// - fHPowerProbability
// - fFGausFit 
// - fFExpFit
// - fPowerSpectrum     
// - fGraphEvents
// - fGraphPowerSpectrum
// 
MHGausEvents::~MHGausEvents()
{

  // delete histograms
  if (fHPowerProbability)
    delete fHPowerProbability;
  
  // delete fits
  if (fFGausFit)
    delete fFGausFit; 
  if (fFExpFit)
    delete fFExpFit;
  
  // delete arrays
  if (fPowerSpectrum)  
    delete fPowerSpectrum;     

  // delete graphs
  if (fGraphEvents)
    delete fGraphEvents;
  if (fGraphPowerSpectrum)
    delete fGraphPowerSpectrum;
}
      
// --------------------------------------------------------------------------
//
// Sets:
// - all other pointers to NULL
// - all variables to 0., except fEventFrequency
// - all flags to kFALSE
// 
// Deletes (if not NULL):
// - all pointers
//
void MHGausEvents::Clear(Option_t *o)
{

  SetGausFitOK        ( kFALSE );
  SetExpFitOK         ( kFALSE );
  SetFourierSpectrumOK( kFALSE );

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

  if (fHPowerProbability)
    {
      delete fHPowerProbability;
      fHPowerProbability = NULL;
    }
  
  // delete fits
  if (fFGausFit)
    {
      delete fFGausFit; 
      fFGausFit = NULL;
    }
  
  if (fFExpFit)
    {
      delete fFExpFit;
      fFExpFit = NULL;
    }
  
  // delete arrays
  if (fPowerSpectrum)  
    {
      delete fPowerSpectrum;     
      fPowerSpectrum = NULL;
    }

  // delete graphs
  if (fGraphEvents)
    {
      delete fGraphEvents;
      fGraphEvents = NULL;
    }

  if (fGraphPowerSpectrum)
    {
      delete fGraphPowerSpectrum;
      fGraphPowerSpectrum = NULL;
    }
}

// --------------------------------------------------------------------------
//
// Executes:
// - Clear()
// - fHGausHist.Reset()
// - fEvents.Set(0)
//
void MHGausEvents::Reset()
{

  Clear();
  fHGausHist.Reset();
  fEvents.Set(0);

}

// --------------------------------------------------------------------------
//
// Executes:
// - FillArray()
// - FillHist()
//
Bool_t MHGausEvents::FillHistAndArray(const Float_t f)
{

  FillArray(f);
  return FillHist(f);
}

// --------------------------------------------------------------------------
//
// Fills fHGausHist with f
// Returns kFALSE, if overflow or underflow occurred, else kTRUE
//
Bool_t MHGausEvents::FillHist(const Float_t f)
{
  return fHGausHist.Fill(f) > -1;
}

// --------------------------------------------------------------------------
//
// Fill fEvents with f
// If size of fEvents is 0, initializes it to 512
// If size of fEvents is smaller than fCurrentSize, double the size
// Increase fCurrentSize by 1
//
void MHGausEvents::FillArray(const Float_t f)
{
  if (fEvents.GetSize() == 0)
    fEvents.Set(512);

  if (fCurrentSize >= fEvents.GetSize())
    fEvents.Set(fEvents.GetSize()*2);
  
  fEvents.AddAt(f,fCurrentSize++);
}


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

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


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


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


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


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


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


const Bool_t MHGausEvents::IsEmpty() const
{
    return !(fHGausHist.GetEntries());
}


const Bool_t MHGausEvents::IsFourierSpectrumOK() const 
{
  return TESTBIT(fFlags,kFourierSpectrumOK);
}


const Bool_t MHGausEvents::IsGausFitOK() const 
{
  return TESTBIT(fFlags,kGausFitOK);
}


const Bool_t MHGausEvents::IsExpFitOK() const 
{
  return TESTBIT(fFlags,kExpFitOK);
}


// -------------------------------------------------------------------
//
// The flag setters are to be used ONLY for Monte-Carlo!!
//
void  MHGausEvents::SetGausFitOK(const Bool_t b)
{
  b ? SETBIT(fFlags,kGausFitOK) : CLRBIT(fFlags,kGausFitOK);

}
// -------------------------------------------------------------------
//
// The flag setters are to be used ONLY for Monte-Carlo!!
//
void  MHGausEvents::SetExpFitOK(const Bool_t b)
{
  
  b ? SETBIT(fFlags,kExpFitOK) : CLRBIT(fFlags,kExpFitOK);  
}

// -------------------------------------------------------------------
//
// The flag setters are to be used ONLY for Monte-Carlo!!
//
void  MHGausEvents::SetFourierSpectrumOK(const Bool_t b)
{

  b ? SETBIT(fFlags,kFourierSpectrumOK) : CLRBIT(fFlags,kFourierSpectrumOK);    
}


// -------------------------------------------------------------------
//
// Create the fourier spectrum using the class MFFT.
// The result is projected into a histogram and fitted by an exponential
// 
void MHGausEvents::CreateFourierSpectrum()
{

  if (fFExpFit)
    return;

  //
  // 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");
  fHPowerProbability->SetXTitle("P(f)");
  fHPowerProbability->SetDirectory(NULL);
  //
  // 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)
    SetExpFitOK(kTRUE);
  
  //
  // For the moment, this is the only check, later we can add more...
  //
  SetFourierSpectrumOK(IsExpFitOK());

  return;
}

// -------------------------------------------------------------------
//
// Fit fGausHist with a Gaussian after stripping zeros from both ends 
// and rebinned to the number of bins specified in fBinsAfterStripping
//
// The fit results are retrieved and stored in class-own variables.  
//
// A flag IsGausFitOK() is set according to whether the fit probability 
// is smaller or bigger than fProbLimit, whether the NDF is bigger than 
// fNDFLimit and whether results are NaNs.
//
Bool_t MHGausEvents::FitGaus(Option_t *option, const Double_t xmin, const Double_t xmax)
{

  if (IsGausFitOK())
    return kTRUE;

  //
  // First, strip the zeros from the edges which contain only zeros and rebin 
  // to about fBinsAfterStripping bins. 
  //
  // (ATTENTION: The Chisquare method is more sensitive, 
  // the _less_ bins, you have!)
  //
  StripZeros(&fHGausHist,fBinsAfterStripping);
  
  //
  // Get the fitting ranges
  //
  Axis_t rmin = (xmin==0.) && (xmax==0.) ? fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()) : xmin;
  Axis_t rmax = (xmin==0.) && (xmax==0.) ? fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast())  : xmax;

  //
  // 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 = fHGausHist.GetRMS();
  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.,area_guess*1.5);
  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) 
      || fFGausFit->GetNDF() < fNDFLimit 
      || fProb < fProbLimit )
    return kFALSE;
  
  SetGausFitOK(kTRUE);
  return kTRUE;
}

// -----------------------------------------------------------------------------
// 
// Bypasses the Gauss fit by taking mean and RMS from the histogram
//
// Errors are determined in the following way:
// MeanErr  = RMS / Sqrt(entries)
// SigmaErr = RMS / (2.*Sqrt(entries) )
//
void MHGausEvents::BypassFit()
{

  const Stat_t entries = fHGausHist.GetEntries();
  
  if (entries <= 0.)
    {
      *fLog << warn << GetDescriptor() << ": Cannot bypass fit. Number of entries smaller or equal 0" << endl;
      return;
    }
  
  SetMean     ( fHGausHist.GetMean() );
  SetMeanErr  ( fHGausHist.GetRMS() / TMath::Sqrt(entries) );
  SetSigma    ( fHGausHist.GetRMS() );
  SetSigmaErr ( fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2. );
}

// -----------------------------------------------------------------------------------
// 
// A default print
//
void MHGausEvents::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;
  
}

// ----------------------------------------------------------------------------------
//
// Create a graph to display the array fEvents
// If the variable fEventFrequency is set, the x-axis is transformed into real time.
//
void MHGausEvents::CreateGraphEvents()
{

  MArray::StripZeros(fEvents);

  const Int_t n = fEvents.GetSize();

  fGraphEvents = new TGraph(n,CreateEventXaxis(n),fEvents.GetArray());  
  fGraphEvents->SetTitle("Evolution of Events with time");
  fGraphEvents->GetXaxis()->SetTitle((fEventFrequency) ? "Time [s]" : "Event Nr.");
}

// ----------------------------------------------------------------------------------
//
// Create a graph to display the array fPowerSpectrum
// If the variable fEventFrequency is set, the x-axis is transformed into real frequency.
//
void MHGausEvents::CreateGraphPowerSpectrum()
{

  MArray::StripZeros(*fPowerSpectrum);

  const Int_t n = fPowerSpectrum->GetSize();

  fGraphPowerSpectrum = new TGraph(n,CreatePSDXaxis(n),fPowerSpectrum->GetArray());  
  fGraphPowerSpectrum->SetTitle("Power Spectrum Density");
  fGraphPowerSpectrum->GetXaxis()->SetTitle((fEventFrequency) ? "Frequency [Hz]" : "Frequency");
  fGraphPowerSpectrum->GetYaxis()->SetTitle("P(f)");
}

// -----------------------------------------------------------------------------
// 
// Create the x-axis for the event graph
//
Float_t *MHGausEvents::CreateEventXaxis(Int_t n)
{

  Float_t *xaxis = new Float_t[n];

  if (fEventFrequency)
    for (Int_t i=0;i<n;i++)
      xaxis[i] = (Float_t)i/fEventFrequency;
  else
    for (Int_t i=0;i<n;i++)
      xaxis[i] = (Float_t)i;

  return xaxis;
                 
}

// -----------------------------------------------------------------------------
// 
// Create the x-axis for the event graph
//
Float_t *MHGausEvents::CreatePSDXaxis(Int_t n)
{

  Float_t *xaxis = new Float_t[n];

  if (fEventFrequency)
    for (Int_t i=0;i<n;i++)
      xaxis[i] = 0.5*(Float_t)i*fEventFrequency/n;
  else
    for (Int_t i=0;i<n;i++)
      xaxis[i] = 0.5*(Float_t)i/n;

  return xaxis;
                 
}
  
// -----------------------------------------------------------------------------
// 
// Default draw:
//
// The following options can be chosen:
//
// "EVENTS": displays a TGraph of the array fEvents
// "FOURIER": display a TGraph of the fourier transform of fEvents 
//            displays the projection of the fourier transform with the fit
//
void MHGausEvents::Draw(const Option_t *opt)
{

  TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this,600, 900);

  TString option(opt);
  option.ToLower();
  
  Int_t win = 1;

  if (option.Contains("events"))
    {
      option.ReplaceAll("events","");
      win += 1;
    }
  if (option.Contains("fourier"))
    {
      option.ReplaceAll("fourier","");
      win += 2;
    }
  
  pad->SetTicks();
  pad->SetBorderMode(0);
  pad->Divide(1,win);
  pad->cd(1);

  if (!IsEmpty())
    gPad->SetLogy();

  fHGausHist.Draw(opt);

  if (fFGausFit)
    {
      fFGausFit->SetLineColor(IsGausFitOK() ? kGreen : kRed);
      fFGausFit->Draw("same");
    }
  switch (win)
    {
    case 2:
      pad->cd(2);
      DrawEvents();
      break;
    case 3:
      pad->cd(2);
      DrawPowerSpectrum(*pad,3);
      break;
    case 4:
      pad->cd(2);
      DrawEvents();
      pad->cd(3);
      DrawPowerSpectrum(*pad,4);
      break;
    }
}

void MHGausEvents::DrawEvents()
{
  
  if (!fGraphEvents)
    CreateGraphEvents();

  fGraphEvents->SetBit(kCanDelete);
  fGraphEvents->SetTitle("Events with time");
  fGraphEvents->Draw("AL");
  
}


void MHGausEvents::DrawPowerSpectrum(TVirtualPad &pad, Int_t i)
{
  
  if (fPowerSpectrum)
    {
      if (!fGraphPowerSpectrum)
        CreateGraphPowerSpectrum();
      
      fGraphPowerSpectrum->Draw("AL");          
      fGraphPowerSpectrum->SetBit(kCanDelete);
    }
  
  pad.cd(i);

  if (fHPowerProbability && fHPowerProbability->GetEntries() > 0)
    {
      gPad->SetLogy();
      fHPowerProbability->Draw();
      if (fFExpFit)
        {
          fFExpFit->SetLineColor(IsExpFitOK() ? kGreen : kRed);
          fFExpFit->Draw("same");
        }
    }
}