source: trunk/MagicSoft/Mars/manalysis/MHExtractedSignalPix.cc@ 3038

/* ======================================================================== *\
!
! *
! * 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>
!              Thomas Bretz 02/2004 <mailto:bretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MHExtractedSignalPix
//
//////////////////////////////////////////////////////////////////////////////
#include "MHExtractedSignalPix.h"

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

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

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

#include "MH.h"

ClassImp(MHExtractedSignalPix);

using namespace std;

const Float_t gkSq2Pi     = TMath::Sqrt(TMath::TwoPi());
const Float_t gkProbLimit = 0.01;
const Int_t   MHExtractedSignalPix::gkChargeNbins    = 500 ;
const Int_t   MHExtractedSignalPix::gkChargevsNbins  = 1000;
const Axis_t  MHExtractedSignalPix::gkChargevsNFirst = -0.5;
const Axis_t  MHExtractedSignalPix::gkChargevsNLast  = 999.5;
const Axis_t  MHExtractedSignalPix::gkChargeFirst    = -0.5;
const Axis_t  MHExtractedSignalPix::gkChargeLast     = 99.5;

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHExtractedSignalPix::MHExtractedSignalPix()
      : fPixId(-1),
        fGausFit(NULL)
{ 

    //fName  = name  ? name  : "MHExtractedSignalPix";
    //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 ",
                                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();
}

MHExtractedSignalPix::~MHExtractedSignalPix()
{

  delete fHPedestalCharge;
  delete fHPedestalChargevsN;

  if (fGausFit)
    delete fGausFit;
}


void MHExtractedSignalPix::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;

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


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


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

Bool_t MHExtractedSignalPix::IsFitted() const 
{
  return TESTBIT(fFlags,kFitted);
}

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


Bool_t MHExtractedSignalPix::IsOscillating() const 
{
  return TESTBIT(fFlags,kOscillating);
}


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

Bool_t MHExtractedSignalPix::FillChargevsN(Float_t q)
{
    fTotalEntries++;
    return (fHPedestalChargevsN->Fill(fTotalEntries,q) > -1);
}

void MHExtractedSignalPix::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 *MHExtractedSignalPix::DrawClone(Option_t *option) const
{

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

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


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

// -------------------------------------------------------------------------
//
// Draw the histogram
//
void MHExtractedSignalPix::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();

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


Bool_t MHExtractedSignalPix::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("width");
  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 MHExtractedSignalPix::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 MHExtractedSignalPix::Print(const Option_t *o) const
{

  gLog << 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;
}