/* ======================================================================== *\
!
! *
! * 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 "MCalibrationFits.h"

#include <TStyle.h>
#include <TMath.h>

#include <TFitter.h>
#include <TGraph.h>
#include <TAxis.h>

#include <TF1.h>
#include <TH2.h>
#include <TProfile.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TPaveText.h>

#include "MParList.h"

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

ClassImp(MHCalibrationPixel);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationPixel::MHCalibrationPixel(const char *name, const char *title)
      : fPixId(-1),
	fTotalEntries(0),
	fHiGains(NULL),
	fLoGains(NULL),
        fChargeGausFit(NULL), 
	fTimeGausFit(NULL), 
	fHivsLoGain(NULL),
	fFitLegend(NULL),
	fLowerFitRange(-2000.),
	fChargeFirstHiGain(-2000.5),
	fChargeLastHiGain(9999.5),
	fChargeNbinsHiGain(12000),
	fChargeFirstLoGain(-2000.5),
	fChargeLastLoGain(9999.5),
	fChargeNbinsLoGain(1200),
        fFitOK(kFALSE),
        fChargeChisquare(-1.),
        fChargeProb(-1.),
        fChargeNdf(-1),
        fTimeChisquare(-1.),
        fTimeProb(-1.),
        fTimeNdf(-1),
        fTimeMean(-1.),
        fTimeSigma(-1.),
        fTimeLowerFitRange(3.),
        fTimeUpperFitRange(9.),
	fUseLoGain(kFALSE),
	fOffset(0.),
	fSlope(0.)
{ 

    fName  = name  ? name  : "MHCalibrationPixel";
    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
    fHChargeHiGain = new TH1F("HChargeHiGain","Distribution of Summed FADC Hi Gain Slices Pixel ",
			      fChargeNbinsHiGain,fChargeFirstHiGain,fChargeLastHiGain);
    fHChargeHiGain->SetXTitle("Sum FADC Slices (Hi Gain)");
    fHChargeHiGain->SetYTitle("Nr. of events");
    fHChargeHiGain->Sumw2();

    fHChargeHiGain->SetDirectory(NULL);

    fHChargeLoGain = new TH1F("HChargeLoGain","Distribution of Summed FADC Lo Gain Slices Pixel ",
			      fChargeNbinsLoGain,fChargeFirstLoGain,fChargeLastLoGain);
    fHChargeLoGain->SetXTitle("Sum FADC Slices (Lo Gain)");
    fHChargeLoGain->SetYTitle("Nr. of events");
    fHChargeLoGain->Sumw2();

    fHChargeLoGain->SetDirectory(NULL);

    Axis_t tfirst = -0.5;
    Axis_t tlast  = 15.5;
    Int_t  ntbins = 16;

    fHTimeHiGain = new TH1I("HTimeHiGain","Distribution of Mean Arrival Hi Gain Times Pixel ",
			    ntbins,tfirst,tlast);
    fHTimeHiGain->SetXTitle("Mean Arrival Times [Hi Gain FADC slice nr]");
    fHTimeHiGain->SetYTitle("Nr. of events");
    // fHTimeHiGain->Sumw2();

    fHTimeHiGain->SetDirectory(NULL);

    fHTimeLoGain = new TH1I("HTimeLoGain","Distribution of Mean Arrival Lo Gain Times Pixel ",
			    ntbins,tfirst,tlast);
    fHTimeLoGain->SetXTitle("Mean Arrival Times [Lo Gain FADC slice nr]");
    fHTimeLoGain->SetYTitle("Nr. of events");
    //    fHTimeLoGain->Sumw2();

    fHTimeLoGain->SetDirectory(NULL);

    // We define a reasonable number and later enlarge it if necessary
    Int_t  nqbins = 20000;
    Axis_t nfirst = -0.5;
    Axis_t nlast  = (Axis_t)nqbins - 0.5;

    fHChargevsNHiGain = new TH1I("HChargevsNHiGain","Sum of Hi Gain Charges vs. Event Number Pixel ",
				 nqbins,nfirst,nlast);
    fHChargevsNHiGain->SetXTitle("Event Nr.");
    fHChargevsNHiGain->SetYTitle("Sum of Hi Gain FADC slices");

    fHChargevsNHiGain->SetDirectory(NULL);

    fHChargevsNLoGain = new TH1I("HChargevsNLoGain","Sum of Lo Gain Charges vs. Event Number Pixel ",
				 nqbins,nfirst,nlast);
    fHChargevsNLoGain->SetXTitle("Event Nr.");
    fHChargevsNLoGain->SetYTitle("Sum of Lo Gain FADC slices");

    fHChargevsNLoGain->SetDirectory(NULL);

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

}

MHCalibrationPixel::~MHCalibrationPixel()
{

  delete fHChargeHiGain;
  delete fHTimeHiGain;
  delete fHChargevsNHiGain;

  delete fHChargeLoGain;
  delete fHTimeLoGain;
  delete fHChargevsNLoGain;

  delete fHiGains;
  delete fLoGains;

  if (fChargeGausFit)
    delete fChargeGausFit;
  if (fTimeGausFit)
    delete fTimeGausFit;
  if (fFitLegend)
    delete fFitLegend;
  if (fHivsLoGain)
    delete fHivsLoGain;


}


void MHCalibrationPixel::ChangeHistId(Int_t id)
{

  fPixId = id;

  //
  // Names Hi gain Histograms
  //
  TString nameQHiGain = TString(fHChargeHiGain->GetName());
  nameQHiGain += id;
  fHChargeHiGain->SetName(nameQHiGain.Data());

  TString nameTHiGain = TString(fHTimeHiGain->GetName());
  nameTHiGain += id;
  fHTimeHiGain->SetName(nameTHiGain.Data());

  TString nameQvsNHiGain  = TString(fHChargevsNHiGain->GetName());
  nameQvsNHiGain += id;
  fHChargevsNHiGain->SetName(nameQvsNHiGain.Data());

  //
  // Title Hi gain Histograms
  //
  TString titleQHiGain = TString(fHChargeHiGain->GetTitle());
  titleQHiGain += id;
  fHChargeHiGain->SetTitle(titleQHiGain.Data());

  TString titleTHiGain = TString(fHTimeHiGain->GetTitle());
  titleTHiGain += id;
  fHTimeHiGain->SetTitle(titleTHiGain.Data());

  TString titleQvsNHiGain  = TString(fHChargevsNHiGain->GetTitle());
  titleQvsNHiGain += id;
  fHChargevsNHiGain->SetTitle(titleQvsNHiGain.Data());

  //
  // Names Low Gain Histograms
  //
  TString nameQLoGain = TString(fHChargeLoGain->GetName());
  nameQLoGain += id;
  fHChargeLoGain->SetName(nameQLoGain.Data());

  TString nameTLoGain = TString(fHTimeLoGain->GetName());
  nameTLoGain += id;
  fHTimeLoGain->SetName(nameTLoGain.Data());

  TString nameQvsNLoGain  = TString(fHChargevsNLoGain->GetName());
  nameQvsNLoGain += id;
  fHChargevsNLoGain->SetName(nameQvsNLoGain.Data());

  //
  // Titles Low Gain Histograms
  //
  TString titleQLoGain = TString(fHChargeLoGain->GetTitle());
  titleQLoGain += id;
  fHChargeLoGain->SetTitle(titleQLoGain.Data());

  TString titleTLoGain = TString(fHTimeLoGain->GetTitle());
  titleTLoGain += id;
  fHTimeLoGain->SetTitle(titleTLoGain.Data());

  TString titleQvsNLoGain  = TString(fHChargevsNLoGain->GetTitle());
  titleQvsNLoGain += id;
  fHChargevsNLoGain->SetTitle(titleQvsNLoGain.Data());
}


void MHCalibrationPixel::Reset()
{
  
  for (Int_t i = fHChargeHiGain->FindBin(fChargeFirstHiGain); 
       i <= fHChargeHiGain->FindBin(fChargeLastHiGain); i++)
    fHChargeHiGain->SetBinContent(i, 1.e-20);

  for (Int_t i = 0; i < 16; i++)
      fHTimeHiGain->SetBinContent(i, 1.e-20);
  
  fChargeLastHiGain     = 9999.5;

  fHChargeHiGain->GetXaxis()->SetRangeUser(0.,fChargeLastHiGain);

 for (Int_t i = fHChargeLoGain->FindBin(fChargeFirstLoGain); 
       i <= fHChargeLoGain->FindBin(fChargeLastLoGain); i++)
    fHChargeLoGain->SetBinContent(i, 1.e-20);

  for (Int_t i = 0; i < 16; i++)
      fHTimeLoGain->SetBinContent(i, 1.e-20);
  
  fChargeLastLoGain     = 9999.5;

  fHChargeLoGain->GetXaxis()->SetRangeUser(0.,fChargeLastLoGain);

  return;
}


// -------------------------------------------------------------------------
//
// Set the binnings and prepare the filling of the histograms
//
Bool_t MHCalibrationPixel::SetupFill(const MParList *plist)
{

  fHChargeHiGain->Reset();
  fHTimeHiGain->Reset();
  fHChargeLoGain->Reset();
  fHTimeLoGain->Reset();

  return kTRUE;
}


Bool_t MHCalibrationPixel::UseLoGain()
{

  if (fHChargeHiGain->Integral() > fHChargeLoGain->Integral())
    {
      fUseLoGain = kFALSE;
      return kFALSE;
    }
  else
    {
      fUseLoGain = kTRUE;
      return kTRUE;
    }
}

void MHCalibrationPixel::SetPointInGraph(Float_t qhi,Float_t qlo)
{

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

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

}


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

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

  if (fFitOK) 
      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);

  fFitLegend->AddText(line1);


  const TString line4 =
    Form("Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fChargeSigma,fChargeSigmaErr);

  fFitLegend->AddText(line4);


  const TString line7 =
    Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChargeChisquare,fChargeNdf);

  fFitLegend->AddText(line7);


  const TString line8 =
    Form("Probability: %4.3f ",fChargeProb);

  fFitLegend->AddText(line8);

  if (fFitOK)
    fFitLegend->AddText("Result of the Fit: OK");
  else
    fFitLegend->AddText("Result of the Fit: NOT OK");

  fFitLegend->SetBit(kCanDelete);
  fFitLegend->Draw();

}


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

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

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

  fHChargeHiGain->DrawCopy(opt);
  
  c->Modified();
  c->Update();
  
  if (fUseLoGain)
    {

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

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

      fHChargeLoGain->DrawCopy(opt);
      
      if (fChargeGausFit)
        {
          if (fFitOK)
            fChargeGausFit->SetLineColor(kGreen);          
          else
            fChargeGausFit->SetLineColor(kRed);
          
          fChargeGausFit->DrawCopy("same");
        }
      c->Modified();
      c->Update();
      
      c->cd(3);
      gROOT->SetSelectedPad(NULL);
      gStyle->SetOptFit();
      fHivsLoGain->Draw("prof");

      gPad->Modified();
      gPad->Update();
      
      c->cd(4);
      DrawLegend();
      
    }
  else
    {
      if (fChargeGausFit)
        {
          if (fFitOK)
            fChargeGausFit->SetLineColor(kGreen);          
          else
            fChargeGausFit->SetLineColor(kRed);
          
          fChargeGausFit->DrawCopy("same");
        }

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

      if (fHChargeLoGain->Integral() > 0)
        gPad->SetLogy(1);
      else
        gPad->SetLogy(0);
      
      fHChargeLoGain->DrawCopy(opt);
      c->Modified();
      c->Update();
      
      c->cd(3);
      DrawLegend();
      
      c->cd(4);
      
      gROOT->SetSelectedPad(NULL);
      gStyle->SetOptFit();
      fHivsLoGain->Draw("prof");
      gPad->Modified();
      gPad->Update();
    }	
  
  c->Modified();
  c->Update();
  
  c->cd(5);
  gStyle->SetOptStat(1111111);

  gPad->SetTicks();  
  gPad->SetLogy(0);
  fHTimeHiGain->DrawCopy(opt);
  c->Modified();
  c->Update();
  
  if (fUseLoGain)
      {
        
	c->cd(6);
        gPad->SetTicks();
	gPad->SetLogy(0);
	fHTimeLoGain->DrawCopy(opt);
	c->Modified();
	c->Update();
	
	if (fTimeGausFit)
	  {
	    if (fTimeChisquare > 20.) 
	      fTimeGausFit->SetLineColor(kRed);
	    else
	      fTimeGausFit->SetLineColor(kGreen);
	    
	    fTimeGausFit->DrawCopy("same");
	    c->Modified();
	    c->Update();
	  }
      }
    else
      {  
	if (fTimeGausFit)
	  {
	    if (fTimeChisquare > 20.) 
	      fTimeGausFit->SetLineColor(kRed);
	    else
	      fTimeGausFit->SetLineColor(kGreen);
	    
	    fTimeGausFit->DrawCopy("same");
	    c->Modified();
	    c->Update();
	  }

	c->cd(6);
        gPad->SetTicks();
	gPad->SetLogy(0);
	fHTimeLoGain->DrawCopy(opt);	
	c->Modified();
	c->Update();
    
      }
    c->Modified();
    c->Update();
    
    c->cd(7);
    gPad->SetTicks();
    fHChargevsNHiGain->DrawCopy(opt);
    c->Modified();
    c->Update();

    c->cd(8);
    gPad->SetTicks();
    fHChargevsNLoGain->DrawCopy(opt);
    c->Modified();
    c->Update();


}

void MHCalibrationPixel::FitHiGainvsLoGain()
{

  if (fHivsLoGain)
    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.*fHivsLoGain->GetRMS());
  
  fOffset = fHivsLoGain->GetFunction("pol1")->GetParameter(0);
  fSlope  = fHivsLoGain->GetFunction("pol1")->GetParameter(1);
  
}


Bool_t MHCalibrationPixel::FitTimeHiGain(Axis_t rmin, Axis_t rmax, Option_t *option)
{

  if (fTimeGausFit)
    return kFALSE;

  rmin = (rmin != 0.) ? rmin : fTimeLowerFitRange;
  rmax = (rmax != 0.) ? rmax : fTimeUpperFitRange;

  const Stat_t entries       = fHTimeHiGain->Integral();
  const Double_t mu_guess    = fHTimeHiGain->GetBinCenter(fHTimeHiGain->GetMaximumBin());
  const Double_t sigma_guess = (rmax - rmin)/2.;
  const Double_t area_guess  = entries/gkSq2Pi;

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

  if (!fTimeGausFit) 
    {
    *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
    return kFALSE;
    }

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

  fHTimeHiGain->Fit(fTimeGausFit,option);

  rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2);
  rmax = fTimeGausFit->GetParameter(1) + 3.*fTimeGausFit->GetParameter(2);
  fTimeGausFit->SetRange(rmin,rmax);  

  fHTimeHiGain->Fit(fTimeGausFit,option);

  fTimeChisquare = fTimeGausFit->GetChisquare();
  fTimeNdf       = fTimeGausFit->GetNDF();
  fTimeProb      = fTimeGausFit->GetProb();

  fTimeMean      = fTimeGausFit->GetParameter(1);
  fTimeSigma     = fTimeGausFit->GetParameter(2);

  if (fTimeChisquare > 20.)  // Cannot use Probability because Ndf is sometimes < 1
    {
      *fLog << warn << "Fit of the Arrival times failed ! " << endl;
      return kFALSE;
    }
  
  return kTRUE;

}

Bool_t MHCalibrationPixel::FitTimeLoGain(Axis_t rmin, Axis_t rmax, Option_t *option)
{

  if (fTimeGausFit)
    return kFALSE;

  rmin = (rmin != 0.) ? rmin : fTimeLowerFitRange;
  rmax = (rmax != 0.) ? rmax : fTimeUpperFitRange;

  const Stat_t   entries     = fHTimeLoGain->Integral();
  const Double_t mu_guess    = fHTimeLoGain->GetBinCenter(fHTimeLoGain->GetMaximumBin());
  const Double_t sigma_guess = (rmax - rmin)/2.;
  const Double_t area_guess  = entries/gkSq2Pi;

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

  if (!fTimeGausFit) 
    {
    *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
    return kFALSE;
    }

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

  fHTimeLoGain->Fit(fTimeGausFit,option);

  rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2);
  rmax = fTimeGausFit->GetParameter(1) + 3.*fTimeGausFit->GetParameter(2);
  fTimeGausFit->SetRange(rmin,rmax);  

  fHTimeLoGain->Fit(fTimeGausFit,option);

  fTimeChisquare = fTimeGausFit->GetChisquare();
  fTimeNdf       = fTimeGausFit->GetNDF();
  fTimeProb      = fTimeGausFit->GetProb();

  fTimeMean      = fTimeGausFit->GetParameter(1);
  fTimeSigma     = fTimeGausFit->GetParameter(2);

  if (fTimeChisquare > 20.)  // Cannot use Probability because Ndf is sometimes < 1
    {
      *fLog << warn << "Fit of the Arrival times failed ! " << endl;
      return kFALSE;
    }
  
  return kTRUE;

}

Bool_t MHCalibrationPixel::FitChargeHiGain(Option_t *option)
{

  if (fChargeGausFit)
    return kFALSE;

  //
  // Get the fitting ranges
  //
  Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirstHiGain;
  Axis_t rmax = 0.;

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

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

  fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLastHiGain);

  if (!fChargeGausFit) 
    {
    *fLog << err << dbginf << "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,fChargeLastHiGain);
  fChargeGausFit->SetParLimits(2,0.,fChargeLastHiGain-rmin);
  fChargeGausFit->SetRange(rmin,fChargeLastHiGain);

  fHChargeHiGain->Fit(fChargeGausFit,option);
  
  Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.0*fChargeGausFit->GetParameter(2);
  
  rmin = (rtry < rmin ? rmin : rtry);
  rmax = fChargeGausFit->GetParameter(1) + 3.5*fChargeGausFit->GetParameter(2);
  fChargeGausFit->SetRange(rmin,rmax);  

  fHChargeHiGain->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
  // The Probability is greater than gkProbLimit (default 0.001 == 99.9%)
  //
  if (fChargeProb < gkProbLimit) 
    {
      *fLog << warn << "Prob: " << fChargeProb 
            << " is smaller than the allowed value: " << gkProbLimit << endl;
      fFitOK = kFALSE;
      return kFALSE;
    }
  
  fFitOK = kTRUE;
  return kTRUE;
}


Bool_t MHCalibrationPixel::FitChargeLoGain(Option_t *option)
{

  if (fChargeGausFit)
    return kFALSE;

  //
  // Get the fitting ranges
  //
  Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirstLoGain;
  Axis_t rmax = 0.;

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

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

  fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLastLoGain);

  if (!fChargeGausFit) 
    {
    *fLog << err << dbginf << "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,fChargeLastLoGain);
  fChargeGausFit->SetParLimits(2,0.,fChargeLastLoGain-rmin);
  fChargeGausFit->SetRange(rmin,fChargeLastLoGain);

  fHChargeLoGain->Fit(fChargeGausFit,option);
  
  Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.*fChargeGausFit->GetParameter(2);
  
  rmin = (rtry < rmin ? rmin : rtry);
  rmax = fChargeGausFit->GetParameter(1) + 3.5*fChargeGausFit->GetParameter(2);
  fChargeGausFit->SetRange(rmin,rmax);  

  fHChargeLoGain->Fit(fChargeGausFit,option);

  //  rmin = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2);
  //  rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2);
  //  fChargeGausFit->SetRange(rmin,rmax);  

  // fHChargeLoGain->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
  // The Probability is greater than gkProbLimit (default 0.01 == 99%)
  //
  if (fChargeProb < gkProbLimit) 
    {
      *fLog << warn << "Prob: " << fChargeProb << " is smaller than the allowed value: " << gkProbLimit << endl;
      fFitOK = kFALSE;
      return kFALSE;
    }
  
  
  fFitOK = kTRUE;
    
  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);
  fChargeNbinsHiGain = nbins;

  CutEdges(fHChargeLoGain,nbins);

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

  CutEdges(fHChargevsNHiGain,0);
  CutEdges(fHChargevsNLoGain,0);

}

void MHCalibrationPixel::PrintChargeFitResult()
{
  
  *fLog << "Results of the Summed Charges Fit: "                      << endl;
  *fLog << "Chisquare: "        << fChargeChisquare                   << endl;
  *fLog << "DoF: "              << fChargeNdf                         << endl;
  *fLog << "Probability: "      << fChargeProb                        << endl;
  *fLog                                                               << endl;
  
}

void MHCalibrationPixel::PrintTimeFitResult()
{

  *fLog << "Results of the Time Slices Fit: "                        << endl;
  *fLog << "Chisquare: "   << fTimeChisquare                         << endl;
  *fLog << "Ndf: "         << fTimeNdf                               << endl;
  *fLog << "Probability: " << fTimeProb                              << endl;
  *fLog                                                              << endl;

}