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

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationPixel::MHCalibrationPixel(const char *name, const char *title)
      : fPixId(-1),
        fChargeNbinsHiGain(2100),
        fChargeNbinsLoGain(1010),
        fTimeNbins(32),
        fChargevsNbins(1000),
        fTimeFirst(-0.25),
        fTimeLast(15.75),
	fHivsLoGain(NULL),
        fHPSD(NULL),
        fChargeGausFit(NULL),
	fTimeGausFit(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();

    fHTimeHiGain = new TH1F("HTimeHiGain","Distribution of Mean Arrival Hi Gain Times Pixel ",
			    fTimeNbins,fTimeFirst,fTimeLast);
    fHTimeLoGain = new TH1F("HTimeLoGain","Distribution of Mean Arrival Lo Gain Times Pixel ",
			    fTimeNbins,fTimeFirst,fTimeLast);

    fHTimeHiGain->SetXTitle("Mean Arrival Times [Hi Gain FADC slice nr]");
    fHTimeLoGain->SetXTitle("Mean Arrival Times [Lo Gain FADC slice nr]");

    fHTimeHiGain->SetYTitle("Nr. of events");
    fHTimeLoGain->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);
    fHTimeHiGain->SetDirectory(NULL);
    fHTimeLoGain->SetDirectory(NULL);
    fHChargevsNHiGain->SetDirectory(NULL);
    fHChargevsNLoGain->SetDirectory(NULL);

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

    Clear();

}

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::Clear(Option_t *o)
{
  
  fTotalEntries            = 0;

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

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

  fTimeChisquare           = -1.;
  fTimeProb                = -1.;
  fTimeNdf                 = -1;
  fTimeMean                = -1.;
  fTimeSigma               = -1.;

  fTimeLowerFitRangeHiGain = 0;
  fTimeUpperFitRangeHiGain = 0;
  fTimeLowerFitRangeLoGain = 0;
  fTimeUpperFitRangeLoGain = 0;

  fOffset = 0.;
  fSlope  = 0.;

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

  return;
}


void MHCalibrationPixel::Reset()
{
  
  Clear();
  
  fHChargeHiGain->Reset();
  fHChargeLoGain->Reset();
  fHTimeHiGain->Reset();
  fHTimeLoGain->Reset();
  fHChargevsNHiGain->Reset();
  fHChargevsNLoGain->Reset();

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

}



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

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

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

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

Bool_t MHCalibrationPixel::FillTimeLoGain(Float_t t)
{
    return (fHTimeLoGain->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::FillTimeHiGain(Float_t t)
{
    return (fHTimeHiGain->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)
{

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


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

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

  fHChargeHiGain->Draw(opt);
  
  c->Modified();
  c->Update();
  
  if (IsUseLoGain())
    {

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

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

      fHChargeLoGain->Draw(opt);
      
      if (fChargeGausFit)
        {
          if (IsFitOK())
            fChargeGausFit->SetLineColor(kGreen);          
          else
            fChargeGausFit->SetLineColor(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)
        {
          if (IsFitOK())
            fChargeGausFit->SetLineColor(kGreen);          
          else
            fChargeGausFit->SetLineColor(kRed);
          
          fChargeGausFit->Draw("same");
        }

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

      if (fHChargeLoGain->Integral() > 0)
        gPad->SetLogy(1);
      else
        gPad->SetLogy(0);
      
      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);
  fHTimeHiGain->Draw(opt);
  c->Modified();
  c->Update();
  
  if (IsUseLoGain())
      {
        
	c->cd(6);
        gPad->SetTicks();
	gPad->SetLogy(0);
	fHTimeLoGain->Draw(opt);
	c->Modified();
	c->Update();
	
	if (fTimeGausFit)
	  {
	    if (fTimeChisquare > 20.) 
	      fTimeGausFit->SetLineColor(kRed);
	    else
	      fTimeGausFit->SetLineColor(kGreen);
	    
	    fTimeGausFit->Draw("same");
	    c->Modified();
	    c->Update();
	  }
      }
    else
      {  
	if (fTimeGausFit)
	  {
	    if (fTimeChisquare > 20.) 
	      fTimeGausFit->SetLineColor(kRed);
	    else
	      fTimeGausFit->SetLineColor(kGreen);
	    
	    fTimeGausFit->Draw("same");
	    c->Modified();
	    c->Update();
	  }

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

    c->cd(8);
    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.0*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 : (Axis_t)fTimeLowerFitRangeHiGain;
  rmax = (rmax != 0.) ? rmax : (Axis_t)fTimeUpperFitRangeHiGain;

  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 << warn << dbginf << "WARNING: Could not create fit function for Time 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 (TMath::IsNaN(fTimeMean) || TMath::IsNaN(fTimeSigma))
    return kFALSE;

  if (TMath::IsNaN(fTimeChisquare) || fTimeChisquare > 20.)  // Cannot use Probability because Ndf is sometimes < 1
    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 : (Axis_t)fTimeLowerFitRangeLoGain;
  rmax = (rmax != 0.) ? rmax : (Axis_t)fTimeUpperFitRangeLoGain;

  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 << warn << dbginf << "WARNING: Could not create fit function for Time 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 << "WARNING: Fit of the Arrival times failed ! " << endl;
      return kFALSE;
    }
  
  return kTRUE;

}

Bool_t MHCalibrationPixel::FitCharge(Option_t *option)
{

  if (fChargeGausFit)
    return kFALSE;

  //
  // Get the fitting ranges
  //
  Axis_t rmin = fChargeFirstHiGain;
  if (TESTBIT(fFlags,kUseLoGain))
    rmin = fChargeFirstLoGain;

  Axis_t rmax = fChargeLastHiGain;
  if (TESTBIT(fFlags,kUseLoGain))
    rmin = fChargeFirstLoGain;

  TH1F *hist = fHChargeHiGain;
  if (TESTBIT(fFlags,kUseLoGain))
    hist = fHChargeLoGain;

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

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

  fChargeGausFit = new TF1(name.Data(),"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() < gkProbLimit)
    {

      Axis_t rtry = fChargeGausFit->GetParameter(1) - 3.0*fChargeGausFit->GetParameter(2);
      rmin        = (rtry < rmin ? rmin : rtry);
      rmax        = fChargeGausFit->GetParameter(1) + 3.0*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 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 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(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: "   << fTimeChisquare                         << endl;
  *fLog << all << "Ndf: "         << fTimeNdf                               << endl;
  *fLog << all << "Probability: " << fTimeProb                              << endl;
  *fLog << all                                                              << endl;

}