source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationPix.cc@ 4985

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

//////////////////////////////////////////////////////////////////////////////
//
//  MHCalibrationPix
//
//  A base class for events which are believed to follow a Gaussian distribution 
//  with time, e.g. calibration events, observables containing white noise, ...
//
//  MHCalibrationPix derives from MHGausEvents, thus all features of 
//  MHGausEvents can be used by a class deriving from MHCalibrationPix
//
//  See also: MHGausEvents
//
//////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationPix.h"

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

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

ClassImp(MHCalibrationPix);

using namespace std;

const Float_t  MHCalibrationPix::fgBlackoutLimit        = 5.;
const Float_t  MHCalibrationPix::fgPickupLimit          = 5.;
// --------------------------------------------------------------------------
//
// Default Constructor. 
// Sets: 
// - the default number for fPickupLimit           (fgPickupLimit)
// - the default number for fBlackoutLimit         (fgBlackoutLimit)
//
// Initializes:
// - all variables to 0.
//
MHCalibrationPix::MHCalibrationPix(const char *name, const char *title)
    : fEventFrequency(0), fPixId(-1)
{ 

  fName  = name  ? name  : "MHCalibrationPix";
  fTitle = title ? title : "Calibration histogram events";

  Clear();
  
  SetBlackoutLimit();
  SetPickupLimit();
}


      
// --------------------------------------------------------------------------
//
// Default Clear(), can be overloaded.
//
// Sets:
// - all other pointers to NULL
// - all variables to 0., except fPixId to -1 and keep fEventFrequency
// - all flags to kFALSE
// 
// Deletes (if not NULL):
// - all pointers
//
void MHCalibrationPix::Clear(Option_t *o)
{

  MHGausEvents::Clear();
  fSaturated         = 0;
}


#if 0
// --------------------------------------------------------------------------
//
// ATTENTION: This nasty Clone function is necessary since the ROOT cloning 
//            lead to crashes on SOME machines (unfortunately not mine...).
//            This function is a workaround in order to achieve the correct 
//            DrawClone() behaviour.
// 
TObject *MHCalibrationPix::Clone(const char *name) const
{

  MHCalibrationPix &pix = *new MHCalibrationPix(name ? name : fName.Data(),fTitle.Data());
  
  // 
  // Copy MHGausEvents data members
  //
  pix.fBinsAfterStripping   = fBinsAfterStripping;
  pix.fCurrentSize          = fCurrentSize;
  pix.fFlags                = fFlags;
  pix.fPowerProbabilityBins = fPowerProbabilityBins;

  if (fHPowerProbability)
    pix.fHPowerProbability=(TH1I*)fHPowerProbability->Clone();
  
  if (fPowerSpectrum)
    pix.fPowerSpectrum = new TArrayF(*fPowerSpectrum);

  pix.fEvents = fEvents;

  if (fFGausFit)
    pix.fFGausFit=(TF1*)fFGausFit->Clone();
  if (fFExpFit)
    pix.fFExpFit=(TF1*)fFExpFit->Clone();

  pix.fFirst = fFirst;

  if (fGraphEvents)
    pix.fGraphEvents=(TGraph*)fGraphEvents->Clone();
  if (fGraphPowerSpectrum)
    pix.fGraphPowerSpectrum=(TGraph*)fGraphPowerSpectrum->Clone();

  pix.fHGausHist = fHGausHist;

  pix.fLast      = fLast;
  pix.fMean      = fMean;
  pix.fMeanErr   = fMeanErr;
  pix.fNbins     = fNbins;
  pix.fNDFLimit  = fNDFLimit;
  pix.fSigma     = fSigma;
  pix.fSigmaErr  = fSigmaErr;
  pix.fProb      = fProb;
  pix.fProbLimit = fProbLimit;

  // 
  // Copy data members
  //
  pix.fEventFrequency       = fEventFrequency;
  pix.fBlackoutLimit        = fBlackoutLimit; 
  pix.fSaturated            = fSaturated;     
  pix.fPickupLimit          = fPickupLimit;   
  pix.fPixId                = fPixId;        

  return &pix;
}
#endif

