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

/* ======================================================================== *\
!
! *
! * 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
//
//////////////////////////////////////////////////////////////////////////////
#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.;
const Int_t   MHCalibrationChargePix::fgPulserFrequency = 500;
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationChargePix::MHCalibrationChargePix(const char *name, const char *title)
    : fHAbsTime(), fFlags(0)
{ 
  
  fName  = name  ? name  : "MHCalibrationChargePix";
  fTitle = title ? title : "Fill the FADC sums of calibration events and perform the fits";

  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();
  
  SetPulserFrequency();
}

void MHCalibrationChargePix::Init()
{

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


void MHCalibrationChargePix::Clear(Option_t *o)
{

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

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


void MHCalibrationChargePix::Reset()
{
}


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::SetPulserFrequency(Float_t f)
{
  SetEventFrequency(f);
}


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

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();
}

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

Bool_t MHCalibrationChargePix::FillAbsTime(Float_t t)
{
  return fHAbsTime.Fill(t) > -1;
}


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");

}

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;

}

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() / fHGausHist.GetEntries() );
  SetSigma    ( fHGausHist.GetRMS() );
  SetSigmaErr ( fHGausHist.GetRMS() / fHGausHist.GetEntries() / 2. );
}

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