source: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.cc@ 3622

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */
//////////////////////////////////////////////////////////////////////////////
//
//  MHCalibrationChargePix
//
//  Histogram class for charge calibration analysis. Holds the histogrammed
//  summed FADC slices and some rough absolute timing. Calculates the mean 
//  sum of FADC slices and its sigma and perform a Fourier analysis.
//
//////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationChargePix.h"

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

#include <TVirtualPad.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TGraph.h>

#include "MH.h"

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

ClassImp(MHCalibrationChargePix);

using namespace std;

const Int_t   MHCalibrationChargePix::fgChargeNbins     = 2000;
const Axis_t  MHCalibrationChargePix::fgChargeFirst     = -0.5;
const Axis_t  MHCalibrationChargePix::fgChargeLast      = 1999.5;
const Int_t   MHCalibrationChargePix::fgAbsTimeNbins    = 15;
const Axis_t  MHCalibrationChargePix::fgAbsTimeFirst    = -0.5;
const Axis_t  MHCalibrationChargePix::fgAbsTimeLast     = 14.5;
const Float_t MHCalibrationChargePix::fgPickupLimit     = 5.;
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
// Sets: 
// - the default number for fChargeNbins  (fgChargeNbins)
// - the default number for fChargeFirst  (fgChargeFirst)
// - the default number for fChargeLast   (fgChargeLast)
// - the default number for fAbsTimeNbins (fgAbstTimeNbins)
// - the default number for fAbsTimeFirst (fgAbsTimeFirst)
// - the default number for fAbsTimeLast  (fgAbsTimeLast)
// - the default number for fPickupLimit  (fgPickupLimit)
//
// - the default name of the  fHGausHist ("HCalibrationCharge")
// - the default title of the fHGausHist ("Distribution of Summed FADC slices Pixel ")
// - the default x-axis title for fHGausHist ("Sum FADC Slices")
// - the default y-axis title for fHGausHist ("Nr. of events")
//
// - the default name of the  fHAbsTime ("HAbsTimePixel")
// - the default title of the fHAbsTime ("Distribution of Absolute Arrival Times Pixel ")
// - the default x-axis title for fHAbsTime ("Absolute Arrival Time [FADC slice nr]")
// - the default y-axis title for fHAbsTime ("Nr. of events"); 
// - the default directory of the fHAbsTime (NULL)
//
// Initializes:
// - fHAbsTime()
// - fPixId to -1
// - all variables to 0.
// - all flags to kFALSE
//
MHCalibrationChargePix::MHCalibrationChargePix(const char *name, const char *title)
    : fHAbsTime(), fFlags(0)
{ 
  
  fName  = name  ? name  : "MHCalibrationChargePix";
  fTitle = title ? title : "Statistics of the FADC sums of calibration events";

  SetChargeNbins();
  SetChargeFirst();
  SetChargeLast();

  SetAbsTimeNbins();
  SetAbsTimeFirst();
  SetAbsTimeLast();

  SetPickupLimit();

  fHGausHist.SetName("HCalibrationCharge");
  fHGausHist.SetTitle("Distribution of Summed FADC slices Pixel");  
  fHGausHist.SetXTitle("Sum FADC Slices");
  fHGausHist.SetYTitle("Nr. of events");

  fHAbsTime.SetName("HAbsTimePixel");
  fHAbsTime.SetTitle("Distribution of Absolute Arrival Times Pixel ");  
  fHAbsTime.SetXTitle("Absolute Arrival Time [FADC slice nr]");
  fHAbsTime.SetYTitle("Nr. of events");

  fHAbsTime.UseCurrentStyle();
  fHAbsTime.SetDirectory(NULL); 

  Clear();
  
}

// --------------------------------------------------------------------------
//
// Sets:
// - fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast);
// - fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast);
//
void MHCalibrationChargePix::Init()
{

  fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast);
  fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast);
}

// --------------------------------------------------------------------------
//
// Sets:
// - all variables to 0., except fPixId to -1
// - all flags to kFALSE
// 
// - executes MHGausEvents::Clear()
//
void MHCalibrationChargePix::Clear(Option_t *o)
{

  fPixId     = -1;
  fSaturated = 0.;
  fPickup    = 0.;

  SetExcluded ( kFALSE );
  
  MHGausEvents::Clear();
  return;
}

// --------------------------------------------------------------------------
//
// Empty function to overload MHGausEvents::Reset()
//
void MHCalibrationChargePix::Reset()
{
}

// --------------------------------------------------------------------------
//
// - Set fPixId to id
// - Add id to name and title of this
// - Add id to names of fHGausHist and fHAbsTime
// - Add id to the titles of fHGausHist and fHAbsTime
//
void MHCalibrationChargePix::ChangeHistId(Int_t id)
{

  fPixId = id;

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

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

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


void MHCalibrationChargePix::SetExcluded(const Bool_t b)
{
    b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded);
}

// --------------------------------------------------------------------------
//
// returns fHGausHist.Integral("width")
//
const Float_t MHCalibrationChargePix::GetIntegral() const 
{ 
   return fHGausHist.Integral("width");  
}

const Float_t MHCalibrationChargePix::GetAbsTimeMean() const 
{
  return fHAbsTime.GetMean();
}

