/* ======================================================================== *\
!
! *
! * 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 "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::fChargevsNbins     = 5000;
const Int_t   MHCalibrationPixel::fAbsTimeNbins      = 32;
const Axis_t  MHCalibrationPixel::fAbsTimeFirst      = - 0.25;
const Axis_t  MHCalibrationPixel::fAbsTimeLast       =  15.75;
const Int_t   MHCalibrationPixel::fRelTimeNbins      = 900;
const Axis_t  MHCalibrationPixel::fRelTimeFirst      = -13.;
const Axis_t  MHCalibrationPixel::fRelTimeLast       =  13.;
const Float_t MHCalibrationPixel::fProbLimit         = 0.001;
const Int_t   MHCalibrationPixel::fNDFLimit          = 5;

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationPixel::MHCalibrationPixel(const char *name, const char *title)
      : fPixId(-1),
	fHivsLoGain(NULL),
        fHPSD(NULL),
        fChargeGausFit(NULL),
	fRelTimeGausFit(NULL), 
	fFitLegend(NULL)
{ 

    fName  = name  ? name  : "MHCalibrationPixel";
    fTitle = title ? title : "Fill the accumulated charges and times 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");

    // Relative Times
    fHRelTimeHiGain = new TH1F("HRelTimeHiGain","Distribution of Relative Arrival Times High Gain Pixel ",
                               fRelTimeNbins,fRelTimeFirst,fRelTimeLast);
    fHRelTimeLoGain = new TH1F("HRelTimeLoGain","Distribution of Relative Arrival Time Low Gain Pixel ",
                               fRelTimeNbins,fRelTimeFirst,fRelTimeLast);

    fHRelTimeHiGain->SetXTitle("Relative Arrival Times [Hi Gain FADC slice nr]");
    fHRelTimeLoGain->SetXTitle("Relative Arrival Times [Lo Gain FADC slice nr]");

    fHRelTimeHiGain->SetYTitle("Nr. of events");
    fHRelTimeLoGain->SetYTitle("Nr. of events");

    // We define a reasonable number and later enlarge it if necessary
    fHChargevsNHiGain = new TH1I("HChargevsNHiGain","Sum of Hi Gain Charges vs. Event Number Pixel ",
				 fChargevsNbins,-0.5,(Axis_t)fChargevsNbins - 0.5);
    fHChargevsNLoGain = new TH1I("HChargevsNLoGain","Sum of Lo Gain Charges vs. Event Number Pixel ",
				 fChargevsNbins,-0.5,(Axis_t)fChargevsNbins - 0.5);

    fHChargevsNHiGain->SetXTitle("Event Nr.");
    fHChargevsNLoGain->SetXTitle("Event Nr.");

    fHChargevsNHiGain->SetYTitle("Sum of Hi Gain FADC slices");
    fHChargevsNLoGain->SetYTitle("Sum of Lo Gain FADC slices");

    fHChargeHiGain->SetDirectory(NULL);
    fHChargeLoGain->SetDirectory(NULL);
    fHAbsTimeHiGain->SetDirectory(NULL);
    fHAbsTimeLoGain->SetDirectory(NULL);
    fHRelTimeHiGain->SetDirectory(NULL);
    fHRelTimeLoGain->SetDirectory(NULL);
    fHChargevsNHiGain->SetDirectory(NULL);
    fHChargevsNLoGain->SetDirectory(NULL);

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

    Clear();

}

MHCalibrationPixel::~MHCalibrationPixel()
{

  delete fHChargeHiGain;
  delete fHAbsTimeHiGain;
  delete fHRelTimeHiGain;
  delete fHChargevsNHiGain;

  delete fHChargeLoGain;
  delete fHAbsTimeLoGain;
  delete fHRelTimeLoGain;
  delete fHChargevsNLoGain;

  delete fHiGains;
  delete fLoGains;

  if (fChargeGausFit)
    delete fChargeGausFit;
  if (fRelTimeGausFit)
    delete fRelTimeGausFit;
  if (fFitLegend)
    delete fFitLegend;
  if (fHivsLoGain)
    delete fHivsLoGain;
}


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

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

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

  fRelTimeChisquare        = -1.;
  fRelTimeProb             = -1.;
  fRelTimeNdf              = -1;
  fRelTimeMean             = -1.;
  fRelTimeSigma            = -1.;

  fRelTimeLowerFitRangeHiGain = -99.;
  fRelTimeUpperFitRangeHiGain = -99.;
  fRelTimeLowerFitRangeLoGain = -99.;
  fRelTimeUpperFitRangeLoGain = -99.;

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

  fOffset = 0.;
  fSlope  = 0.;

  if (fChargeGausFit)
    delete fChargeGausFit;
  if (fRelTimeGausFit)
    delete fRelTimeGausFit;
  if (fFitLegend)
    delete fFitLegend;
  if (fHivsLoGain)
    delete fHivsLoGain;
  if (fHPSD)
    delete fHPSD;

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

  return;
}


void MHCalibrationPixel::Reset()
{
  
  Clear();
  
  fHChargeHiGain->Reset();
  fHChargeLoGain->Reset();
  fHAbsTimeHiGain->Reset();
  fHAbsTimeLoGain->Reset();
  fHRelTimeHiGain->Reset();
  fHRelTimeLoGain->Reset();
  fHChargevsNHiGain->Reset();
  fHChargevsNLoGain->Reset();

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

}

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


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::IsTimeFitOK() const 
{
    return TESTBIT(fFlags,kTimeFitOK);
}

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::FillRelTimeLoGain(Float_t t)
{
    return (fHRelTimeLoGain->Fill(t) > -1);
}

Bool_t MHCalibrationPixel::FillChargevsNLoGain(Float_t q, Int_t n)
{
    return (fHChargevsNLoGain->Fill(n,q) > -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);
}

Bool_t MHCalibrationPixel::FillRelTimeHiGain(Float_t t)
{
    return (fHRelTimeHiGain->Fill(t) > -1);
}

Bool_t MHCalibrationPixel::FillChargevsNHiGain(Float_t q, Int_t n)
{
    return (fHChargevsNHiGain->Fill(n,q) > -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));
      fHRelTimeHiGain->SetName(  Form("%s%d",fHRelTimeHiGain->GetName(),  id));
      fHChargevsNHiGain->SetName(Form("%s%d",fHChargevsNHiGain->GetName(),id));
      
      //
      // Title Hi gain Histograms
      //
      fHChargeHiGain->SetTitle(    Form("%s%d",fHChargeHiGain->GetTitle(),   id));
      fHAbsTimeHiGain->SetTitle(   Form("%s%d",fHAbsTimeHiGain->GetTitle(),  id));
      fHRelTimeHiGain->SetTitle(   Form("%s%d",fHRelTimeHiGain->GetTitle(),  id));
      fHChargevsNHiGain->SetTitle( Form("%s%d",fHChargevsNHiGain->GetTitle(),id));
  
      //
      // Names Low Gain Histograms
      //
      fHChargeLoGain->SetName( Form("%s%d",fHChargeLoGain->GetName(),id));
      fHAbsTimeLoGain->SetName( Form("%s%d",fHAbsTimeLoGain->GetName(),id));
      fHRelTimeLoGain->SetName( Form("%s%d",fHRelTimeLoGain->GetName(),id));
      fHChargevsNLoGain->SetName( Form("%s%d",fHChargevsNLoGain->GetName(),id));

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


Bool_t MHCalibrationPixel::SetupFill(const MParList *plist)
{
 
  Reset();
  
  return kTRUE;
}


Bool_t MHCalibrationPixel::UseLoGain()
{

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

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

  fHiGains->Set(++fTotalEntries);
  fLoGains->Set(fTotalEntries);

  fHiGains->AddAt(qhi,fTotalEntries-1);
  fLoGains->AddAt(qlo,fTotalEntries-1);

  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->Divide(2,5);
  
  c->cd(1);
  gPad->SetBorderMode(0);
  gPad->SetTicks();

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

  fHChargeHiGain->Draw(opt);
  
  c->Modified();
  c->Update();
  
  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->Modified();
      c->Update();
      
      c->cd(3);
      gROOT->SetSelectedPad(NULL);
      gStyle->SetOptFit();
      if (fHivsLoGain)
        fHivsLoGain->Draw("prof");
      gPad->Modified();
      gPad->Update();
      
      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->Modified();
      c->Update();
      
      c->cd(3);
      DrawLegend();
      
      c->cd(4);
      
      gROOT->SetSelectedPad(NULL);
      gStyle->SetOptFit();
      if (fHivsLoGain)
        fHivsLoGain->Draw("prof");
      gPad->Modified();
      gPad->Update();
    }	
  
  c->Modified();
  c->Update();
  
  c->cd(5);
  gStyle->SetOptStat(1111111);

  gPad->SetTicks();  
  gPad->SetLogy(0);
  fHRelTimeHiGain->Draw(opt);
  c->Modified();
  c->Update();
  
  if (IsUseLoGain())
      {
        
	c->cd(6);
        gPad->SetTicks();
	gPad->SetLogy(0);
	fHRelTimeLoGain->Draw(opt);
	c->Modified();
	c->Update();
	
	if (fRelTimeGausFit)
	  {
            fRelTimeGausFit->SetLineColor(IsTimeFitOK() ? kGreen : kRed);
	    fRelTimeGausFit->Draw("same");
	  }

        c->Modified();
        c->Update();
      }
    else
      {  
	if (fRelTimeGausFit)
	  {
            fRelTimeGausFit->SetLineColor(IsTimeFitOK() ? kGreen : kRed);
	    fRelTimeGausFit->Draw("same");
	  }
        
        c->Modified();
        c->Update();
	
        c->cd(6);
        gPad->SetTicks();
	gPad->SetLogy(0);
	fHRelTimeLoGain->Draw(opt);	
	c->Modified();
	c->Update();
    
      }
    c->Modified();
    c->Update();
    
    c->cd(7);
    gPad->SetTicks();
    fHAbsTimeHiGain->Draw(opt);
    c->Modified();
    c->Update();

    c->cd(8);
    gPad->SetTicks();
    fHAbsTimeLoGain->Draw(opt);
    c->Modified();
    c->Update();


    c->cd(9);
    gPad->SetTicks();
    fHChargevsNHiGain->Draw(opt);
    c->Modified();
    c->Update();

    c->cd(10);
    gPad->SetTicks();
    fHChargevsNLoGain->Draw(opt);
    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");

  TString titleHiLoGain  = TString(fHivsLoGain->GetName());
  titleHiLoGain += fPixId;
  fHivsLoGain->SetName(titleHiLoGain.Data());
  
  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::FitTime(Option_t *option)
{

  if (fRelTimeGausFit)
    return kFALSE;

  //
  // 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());
    }
  
  //
  // Get the fitting ranges
  //
  Axis_t rmin = fRelTimeLowerFitRangeHiGain;
  Axis_t rmax = fRelTimeUpperFitRangeHiGain;
  TH1F *hist = fHRelTimeHiGain;

  if (TESTBIT(fFlags,kUseLoGain))
    {
      rmin = fRelTimeLowerFitRangeLoGain;
      rmax = fRelTimeUpperFitRangeLoGain;
      hist = fHRelTimeLoGain;
    }
  
  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  = 2.*entries/gkSq2Pi/sigma_guess;

  TString name = TString("GausRelTime");
  name += fPixId;
  fRelTimeGausFit = new TF1(name.Data(),"gaus",rmin,rmax);  

  if (!fRelTimeGausFit) 
    {
    *fLog << warn << dbginf << "WARNING: Could not create fit function for RelTime fit" << endl;
    return kFALSE;
    }

  fRelTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
  fRelTimeGausFit->SetParNames("Area","#mu","#sigma");
  fRelTimeGausFit->SetParLimits(0,0.,5.*area_guess);
  fRelTimeGausFit->SetParLimits(1,rmin,rmax);
  fRelTimeGausFit->SetParLimits(2,0.,(rmax-rmin));
  fRelTimeGausFit->SetRange(rmin,rmax);

  hist->Fit(fRelTimeGausFit,option);

  //
  // If the fit does not converge, try another one with smaller bounderies
  //
  if (fRelTimeGausFit->GetProb() < 0.001)
    {
      rmin = fRelTimeGausFit->GetParameter(1) - 1.5*fRelTimeGausFit->GetParameter(2);
      rmax = fRelTimeGausFit->GetParameter(1) + 1.5*fRelTimeGausFit->GetParameter(2);
      fRelTimeGausFit->SetRange(rmin,rmax);  
      hist->Fit(fRelTimeGausFit,option);
    }
  
  fRelTimeChisquare = fRelTimeGausFit->GetChisquare();
  fRelTimeNdf       = fRelTimeGausFit->GetNDF();
  fRelTimeProb      = fRelTimeGausFit->GetProb();

  fRelTimeMean      = fRelTimeGausFit->GetParameter(1);
  fRelTimeSigma     = fRelTimeGausFit->GetParameter(2);

  fRelTimeMeanErr   = fRelTimeGausFit->GetParError(1);

  if (TMath::IsNaN(fRelTimeMean) || TMath::IsNaN(fRelTimeSigma))
    {
      CLRBIT(fFlags,kTimeFitOK);
      return kFALSE;
    }
  
  if (TMath::IsNaN(fRelTimeChisquare) || (fRelTimeProb < fProbLimit))
    {
      CLRBIT(fFlags,kTimeFitOK);
      return kFALSE;
    }
  
  SETBIT(fFlags,kTimeFitOK);

  return kTRUE;

}


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

  //
  // 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 = 50;

  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(fHRelTimeHiGain,0);

  fRelTimeLowerFitRangeHiGain = fHRelTimeHiGain->GetBinLowEdge(fHRelTimeHiGain->GetXaxis()->GetFirst());
  fRelTimeUpperFitRangeHiGain = fHRelTimeHiGain->GetBinLowEdge(fHRelTimeHiGain->GetXaxis()->GetLast())
                               +fHRelTimeHiGain->GetBinWidth(0);

  CutEdges(fHRelTimeLoGain,0);

  fRelTimeLowerFitRangeLoGain = fHRelTimeLoGain->GetBinLowEdge(fHRelTimeLoGain->GetXaxis()->GetFirst());
  fRelTimeUpperFitRangeLoGain = fHRelTimeLoGain->GetBinLowEdge(fHRelTimeLoGain->GetXaxis()->GetLast())
                               +fHRelTimeLoGain->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);

  CutEdges(fHChargevsNHiGain,0);
  CutEdges(fHChargevsNLoGain,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;
  
}

void MHCalibrationPixel::PrintTimeFitResult()
{

  *fLog << all << "Results of the Time Slices Fit: "                        << endl;
  *fLog << all << "Chisquare: "   << fRelTimeChisquare                      << endl;
  *fLog << all << "Ndf: "         << fRelTimeNdf                            << endl;
  *fLog << all << "Probability: " << fRelTimeProb                           << endl;
  *fLog << all                                                              << endl;

}