/* ======================================================================== *\
!
! *
! * 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 11/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

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

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

#include <TStyle.h>
#include <TCanvas.h>
#include <TPaveText.h>
#include <TText.h>
#include <TGraph.h>

#include "MFFT.h"

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

ClassImp(MHCalibrationPixel);

using namespace std;

const Int_t   MHCalibrationPixel::fChargeNbinsHiGain = 2100;
const Int_t   MHCalibrationPixel::fChargeNbinsLoGain = 1010;
const Int_t   MHCalibrationPixel::fAbsTimeNbins      = 16;
const Axis_t  MHCalibrationPixel::fAbsTimeFirst      = - 0.25;
const Axis_t  MHCalibrationPixel::fAbsTimeLast       =  15.75;
const Float_t MHCalibrationPixel::fProbLimit         = 0.001;
const Int_t   MHCalibrationPixel::fNDFLimit          = 5;

const Int_t   MHCalibrationPixel::fPSDNbins          = 30;
  
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationPixel::MHCalibrationPixel(const char *name, const char *title)
      : fPixId(-1),
	fHivsLoGain(NULL),
        fPSDHiGain(NULL),
        fPSDLoGain(NULL),
        fChargeGausFit(NULL),
        fHPSD(NULL),
        fPSDExpFit(NULL),
        fChargeXaxis(NULL),
        fPSDXaxis(NULL),
	fFitLegend(NULL),
        fCurrentSize(1024)
{ 

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

    fChargeFirstHiGain = -100.5;
    fChargeLastHiGain  = 1999.5;
    fChargeFirstLoGain = -100.5;
    fChargeLastLoGain  = 9999.5;

    // Create a large number of bins, later we will rebin
    fHChargeHiGain = new TH1F("HChargeHiGain","Distribution of Summed FADC Hi Gain Slices Pixel ",
			      fChargeNbinsHiGain,fChargeFirstHiGain,fChargeLastHiGain);
    fHChargeLoGain = new TH1F("HChargeLoGain","Distribution of Summed FADC Lo Gain Slices Pixel ",
			      fChargeNbinsLoGain,fChargeFirstLoGain,fChargeLastLoGain);

    fHChargeLoGain->SetXTitle("Sum FADC Slices (Lo Gain)");
    fHChargeHiGain->SetXTitle("Sum FADC Slices (Hi Gain)");

    fHChargeLoGain->SetYTitle("Nr. of events");
    fHChargeHiGain->SetYTitle("Nr. of events");

    fHChargeHiGain->Sumw2();
    fHChargeLoGain->Sumw2();

    // Absolute Times
    fHAbsTimeHiGain = new TH1F("HAbsTimeHiGain","Distribution of Absolute Arrival Hi Gain Times Pixel ",
                               fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast);
    fHAbsTimeLoGain = new TH1F("HAbsTimeLoGain","Distribution of Absolute Arrival Lo Gain Times Pixel ",
                               fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast);

    fHAbsTimeHiGain->SetXTitle("Absolute Arrival Time [Hi Gain FADC slice nr]");
    fHAbsTimeLoGain->SetXTitle("Absolute Arrival Time [Lo Gain FADC slice nr]");

    fHAbsTimeHiGain->SetYTitle("Nr. of events");
    fHAbsTimeLoGain->SetYTitle("Nr. of events");


    fHChargeHiGain->SetDirectory(NULL);
    fHChargeLoGain->SetDirectory(NULL);
    fHAbsTimeHiGain->SetDirectory(NULL);
    fHAbsTimeLoGain->SetDirectory(NULL);

    fHiGains = new TArrayF(fCurrentSize);
    fLoGains = new TArrayF(fCurrentSize);

    Clear();

}

MHCalibrationPixel::~MHCalibrationPixel()
{

  delete fHChargeHiGain;
  delete fHAbsTimeHiGain;

  delete fHChargeLoGain;
  delete fHAbsTimeLoGain;

  delete fHiGains;
  delete fLoGains;

  if (fChargeGausFit)
    delete fChargeGausFit;
  if (fPSDExpFit)
    delete fPSDExpFit;
  if (fHPSD)
    delete fHPSD;
  if (fFitLegend)
    delete fFitLegend;
  if (fHivsLoGain)
    delete fHivsLoGain;
  if (fChargeXaxis)
    delete fChargeXaxis;
  if (fPSDXaxis)
    delete fPSDXaxis;
}


void MHCalibrationPixel::Clear(Option_t *o)
{
  
  fTotalEntries            = 0;
  fCurrentSize             = 1024;

  fChargeFirstHiGain       = -100.5;
  fChargeLastHiGain        = 1999.5;
  fChargeFirstLoGain       = -100.5;
  fChargeLastLoGain        = 9999.5;

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

  fAbsTimeFirstHiGain     = -1.;
  fAbsTimeFirstLoGain     = -1.;
  fAbsTimeLastHiGain      = -1.;
  fAbsTimeLastLoGain      = -1.;
  
  fAbsTimeMean            = -1.;
  fAbsTimeMeanErr         = -1.;
  fAbsTimeRms             = -1.;

  fOffset = 0.;
  fSlope  = 0.;

  if (fChargeGausFit)
    delete fChargeGausFit;
  if (fPSDExpFit)
    delete fPSDExpFit;
  if (fHPSD)
    delete fHPSD;
  if (fFitLegend)
    delete fFitLegend;
  if (fHivsLoGain)
    delete fHivsLoGain;
  if (fChargeXaxis)
    delete fChargeXaxis;
  if (fPSDXaxis)
    delete fPSDXaxis;
  if (fPSDHiGain)
    delete fPSDHiGain;
  if (fPSDLoGain)
    delete fPSDLoGain;

  CLRBIT(fFlags,kUseLoGain);
  CLRBIT(fFlags,kChargeFitOK);
  CLRBIT(fFlags,kOscillating);

  return;
}


void MHCalibrationPixel::Reset()
{
  
  Clear();
  
  fHChargeHiGain->Reset();
  fHChargeLoGain->Reset();

  fHAbsTimeHiGain->Reset();
  fHAbsTimeLoGain->Reset();

  fHiGains->Set(1024);
  fLoGains->Set(1024);

  fHiGains->Reset(0.);
  fLoGains->Reset(0.);

}

void MHCalibrationPixel::SetUseLoGain(Bool_t b)
{
  if (b)
    SETBIT(fFlags, kUseLoGain) ;
  else
    CLRBIT(fFlags, kUseLoGain);
}


Bool_t MHCalibrationPixel::CheckOscillations()
{

  if (fPSDExpFit)
    return IsOscillating();

  //
  // The number of entries HAS to be a potence of 2, 
  // so we can only cut out from the last potence of 2 to the rest. 
  // Another possibility would be to fill everything with 
  // zeros, but that gives a low frequency peak, which we would 
  // have to cut out later again. 
  //
  // So, we have to live with the possibility that at the end 
  // of the calibration run, something has happened without noticing 
  // it...
  //
  
  // This cuts only the non-used zero's, but MFFT will later cut the rest
  CutArrayBorder(fHiGains);
  CutArrayBorder(fLoGains);  

  MFFT fourier;

  fPSDLoGain = fourier.PowerSpectrumDensity(fLoGains);
  fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains);

  if (IsUseLoGain())
    fHPSD = MH::ProjectArray(fPSDLoGain, fPSDNbins, 
                             Form("%s%d","PSDProjection",fPixId),
                             Form("%s%d","Power Spectrum Density Projection LoGain ",fPixId));
  else
    fHPSD = MH::ProjectArray(fPSDHiGain, fPSDNbins, 
                             Form("%s%d","PSDProjection",fPixId),
                             Form("%s%d","Power Spectrum Density Projection HiGain ",fPixId));

  //
  // First guesses for the fit (should be as close to reality as possible, 
  //
  const Double_t xmax = fHPSD->GetXaxis()->GetXmax();

  fPSDExpFit = new TF1(Form("%s%d","PSDExpFit",fPixId),"exp([0]-[1]*x)",0.,xmax);

  const Double_t slope_guess  = (TMath::Log(fHPSD->GetEntries())+1.)/xmax;
  const Double_t offset_guess = slope_guess*xmax;

  fPSDExpFit->SetParameters(offset_guess, slope_guess);
  fPSDExpFit->SetParNames("Offset","Slope");
  fPSDExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess);
  fPSDExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess);
  fPSDExpFit->SetRange(0.,xmax);

  fHPSD->Fit(fPSDExpFit,"RQL0");
  
  fPSDProb  = fPSDExpFit->GetProb();

  if (fPSDProb < gkProbLimit)
    {
      SETBIT(fFlags,kOscillating);
      return kTRUE;
    }
  
  CLRBIT(fFlags,kOscillating);

  return kFALSE;
}

void MHCalibrationPixel::CreatePSDXaxis(Int_t n)
{
  
  if (fPSDXaxis)
    return;

  fPSDXaxis = new TArrayF(n);

  for (Int_t i=0;i<n;i++)
    fPSDXaxis->AddAt((Float_t)i,i);
}

void MHCalibrationPixel::CreateChargeXaxis(Int_t n)
{
  
  if (!fChargeXaxis)
    {
      fChargeXaxis = new TArrayF(n);
      for (Int_t i=0;i<n;i++)
        fChargeXaxis->AddAt((Float_t)i,i);
      return;
    }

  if (fChargeXaxis->GetSize() == n)
    return;

  const Int_t diff = fChargeXaxis->GetSize()-n;

  fChargeXaxis->Set(n);
  if (diff < 0)
    for (Int_t i=n-1;i>n-diff-1;i--)
      fChargeXaxis->AddAt((Float_t)i,i);
}

void MHCalibrationPixel::CutArrayBorder(TArrayF *array)
{
  
  Int_t i;

  for (i=array->GetSize()-1;i>=0;i--)
    if (array->At(i) != 0)
      {
        array->Set(i+1);
        break;
      }
}



Bool_t MHCalibrationPixel::IsEmpty() const
{
    return !(fHChargeHiGain->GetEntries() || fHChargeLoGain->GetEntries());
}

Bool_t MHCalibrationPixel::IsUseLoGain() const
{
  return TESTBIT(fFlags,kUseLoGain);
}

Bool_t MHCalibrationPixel::IsChargeFitOK() const 
{
    return TESTBIT(fFlags,kChargeFitOK);
}

Bool_t MHCalibrationPixel::IsOscillating()
{

  if (fPSDExpFit)
    return TESTBIT(fFlags,kOscillating);

  return CheckOscillations();

}

Bool_t MHCalibrationPixel::FillChargeLoGain(Float_t q)
{
    return (fHChargeLoGain->Fill(q) > -1);
}

Bool_t MHCalibrationPixel::FillAbsTimeLoGain(Float_t t)
{
    return (fHAbsTimeLoGain->Fill(t)> -1);
}

Bool_t MHCalibrationPixel::FillChargeHiGain(Float_t q)
{
    return (fHChargeHiGain->Fill(q) > -1);
}

Bool_t MHCalibrationPixel::FillAbsTimeHiGain(Float_t t)
{
    return (fHAbsTimeHiGain->Fill(t) > -1);
}

void MHCalibrationPixel::ChangeHistId(Int_t id)
{

  // Change only if the names have not yet been changed
  if (fPixId == -1)
    {
      
      //
      // Names Hi gain Histograms
      //
      fHChargeHiGain->SetName(   Form("%s%d",fHChargeHiGain->GetName(),   id));
      fHAbsTimeHiGain->SetName(  Form("%s%d",fHAbsTimeHiGain->GetName(),  id));
      
      //
      // Title Hi gain Histograms
      //
      fHChargeHiGain->SetTitle(    Form("%s%d",fHChargeHiGain->GetTitle(),   id));
      fHAbsTimeHiGain->SetTitle(   Form("%s%d",fHAbsTimeHiGain->GetTitle(),  id));
  
      //
      // Names Low Gain Histograms
      //
      fHChargeLoGain->SetName( Form("%s%d",fHChargeLoGain->GetName(),id));
      fHAbsTimeLoGain->SetName( Form("%s%d",fHAbsTimeLoGain->GetName(),id));

      //
      // Titles Low Gain Histograms
      //
      fHChargeLoGain->SetTitle( Form("%s%d",fHChargeLoGain->GetTitle(),id));
      fHAbsTimeLoGain->SetTitle( Form("%s%d",fHAbsTimeLoGain->GetTitle(),id));
      
      fPixId = id;
    }
  
}


Bool_t MHCalibrationPixel::UseLoGain()
{

  if (fHChargeHiGain->Integral() > fHChargeLoGain->Integral())
    {
      CLRBIT(fFlags,kUseLoGain);
      return kFALSE;
    }
  else
    {
      SETBIT(fFlags,kUseLoGain);
      return kTRUE;
    }
}

Bool_t MHCalibrationPixel::FillGraphs(Float_t qhi,Float_t qlo)
{

  if (fTotalEntries >= fCurrentSize)
    {
      fCurrentSize *= 2;

      fHiGains->Set(fCurrentSize);
      fLoGains->Set(fCurrentSize);
    }
  
  fHiGains->AddAt(qhi,fTotalEntries);
  fLoGains->AddAt(qlo,fTotalEntries);

  fTotalEntries++;

  return kTRUE;

}


// -------------------------------------------------------------------------
//
// Draw a legend with the fit results
//
void MHCalibrationPixel::DrawLegend()
{

  fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);

  if (IsChargeFitOK()) 
      fFitLegend->SetFillColor(80);
  else
      fFitLegend->SetFillColor(2);    
  
  fFitLegend->SetLabel("Results of the Gauss Fit:");
  fFitLegend->SetTextSize(0.05);

  const TString line1 = 
    Form("Mean: Q_{#mu} = %2.2f #pm %2.2f",fChargeMean,fChargeMeanErr);
  TText *t1 = fFitLegend->AddText(line1);
  t1->SetBit(kCanDelete);

  const TString line4 =
    Form("Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fChargeSigma,fChargeSigmaErr);
  TText *t2 =  fFitLegend->AddText(line4);
  t2->SetBit(kCanDelete);

  const TString line7 =
    Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChargeChisquare,fChargeNdf);
  TText *t3 =  fFitLegend->AddText(line7);
  t3->SetBit(kCanDelete);

  const TString line8 =
    Form("Probability: %4.3f ",fChargeProb);
  TText *t4 =  fFitLegend->AddText(line8);
  t4->SetBit(kCanDelete);

  if (IsChargeFitOK())
    {
      TText *t5 = fFitLegend->AddText("Result of the Fit: OK");
      t5->SetBit(kCanDelete);
    }
  else
    {
      TText *t6 = fFitLegend->AddText("Result of the Fit: NOT OK");
      t6->SetBit(kCanDelete);
    }
  
  fFitLegend->SetBit(kCanDelete);
  fFitLegend->Draw();

}


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

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

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

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


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

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

  if (fHChargeHiGain->Integral() > 0)
    gPad->SetLogy();

  fHChargeHiGain->Draw(opt);
  
  if (IsUseLoGain())
    {

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

      if (fHChargeLoGain->Integral() > 0)
        gPad->SetLogy();

      fHChargeLoGain->Draw(opt);
      
      if (fChargeGausFit)
        {
          fChargeGausFit->SetLineColor(IsChargeFitOK() ? kGreen : kRed);          
          fChargeGausFit->Draw("same");
        }

      c->cd(3);
      gROOT->SetSelectedPad(NULL);
      gStyle->SetOptFit();
      if (fHivsLoGain)
        fHivsLoGain->Draw("prof");

      c->cd(4);
      DrawLegend();
    }
  else
    {
      if (fChargeGausFit)
        {
          fChargeGausFit->SetLineColor(IsChargeFitOK() ? kGreen : kRed);          
          fChargeGausFit->Draw("same");
        }

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

      if (fHChargeLoGain->Integral() > 0)
        gPad->SetLogy(1);
      
      fHChargeLoGain->Draw(opt);

      c->cd(3);
      DrawLegend();
      
      c->cd(4);
      
      gROOT->SetSelectedPad(NULL);
      gStyle->SetOptFit();
      if (fHivsLoGain)
        fHivsLoGain->Draw("prof");
    }	
  
  c->cd(5);
  gPad->SetTicks();
  fHAbsTimeHiGain->Draw(opt);

  c->cd(6);
  gPad->SetTicks();
  fHAbsTimeLoGain->Draw(opt);

  CreateChargeXaxis(fHiGains->GetSize());

  c->cd(7);
  gPad->SetTicks();
  TGraph *gr1 = new TGraph(fChargeXaxis->GetSize(),
                           fChargeXaxis->GetArray(),
                           fHiGains->GetArray());  
  gr1->SetTitle("Evolution of HiGain charges with event number");
  gr1->SetBit(kCanDelete);
  gr1->Draw("AL");
  
  CreateChargeXaxis(fLoGains->GetSize());

  c->cd(8);
  gPad->SetTicks();
  TGraph *gr2 = new TGraph(fChargeXaxis->GetSize(),
                           fChargeXaxis->GetArray(),
                           fLoGains->GetArray());  
  gr2->SetTitle("Evolution of LoGain charges with event number");
  gr2->SetBit(kCanDelete);
  gr2->Draw("AL");
  
  c->Modified();
  c->Update();
  
  c->cd(9);
  
  TArrayF *array;
  TString title = "Power Spectrum Density ";

  if(IsUseLoGain())
    {
      if (!fPSDLoGain)
        return;
      array = fPSDLoGain;
      title += "LoGain";
    }
  else
    {
      if (!fPSDHiGain)
        return;
      array = fPSDHiGain;
      title += "HiGain";
    }
  
  if (!fPSDXaxis)
    CreatePSDXaxis(array->GetSize());

  TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray());
  gr3->SetTitle(title.Data());
  gr3->SetBit(kCanDelete);
  gr3->Draw("AL");

  c->Modified();
  c->Update();
  
  c->cd(10);
  
  gStyle->SetOptStat(111111);
  gStyle->SetOptFit();
  gPad->SetTicks();  

  if (fHPSD->Integral() > 0)
    gPad->SetLogy();

  fHPSD->Draw(opt);

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

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

  return;
}

void MHCalibrationPixel::FitHiGainvsLoGain()
{

  *fLog << inf << "Fit Hi Gain vs Lo Gain " << endl;

  if (fHivsLoGain)
    return;

  if ((!fHiGains) || (!fLoGains) || (fHiGains->GetSize() == 0) || (fLoGains->GetSize() == 0))
    return;

  gStyle->SetOptFit();

  fHivsLoGain = new TProfile("HiGainvsLoGain","Plot the High Gain vs. Low Gain",100,0.,1000.,0.,1000.);
  fHivsLoGain->GetXaxis()->SetTitle("Sum of Charges High Gain");
  fHivsLoGain->GetYaxis()->SetTitle("Sum of Charges Low Gain");
  fHivsLoGain->SetName(Form("%s%d",fHivsLoGain->GetName(),fPixId));
  
  for (Int_t i=0;i<fTotalEntries;i++)
    fHivsLoGain->Fill(fHiGains->At(i),fLoGains->At(i),1);
  
  fHivsLoGain->Fit("pol1","rq","",fHivsLoGain->GetMean()-2.5*fHivsLoGain->GetRMS(),
                                  fHivsLoGain->GetMean()+2.5*fHivsLoGain->GetRMS());
  
  fOffset = fHivsLoGain->GetFunction("pol1")->GetParameter(0);
  fSlope  = fHivsLoGain->GetFunction("pol1")->GetParameter(1);
  
}



Bool_t MHCalibrationPixel::FitCharge(Option_t *option)
{

  if (fChargeGausFit)
    return kFALSE;

  //
  // Get the fitting ranges
  //
  Axis_t rmin = fChargeFirstHiGain;
  Axis_t rmax = fChargeLastHiGain;
  TH1F *hist = fHChargeHiGain;

  if (TESTBIT(fFlags,kUseLoGain))
    {
      rmin = fChargeFirstLoGain;
      rmax = fChargeLastLoGain;
      hist = fHChargeLoGain;
    }
  
  //
  // First guesses for the fit (should be as close to reality as possible, 
  //
  const Stat_t   entries     = hist->Integral("width");
  const Double_t mu_guess    = hist->GetBinCenter(hist->GetMaximumBin());
  const Double_t sigma_guess = (rmax-rmin)/2.;
  const Double_t area_guess  = entries/gkSq2Pi/sigma_guess;

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

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

  hist->Fit(fChargeGausFit,option);
  
  // 
  // If we are not able to fit, try once again
  //
  if (fChargeGausFit->GetProb() < fProbLimit)
    {

      Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2);
      rmin        = (rtry < rmin ? rmin : rtry);
      rmax        = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2);

      fChargeGausFit->SetRange(rmin,rmax);  
      hist->Fit(fChargeGausFit,option);
    }
  
  fChargeChisquare = fChargeGausFit->GetChisquare();
  fChargeNdf       = fChargeGausFit->GetNDF();
  fChargeProb      = fChargeGausFit->GetProb();
  fChargeMean      = fChargeGausFit->GetParameter(1);
  fChargeMeanErr   = fChargeGausFit->GetParError(1);
  fChargeSigma     = fChargeGausFit->GetParameter(2);
  fChargeSigmaErr  = fChargeGausFit->GetParError(2);

  //
  // From the absolute time, we only take the mean and RMS
  //
  fAbsTimeMean    = (Float_t)fHAbsTimeHiGain->GetMean();
  fAbsTimeRms     = (Float_t)fHAbsTimeHiGain->GetRMS();
  fAbsTimeMeanErr = (Float_t)fAbsTimeRms / TMath::Sqrt(fHAbsTimeHiGain->GetEntries());
  
  if (TESTBIT(fFlags,kUseLoGain))
    {
      fAbsTimeMean    = fHAbsTimeLoGain->GetMean();
      fAbsTimeRms     = fHAbsTimeLoGain->GetRMS();
      fAbsTimeMeanErr = fAbsTimeRms / TMath::Sqrt(fHAbsTimeLoGain->GetEntries());
    }

  //
  // 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%)
  //
  // Otherwise take means and RMS of the histograms
  //
  if (    TMath::IsNaN(fChargeMean) 
      || TMath::IsNaN(fChargeMeanErr)
      || TMath::IsNaN(fChargeProb)    
      || TMath::IsNaN(fChargeSigma)
      || TMath::IsNaN(fChargeSigmaErr)
      || (fChargeNdf < fNDFLimit)
      || (fChargeProb < fProbLimit) )
    {

      fChargeMean      = hist->GetMean();
      fChargeMeanErr   = hist->GetRMS()/hist->GetEntries();
      fChargeSigma     = hist->GetRMS();
      fChargeSigmaErr  = fChargeMeanErr/2.;
      
      CLRBIT(fFlags,kChargeFitOK);
      return kFALSE;
    }
  
  SETBIT(fFlags,kChargeFitOK);
  return kTRUE;
}


void MHCalibrationPixel::CutAllEdges()
{

  Int_t nbins = 30;

  CutEdges(fHChargeHiGain,nbins);

  fChargeFirstHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetFirst());
  fChargeLastHiGain  = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetLast())
                      +fHChargeHiGain->GetBinWidth(0);

  CutEdges(fHChargeLoGain,nbins);

  fChargeFirstLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetFirst());
  fChargeLastLoGain  = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetLast())
                      +fHChargeLoGain->GetBinWidth(0);


  CutEdges(fHAbsTimeHiGain,0);

  fAbsTimeFirstHiGain = fHAbsTimeHiGain->GetBinLowEdge(fHAbsTimeHiGain->GetXaxis()->GetFirst());
  fAbsTimeLastHiGain  = fHAbsTimeHiGain->GetBinLowEdge(fHAbsTimeHiGain->GetXaxis()->GetLast())
                       +fHAbsTimeHiGain->GetBinWidth(0);

  CutEdges(fHAbsTimeLoGain,0);

  fAbsTimeFirstLoGain = fHAbsTimeLoGain->GetBinLowEdge(fHAbsTimeLoGain->GetXaxis()->GetFirst());
  fAbsTimeLastLoGain  = fHAbsTimeLoGain->GetBinLowEdge(fHAbsTimeLoGain->GetXaxis()->GetLast())
                       +fHAbsTimeLoGain->GetBinWidth(0);


}

void MHCalibrationPixel::PrintChargeFitResult()
{
  
  *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                                                               << endl;
  
}