const Float_t MHCalibrationChargePix::GetAbsTimeRms()  const 
{
  return fHAbsTime.GetRMS();
}

const Bool_t MHCalibrationChargePix::IsExcluded()      const
{
  return TESTBIT(fFlags,kExcluded);
}

// --------------------------------------------------------------------------
//
// Fills fHAbsTime with t
// Returns kFALSE, if overflow or underflow occurred, else kTRUE
//
Bool_t MHCalibrationChargePix::FillAbsTime(Float_t t)
{
  return fHAbsTime.Fill(t) > -1;
}

// -----------------------------------------------------------------------------
// 
// Default draw:
//
// The following options can be chosen:
//
// "": displays the fHGausHist and the fHAbsTime
// "all": executes additionally MHGausEvents::Draw(), with options
//
void MHCalibrationChargePix::Draw(const Option_t *opt)
{

  TString option(opt);
  option.ToLower();
  
  Int_t win = 1;
  
  TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,600, 600);
  TVirtualPad *pad    = NULL;

  if (option.Contains("all"))
  {
      option.ReplaceAll("all","");
      oldpad->Divide(2,1);
      win = 2;
      oldpad->cd(1);
      TVirtualPad *newpad = gPad;
      pad = newpad;
      pad->Divide(1,2);
      pad->cd(1);
  }
  else
  {
      pad = oldpad;
      pad->Divide(1,2);
      pad->cd(1);
  }

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

  gPad->SetTicks();

  fHGausHist.GetXaxis()->SetLabelSize(0.06);
  fHGausHist.GetYaxis()->SetLabelSize(0.07);
  fHGausHist.GetXaxis()->SetLabelOffset(0.01);
  fHGausHist.GetYaxis()->SetLabelOffset(0.01);
  fHGausHist.GetXaxis()->SetTitleSize(0.065);
  fHGausHist.GetYaxis()->SetTitleSize(0.07);
  fHGausHist.GetXaxis()->SetTitleOffset(0.6);
  fHGausHist.GetYaxis()->SetTitleOffset(0.6);
  fHGausHist.Draw(opt); 
  if (fFGausFit)
  {
      fFGausFit->SetLineColor(IsGausFitOK() ? kGreen : kRed);
      fFGausFit->Draw("same");
  }

  pad->cd(2);
  gPad->SetTicks();

  fHAbsTime.GetXaxis()->SetLabelSize(0.06);
  fHAbsTime.GetYaxis()->SetLabelSize(0.07);
  fHAbsTime.GetXaxis()->SetLabelOffset(0.01);
  fHAbsTime.GetYaxis()->SetLabelOffset(0.01);
  fHAbsTime.GetXaxis()->SetTitleSize(0.065);
  fHAbsTime.GetYaxis()->SetTitleSize(0.07);
  fHAbsTime.GetXaxis()->SetTitleOffset(0.6);
  fHAbsTime.GetYaxis()->SetTitleOffset(0.6);
  fHAbsTime.Draw(opt);

  if (win < 2)
      return;

  oldpad->cd(2);
  MHGausEvents::Draw("fourierevents");

}

// -----------------------------------------------------------------------------
// 
// Repeats the Gauss fit in a smaller range, defined by: 
// 
// min = GetMean() - fPickupLimit * GetSigma();
// max = GetMean() + fPickupLimit * GetSigma();
//
Bool_t MHCalibrationChargePix::RepeatFit(const Option_t *option)
{

  //
  // Get new fitting ranges
  //
  Axis_t rmin = GetMean() - fPickupLimit * GetSigma();
  Axis_t rmax = GetMean() + fPickupLimit * GetSigma();

  GetFGausFit()->SetRange(rmin,rmax);

  GetHGausHist()->Fit(GetFGausFit(),option);

  SetMean     ( GetFGausFit()->GetParameter(1) );
  SetMeanErr  ( GetFGausFit()->GetParameter(2) );
  SetSigma    ( GetFGausFit()->GetParError(1)  ); 
  SetSigmaErr ( GetFGausFit()->GetParError(2)  ); 
  SetProb     ( GetFGausFit()->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 ( GetMean()     ) 
      || TMath::IsNaN ( GetMeanErr()  )
      || TMath::IsNaN ( GetProb()     )    
      || TMath::IsNaN ( GetSigma()    )
      || TMath::IsNaN ( GetSigmaErr() ) 
      || GetFGausFit()->GetNDF() < fNDFLimit 
      || GetProb() < 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 MHCalibrationChargePix::BypassFit()
{

  //
  // In case, we do not obtain reasonable values 
  // with the fit, we take the histogram values
  //
  SetMean     ( fHGausHist.GetMean() );
  SetMeanErr  ( fHGausHist.GetRMS() / TMath::Sqrt(fHGausHist.GetEntries()) );
  SetSigma    ( fHGausHist.GetRMS() );
  SetSigmaErr ( fHGausHist.GetRMS() / TMath::Sqrt(fHGausHist.GetEntries()) / 2. );
}

void MHCalibrationChargePix::CountPickup()
{
   fPickup  = (Int_t)GetHGausHist()->Integral(GetHGausHist()->GetXaxis()->FindBin(GetMean()+fPickupLimit*GetSigma()),
                                               GetHGausHist()->GetXaxis()->GetLast(), 
                                               "width");
}