source: trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.cc@ 3033

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

//////////////////////////////////////////////////////////////////////////////
//
//  MHPedestalPixel
//
//  Performs all the necessary fits to extract the mean number of photons
//              out of the derived light flux
//
//////////////////////////////////////////////////////////////////////////////
#include "MHPedestalPixel.h"

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

#include <TStyle.h>
#include <TCanvas.h>

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

ClassImp(MHPedestalPixel);

using namespace std;
// Square root of 2pi:
const Float_t gkSq2Pi = 2.506628274631;
const Float_t gkProbLimit = 0.01;
const Int_t   MHPedestalPixel::gkChargeNbins    = 500 ;
const Int_t   MHPedestalPixel::gkChargevsNbins  = 1000;
const Axis_t  MHPedestalPixel::gkChargevsNFirst = -0.5;
const Axis_t  MHPedestalPixel::gkChargevsNLast  = 999.5;
const Axis_t  MHPedestalPixel::gkChargeFirst    = -0.5;
const Axis_t  MHPedestalPixel::gkChargeLast     = 499.5;

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHPedestalPixel::MHPedestalPixel()
      : fPixId(-1),
        fGausFit(NULL)
{ 
  
  // Create a large number of bins, later we will rebin
  fHPedestalCharge = new TH1F("HPedestalCharge","Distribution of Summed FADC Pedestal Slices Pixel ",
                              gkChargeNbins,gkChargeFirst,gkChargeLast);
  fHPedestalCharge->SetXTitle("Sum FADC Slices");
  fHPedestalCharge->SetYTitle("Nr. of events");
  fHPedestalCharge->Sumw2();
  
  // We define a reasonable number and later enlarge it if necessary
  fHPedestalChargevsN = new TH1I("HChargevsN","Sum of Charges vs. Event Number Pixel ",
                                 gkChargevsNbins,gkChargevsNFirst,gkChargevsNLast);
  fHPedestalChargevsN->SetXTitle("Event Nr.");
  fHPedestalChargevsN->SetYTitle("Sum of FADC slices");
  
  fHPedestalCharge->SetDirectory(NULL);
  fHPedestalChargevsN->SetDirectory(NULL);
  
  Clear();
}

MHPedestalPixel::~MHPedestalPixel()
{

  delete fHPedestalCharge;
  delete fHPedestalChargevsN;
  
  if (fGausFit)
    delete fGausFit;
}


void MHPedestalPixel::Clear(Option_t *o)
{
  
  fTotalEntries    = 0;

  fChargeMean      = -1.;
  fChargeMeanErr   = -1.;
  fChargeSigma     = -1;
  fChargeSigmaErr  = -1;

  fChargeChisquare = -1.;
  fChargeProb      = -1.;
  fChargeNdf       = -1;

  fChargeFirst     = gkChargeFirst;
  fChargeLast      = gkChargeLast;

  if (fGausFit)
    delete fGausFit;

  CLRBIT(fFlags,kFitted);
  CLRBIT(fFlags,kFitOK);
  CLRBIT(fFlags,kOscillating);

  return;
}


void MHPedestalPixel::Reset()
{
  
  Clear();
  
  fHPedestalCharge->Reset();
  fHPedestalChargevsN->Reset();

}



Bool_t MHPedestalPixel::IsEmpty() const
{
    return !fHPedestalCharge->GetEntries();
}

Bool_t MHPedestalPixel::IsFitOK() const 
{
    return TESTBIT(fFlags,kFitOK);
}

Bool_t MHPedestalPixel::FillCharge(Float_t q)
{
    return (fHPedestalCharge->Fill(q) > -1);
}

Bool_t MHPedestalPixel::FillChargevsN(Float_t q)
{

  fTotalEntries++;
  return (fHPedestalChargevsN->Fill(fTotalEntries,q) > -1);
}

void MHPedestalPixel::ChangeHistId(Int_t id)
{

  fPixId = id;

  fHPedestalCharge->SetName(Form("%s%d", fHPedestalCharge->GetName(), id));
  fHPedestalChargevsN->SetName(Form("%s%d", fHPedestalChargevsN->GetName(), id));
  fHPedestalCharge->SetTitle(Form("%s%d", fHPedestalCharge->GetTitle(), id));
  fHPedestalChargevsN->SetTitle(Form("%s%d", fHPedestalChargevsN->GetTitle(), id));

}

TObject *MHPedestalPixel::DrawClone(Option_t *option) const
{

  gROOT->SetSelectedPad(NULL);
  
  MHPedestalPixel *newobj = (MHPedestalPixel*)Clone();

  if (!newobj) 
    return 0;
  newobj->SetBit(kCanDelete);


  if (strlen(option)) 
    newobj->Draw(option);
  else    
    newobj->Draw(GetDrawOption());
  
  return newobj;
}
  