// -----------------------------------------------------------------------------
// 
// 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 MHCalibrationPix::BypassFit()
{

  const Stat_t entries = fHGausHist.GetEntries();
  
  if (entries <= 0.)
    {
      *fLog << warn << GetDescriptor() 
            << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl;
      return;
    }
  
  fMean     = fHGausHist.GetMean();
  fMeanErr  = fHGausHist.GetRMS() / TMath::Sqrt(entries);
  fSigma    = fHGausHist.GetRMS() ;
  fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.;
}

// --------------------------------------------------------------------------
//
// - Set fPixId to id
//
// Add id to names and titles of:
// - fHGausHist
//
void MHCalibrationPix::ChangeHistId(const Int_t id)
{

  fPixId = id;

  fHGausHist.SetName(  Form("%s%d", fHGausHist.GetName(),  id));
  fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id));

  fName  = Form("%s%d", fName.Data(),  id);
  fTitle = Form("%s%d", fTitle.Data(), id);

}


// -----------------------------------------------------------------------------
// 
// Create the x-axis for the event graph
//
Float_t *MHCalibrationPix::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 *MHCalibrationPix::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;
                 
}

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

  MHGausEvents::CreateGraphEvents();
  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 MHCalibrationPix::CreateGraphPowerSpectrum()
{

  MHGausEvents::CreateGraphPowerSpectrum();
  
  fGraphPowerSpectrum->GetXaxis()->SetTitle((fEventFrequency) ? "Frequency [Hz]" : "Frequency");
}

// -------------------------------------------------------------------------------
//
// Return the number of "blackout" events, which are events with values higher 
// than fBlackoutLimit sigmas from the mean
//
//
const Double_t MHCalibrationPix::GetBlackout() const 
{
  
  if ((fMean == 0.) && (fSigma == 0.))
    return -1.;

  const Int_t first = fHGausHist.GetXaxis()->GetFirst();
  const Int_t last  = fHGausHist.GetXaxis()->FindBin(fMean-fBlackoutLimit*fSigma);

  if (first >= last)
    return 0.;
  
  return fHGausHist.Integral(first, last, "width");
}


// -------------------------------------------------------------------------------
//
// Return the number of "pickup" events, which are events with values higher 
// than fPickupLimit sigmas from the mean
//
//
const Double_t MHCalibrationPix::GetPickup() const 
{
  
  if ((fMean == 0.) && (fSigma == 0.))
    return -1.;

  const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma);
  const Int_t last  = fHGausHist.GetXaxis()->GetLast();

  if (first >= last)
    return 0.;
  
  return fHGausHist.Integral(first, last, "width");
}


// --------------------------------------------------------------------------
//
// Re-normalize the results, has to be overloaded
//
void  MHCalibrationPix::Renorm()
{
}

// -----------------------------------------------------------------------------
//
// If flag IsGausFitOK() is set (histogram already successfully fitted), 
// returns kTRUE
// 
// If both fMean and fSigma are still zero, call FitGaus() 
// 
// Repeats the Gauss fit in a smaller range, defined by: 
// 
// min = GetMean() - fBlackoutLimit * GetSigma();
// max = GetMean() + fPickupLimit   * GetSigma();
//
// 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 MHCalibrationPix::RepeatFit(const Option_t *option)
{

  if (IsGausFitOK())
    return kTRUE;

  if ((fMean == 0.) && (fSigma == 0.))
    return FitGaus();

  //
  // Get new fitting ranges
  //
  Axis_t rmin = fMean - fBlackoutLimit * fSigma;
  Axis_t rmax = fMean + fPickupLimit   * fSigma;

  Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
  Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ;

  fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax);

  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 (default: fgNDFLimit)
  // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
  //
  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;

}