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

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

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHPedestalPixel::MHPedestalPixel(const char *name, const char *title)
      : fPixId(-1),
        fChargeNbins(500),
        fChargevsNbins(1000),
        fChargevsNFirst(-0.5),
        fChargevsNLast(999.5),        
        fChargeFirst(-0.5),
        fChargeLast(99.5),
        fGausFit(NULL)
{ 

    fName  = name  ? name  : "MHPedestalPixel";
    fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits";

    // Create a large number of bins, later we will rebin
    fHPedestalCharge = new TH1F("HPedestalCharge","Distribution of Summed FADC Pedestal Slices Pixel ",
                        fChargeNbins,fChargeFirst,fChargeLast);
    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 ",
                           fChargevsNbins,fChargevsNFirst,fChargevsNLast);
    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     = -0.5;
  fChargeLast      = 99.5;

  if (fGausFit)
    delete fGausFit;

  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;

  TString nameQ = TString(fHPedestalCharge->GetName());
  nameQ += id;
  fHPedestalCharge->SetName(nameQ.Data());

  TString nameQvsN  = TString(fHPedestalChargevsN->GetName());
  nameQvsN += id;
  fHPedestalChargevsN->SetName(nameQvsN.Data());

  TString titleQ = TString(fHPedestalCharge->GetTitle());
  titleQ += id;
  fHPedestalCharge->SetTitle(titleQ.Data());

  TString titleQvsN  = TString(fHPedestalChargevsN->GetTitle());
  titleQvsN += id;
  fHPedestalChargevsN->SetTitle(titleQvsN.Data());

}


Bool_t MHPedestalPixel::SetupFill(const MParList *plist)
{
 
  Reset();

  return kTRUE;
}




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 = MakeDefCanvas(this,600,900); 
  
  c->Divide(1,2);
  
  c->cd(1);
  gPad->SetBorderMode(0);
  gPad->SetTicks();

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

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

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

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

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

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

  return;
}


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 area_guess  = entries/gkSq2Pi;
  const Double_t mu_guess    = fHPedestalCharge->GetBinCenter(fHPedestalCharge->GetMaximumBin());
  const Double_t sigma_guess = mu_guess/15.;

  TString name = TString("GausFit");
  name += fPixId;

  fGausFit = new TF1(name.Data(),"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) - 3.0*fGausFit->GetParameter(2);
      rmin        = (rtry < rmin ? rmin : rtry);
      rmax        = fGausFit->GetParameter(1) + 3.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 = 50;

  CutEdges(fHPedestalCharge,nbins);

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

  CutEdges(fHPedestalChargevsN,0);

}

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

  if (TESTBIT(fFlags,kFitted))
    {
      
      *fLog << all << "Results of the Summed Charges Fit: "                      << endl;
      *fLog << all << "Chisquare: "        << fChargeChisquare                   << endl;
      *fLog << all << "DoF: "              << fChargeNdf                         << endl;
      *fLog << all << "Probability: "      << fChargeProb                        << endl;
      *fLog << all << "Results of fit: "   ;
        if (TESTBIT(fFlags,kFitOK))
          *fLog << inf << "OK" << endl;
        else 
          *fLog << err << "NOT OK" << endl;          
      *fLog << all                                                               << endl;
    }
  else
    {
      *fLog << all << "Pedestal Histogram has not yet been fitted"               << endl;      
      *fLog << all                                                               << endl;      
    }
    
}