// -------------------------------------------------------------------------
//
// Draw the histogram
//
void MHPedestalPixel::Draw(Option_t *opt) 
{

  gStyle->SetOptFit(1);
  gStyle->SetOptStat(111111);
  
  gROOT->SetSelectedPad(NULL);
  
  TCanvas *c = MH::MakeDefCanvas(this,600,900); 
  
  c->Divide(1,2);
  
  c->cd(1);
  gPad->SetBorderMode(0);
  gPad->SetTicks();

  if (fHPedestalCharge->Integral() > 0)
    gPad->SetLogy(1);

  fHPedestalCharge->Draw(opt);
  
  c->Modified();
  c->Update();

  if (fGausFit)
    {
      fGausFit->SetLineColor(IsFitOK() ? kGreen : kRed);
      fGausFit->Draw("same");
    }

  c->Modified();
  c->Update();

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

  fHPedestalChargevsN->Draw(opt);
  c->Modified();
  c->Update();

  return;
}


// --------------------------------------------------------------------------
//
// 1) Return if the charge distribution is already succesfully fitted  
//    or if the histogram is empty
// 2) Cut the histograms empty edges
// 3) Fit the histograms with a Gaussian
// 4) In case of failure print out the fit results
// 5) Retrieve the results and store them in this class
//
Bool_t MHPedestalPixel::FitCharge(Option_t *option)
{

  if (fGausFit)
    return kFALSE;

  //
  // Get the fitting ranges
  //
  Axis_t rmin = fChargeFirst;
  Axis_t rmax = fChargeLast;

  //
  // First guesses for the fit (should be as close to reality as possible, 
  // otherwise the fit goes gaga because of high number of dimensions ...
  //
  const Stat_t   entries     = fHPedestalCharge->Integral();
  const Double_t mu_guess    = fHPedestalCharge->GetBinCenter(fHPedestalCharge->GetMaximumBin());
  const Double_t sigma_guess = mu_guess/15.;
  const Double_t area_guess  = entries/gkSq2Pi/sigma_guess;

  fGausFit = new TF1(Form("%s%d","GausFit ",fPixId),"gaus",rmin,rmax);

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

  fHPedestalCharge->Fit(fGausFit,option);
  
  // 
  // If we are not able to fit, try once again
  //
  if (fGausFit->GetProb() < gkProbLimit)
    {

      Axis_t rtry = fGausFit->GetParameter(1) - 2.0*fGausFit->GetParameter(2);
      rmin        = (rtry < rmin ? rmin : rtry);
      rmax        = fGausFit->GetParameter(1) + 2.0*fGausFit->GetParameter(2);
      fGausFit->SetRange(rmin,rmax);  

      fHPedestalCharge->Fit(fGausFit,option);
    }
  
  fChargeChisquare = fGausFit->GetChisquare();
  fChargeNdf       = fGausFit->GetNDF();
  fChargeProb      = fGausFit->GetProb();
  fChargeMean      = fGausFit->GetParameter(1);
  fChargeMeanErr   = fGausFit->GetParError(1);
  fChargeSigma     = fGausFit->GetParameter(2);
  fChargeSigmaErr  = fGausFit->GetParError(2);

  SETBIT(fFlags,kFitted);
  //
  // The fit result is accepted under condition:
  // The Results are not nan's
  // The Probability is greater than gkProbLimit (default 0.001 == 99.9%)
  //
  if (TMath::IsNaN(fChargeMean) || TMath::IsNaN(fChargeMeanErr))
    {
      CLRBIT(fFlags,kFitOK);
      return kFALSE;
    }
  
  if ((fChargeProb < gkProbLimit) || (TMath::IsNaN(fChargeProb)))
    {
      CLRBIT(fFlags,kFitOK);
      return kFALSE;
    }
  
  SETBIT(fFlags,kFitOK);
  return kTRUE;
}

void MHPedestalPixel::CutAllEdges()
{

  Int_t nbins = 30;

  MH::CutEdges(fHPedestalCharge,nbins);

  fChargeFirst = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetFirst());
  fChargeLast  = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetLast())
                +fHPedestalCharge->GetBinWidth(0);

  MH::CutEdges(fHPedestalChargevsN,0);

}

void MHPedestalPixel::Print(const Option_t *o) const
{
  
  *fLog << all << "Pedestal Fits Pixel: "  << fPixId  << endl;

  if (!TESTBIT(fFlags,kFitted))
    {
      gLog << "Pedestal Histogram has not yet been fitted" << endl;
      return;
    }
  
  gLog << "Results of the Summed Charges Fit: "   << endl;
  gLog << "Chisquare:      " << fChargeChisquare  << endl;
  gLog << "DoF:            " << fChargeNdf        << endl;
  gLog << "Probability:    " << fChargeProb       << endl;
  gLog << "Results of fit: " << (TESTBIT(fFlags,kFitOK) ? "Ok.": "Not OK!") << endl;
    
}