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

//////////////////////////////////////////////////////////////////////////////
//
//  MHCalibrationBlindPixel
//
//  Performs all the Single Photo-Electron Fit to extract
//  the mean number of photons and to derive the light flux
//
// The fit result is accepted under condition that:
// 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
// 2) at least 100 events are in the single Photo-electron peak
//
// Used numbers are the following:
//
// Electronic conversion factor:
//   Assume, we have N_e electrons at the anode, 
//   thus a charge of N_e*e (e = electron charge) Coulomb.
//
//   This charge is AC coupled and runs into a R_pre = 50 Ohm resistency. 
//   The corresponding current is amplified by a gain factor G_pre = 400 
//   (the precision of this value still has to be checked !!!) and again AC coupled to 
//   the output. 
//   The corresponding signal goes through the whole transmission and 
//   amplification chain and is digitized in the FADCs. 
//   The conversion Signal Area to FADC counts (Conv_trans) has been measured 
//   by David and Oscar to be approx. 3.9 pVs^-1
//
//   Thus: Conversion FADC counts to Number of Electrons at Anode: 
//         FADC counts = (1/Conv_tran) * G_pre * R_pre *  e * N_e = 8 * 10^-4 N_e. 
//
//   Also: FADC counts = 8*10^-4 * GAIN * N_phe
//
//   In the blind pixel, there is an additional pre-amplifier with an amplification of 
//   about 10. Therefore, we have for the blind pixel:
//
//   FADC counts (Blind Pixel) = 8*10^-3 * GAIN * N_phe
//
//////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationBlindPixel.h"
#include "MHCalibrationConfig.h"


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

#include <TGraph.h>
#include <TF1.h>
#include <TH1.h>
#include <TRandom.h>

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

ClassImp(MHCalibrationBlindPixel);

using namespace std;

const Int_t    MHCalibrationBlindPixel::fgBlindPixelChargeNbins        = 1000;
const Int_t    MHCalibrationBlindPixel::fgBlindPixelTimeNbins          = 22;
const Axis_t   MHCalibrationBlindPixel::fgBlindPixelTimeFirst          = -9.00;
const Axis_t   MHCalibrationBlindPixel::fgBlindPixelTimeLast           = 12.00;
const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmp      = 0.008;
const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002;
  
const Int_t   MHCalibrationBlindPixel::fPSDNbins          = 50;
const Int_t   MHCalibrationBlindPixel::fPulserFrequency   = 500;  
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title)
    :   fHBlindPixelPSD(NULL), 
        fSinglePheFit(NULL),  
        fTimeGausFit(NULL),
        fSinglePhePedFit(NULL),
        fPSDHiGain(NULL),
        fPSDLoGain(NULL),
        fHPSD(NULL),
        fPSDExpFit(NULL),
        fChargeXaxis(NULL),
        fPSDXaxis(NULL),
        fCurrentSize(1024),
        fFitLegend(NULL)
{

    fName  = name  ? name  : "MHCalibrationBlindPixel";
    fTitle = title ? title : "Fill the accumulated charges and times all Blind Pixel events and perform fits";

    // Create a large number of bins, later we will rebin
    fBlindPixelChargefirst = -400.;
    fBlindPixelChargelast  = 1600.;

    fHBlindPixelCharge = new TH1I("HBlindPixelCharge","Distribution of Summed FADC Slices",
				  fgBlindPixelChargeNbins,fBlindPixelChargefirst,fBlindPixelChargelast);
    fHBlindPixelCharge->SetXTitle("Sum FADC Slices");
    fHBlindPixelCharge->SetYTitle("Nr. of events");
    fHBlindPixelCharge->Sumw2();
    fHBlindPixelCharge->SetDirectory(NULL);

    fHBlindPixelTime = new TH1F("HBlindPixelTime","Distribution of Mean Arrival Times",
                                fgBlindPixelTimeNbins,fgBlindPixelTimeFirst,fgBlindPixelTimeLast);
    fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
    fHBlindPixelTime->SetYTitle("Nr. of events");
    fHBlindPixelTime->Sumw2();
    fHBlindPixelTime->SetDirectory(NULL);

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

    fHSinglePheFADCSlices = new TH1I("HSinglePheFADCSlices","FADC slices Single Phe events",30,0.5,30.5);  
    fHSinglePheFADCSlices->SetXTitle("FADC slice");
    fHSinglePheFADCSlices->SetYTitle("Counts");    
    fHSinglePheFADCSlices->SetDirectory(NULL);    

    fHPedestalFADCSlices  = new TH1I("HPedestalFADCSlices", "FADC slices Pedestal events",30,0.5,30.5);  
    fHPedestalFADCSlices->SetXTitle("FADC slice");
    fHPedestalFADCSlices->SetYTitle("Counts");    
    fHPedestalFADCSlices->SetDirectory(NULL);    

    Clear();
}

MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
{

  if (fFitLegend)
    delete fFitLegend;

  delete fHBlindPixelCharge;
  delete fHBlindPixelTime;

  delete fHiGains;
  delete fLoGains;

  delete fHSinglePheFADCSlices;
  delete fHPedestalFADCSlices;

  if (fHBlindPixelPSD)
    delete fHBlindPixelPSD;
  if (fSinglePheFit)  
    delete fSinglePheFit;
  if (fTimeGausFit)
    delete fTimeGausFit;
  if(fSinglePhePedFit)
    delete fSinglePhePedFit;
  if (fPSDExpFit)
    delete fPSDExpFit;
  if (fHPSD)
    delete fHPSD;
  if (fChargeXaxis)
    delete fChargeXaxis;
  if (fPSDXaxis)
    delete fPSDXaxis;

}

void MHCalibrationBlindPixel::Clear(Option_t *o)
{

  fTotalEntries            = 0;
  fCurrentSize             = 1024;
  
  fBlindPixelChargefirst = -400.;
  fBlindPixelChargelast  = 1600.;

  fLambda    = 0.;
  fMu0       = 0.;
  fMu1       = 0.;
  fSigma0    = 0.;
  fSigma1    = 0.;
  
  fLambdaErr = 0.;
  fMu0Err    = 0.;
  fMu1Err    = 0.;
  fSigma0Err = 0.;
  fSigma1Err = 0.;

  fChisquare = -1.;
  fProb      = -1.;
  fNdf       = -1;

  fMeanTime     = -1.;
  fMeanTimeErr  = -1.;
  fSigmaTime    = -1.;
  fSigmaTimeErr = -1.;
  
  fLambdaCheck    = -1.;
  fLambdaCheckErr = -1.;
  
  fMeanPedestal     = 0.;
  fMeanPedestalErr  = 0.;
  fSigmaPedestal    = 0.;
  fSigmaPedestalErr = 0.;

  fFitFunc = kEPoisson5;

  if (fFitLegend)
    delete fFitLegend;
  if (fHBlindPixelPSD)
    delete fHBlindPixelPSD;
  if (fSinglePheFit)  
    delete fSinglePheFit;
  if (fTimeGausFit)
    delete fTimeGausFit;
  if(fSinglePhePedFit)
    delete fSinglePhePedFit;
  if (fPSDExpFit)
    delete fPSDExpFit;
  if (fHPSD)
    delete fHPSD;
  if (fChargeXaxis)
    delete fChargeXaxis;
  if (fPSDXaxis)
    delete fPSDXaxis;
  if (fPSDHiGain)
    delete fPSDHiGain;
  if (fPSDLoGain)
    delete fPSDLoGain;

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

  return;
}

void MHCalibrationBlindPixel::Reset()
{
  
  Clear();

  fHBlindPixelCharge->Reset();
  fHBlindPixelTime->Reset();
  
  fHiGains->Set(1024);
  fLoGains->Set(1024);

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


}

Bool_t MHCalibrationBlindPixel::CheckOscillations()
{

  if (fPSDExpFit)
    return IsOscillating();

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

  fHPSD = ProjectArray(*fPSDHiGain, fPSDNbins,
                       "PSDProjectionBlindPixel",
                       "Power Spectrum Density Projection HiGain Blind Pixel");

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

  fPSDExpFit = new TF1("PSDExpFit","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 MHCalibrationBlindPixel::CreatePSDXaxis(Int_t n)
{
  
  if (fPSDXaxis)
    return;

  fPSDXaxis = new TArrayF(n);

  for (Int_t i=0;i<n;i++)
  for (Int_t i=0;i<n;i++)
    fPSDXaxis->AddAt((Float_t)(fPulserFrequency*i)/2./n,i);
}

void MHCalibrationBlindPixel::CreateChargeXaxis(Int_t n)
{
  
  if (!fChargeXaxis)
    {
      fChargeXaxis = new TArrayF(n);
      for (Int_t i=0;i<n;i++)
        fChargeXaxis->AddAt((Float_t)i/fPulserFrequency,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;i<n+diff;i++)
      fChargeXaxis->AddAt((Float_t)i/fPulserFrequency,i);
}

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

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

void MHCalibrationBlindPixel::CutArrayBorder(TArrayI *array) const 
{
  
  Int_t i;

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

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

const Bool_t MHCalibrationBlindPixel::IsOscillating()
{

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

  return CheckOscillations();

}

Bool_t MHCalibrationBlindPixel::FillGraphs(const Int_t qhi, const Int_t qlo)
{

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

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

  fTotalEntries++;

  return kTRUE;

}



Bool_t MHCalibrationBlindPixel::FillBlindPixelCharge(const Int_t q)
{
    return fHBlindPixelCharge->Fill(q) > -1;
}

Bool_t MHCalibrationBlindPixel::FillBlindPixelTime(const Float_t t)
{
    return fHBlindPixelTime->Fill(t) > -1;
}


// -------------------------------------------------------------------------
//
// Draw a legend with the fit results
//
void MHCalibrationBlindPixel::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 single PhE Fit (to k=6):");
  fFitLegend->SetTextSize(0.05);

  const TString line1 = 
  Form("Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr());
  TText *t1 = fFitLegend->AddText(line1.Data());
  t1->SetBit(kCanDelete);

  const TString line6 =
  Form("Mean #lambda (check) = %2.2f #pm %2.2f",GetLambdaCheck(),GetLambdaCheckErr());
  TText *t2 = fFitLegend->AddText(line6.Data());
  t2->SetBit(kCanDelete);

  const TString line2 = 
  Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err());
  TText *t3 = fFitLegend->AddText(line2.Data());
  t3->SetBit(kCanDelete);

  const TString line3 =
  Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err());
  TText *t4 = fFitLegend->AddText(line3.Data());
  t4->SetBit(kCanDelete);

  const TString line4 =
  Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err());
  TText *t5 = fFitLegend->AddText(line4.Data());
  t5->SetBit(kCanDelete);

  const TString line5 =
  Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err());
  TText *t6 = fFitLegend->AddText(line5.Data());
  t6->SetBit(kCanDelete);

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

  const TString line8 =
  Form("Probability: %4.2f ",GetProb());
  TText *t8 = fFitLegend->AddText(line8.Data());
  t8->SetBit(kCanDelete);

  if (IsFitOK())
    {
      TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
      t->SetBit(kCanDelete);
    }
   else
     {
       TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
       t->SetBit(kCanDelete);
     }
  
  fFitLegend->SetBit(kCanDelete);
  fFitLegend->Draw();
  
  return;
}

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

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

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

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

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

    gStyle->SetOptFit(1);
    gStyle->SetOptStat(111111);

    TCanvas *c = MakeDefCanvas(this,550,700);

    c->Divide(2,5);

    gROOT->SetSelectedPad(NULL);

    c->cd(1);
    gPad->SetLogx(0);
    gPad->SetLogy(1);
    gPad->SetTicks();

    fHBlindPixelCharge->Draw(opt);
    
    if (IsFitOK())
      fSinglePheFit->SetLineColor(kGreen);          
    else
      fSinglePheFit->SetLineColor(kRed);
    
    fSinglePheFit->Draw("same");
    c->Modified();
    c->Update();

    fSinglePhePedFit->SetLineColor(kBlue);
    fSinglePhePedFit->Draw("same");

    c->cd(2);
    DrawLegend();
    

    c->cd(3);
    gPad->SetLogy(0);
    gPad->SetBorderMode(0);
    //    fHBlindPixelTime->Draw(opt);
    fHSinglePheFADCSlices->Draw(opt);

    c->cd(4);
    gPad->SetLogy(0);
    gPad->SetBorderMode(0);
    //    fHBlindPixelTime->Draw(opt);
    fHPedestalFADCSlices->Draw(opt);
    
    c->Modified();
    c->Update();

    CutArrayBorder(fHiGains);
    CreateChargeXaxis(fHiGains->GetSize());
    
    c->cd(5);
    gPad->SetTicks();
    TGraph *gr1 = new TGraph(fChargeXaxis->GetSize(),
                             fChargeXaxis->GetArray(),
                             fHiGains->GetArray());  
    gr1->SetTitle("Evolution of HiGain charges with time");
    gr1->GetXaxis()->SetTitle("Time [s]");
    gr1->GetYaxis()->SetTitle("Sum FADC slices");  
    gr1->SetBit(kCanDelete);
    gr1->Draw("AL");
    
    CutArrayBorder(fLoGains);  
    CreateChargeXaxis(fHiGains->GetSize());
    
    c->cd(6);
    gPad->SetTicks();
    TGraph *gr2 = new TGraph(fChargeXaxis->GetSize(),
                             fChargeXaxis->GetArray(),
                             fLoGains->GetArray());  
    gr2->SetTitle("Evolution of LoGain charges with time");
    gr2->GetXaxis()->SetTitle("Time [s]");
    gr2->GetYaxis()->SetTitle("Sum FADC slices");  
    gr2->SetBit(kCanDelete);
    gr2->Draw("AL");
  
    c->Modified();
    c->Update();
    
    c->cd(7);
  
    TArrayF *array;
    
    if (!fPSDHiGain)
      return;
    array = fPSDHiGain;
  
    if (!fPSDXaxis)
      CreatePSDXaxis(array->GetSize());

    TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray());
    gr3->SetTitle("Power Spectrum Density");
    gr3->GetXaxis()->SetTitle("Freq. [Hz]");
    gr3->GetYaxis()->SetTitle("P(f)");  
    gr3->SetBit(kCanDelete);
    gr3->Draw("AL");
    
    c->Modified();
    c->Update();
    
    c->cd(8);
    
    gStyle->SetOptStat(111111);
    gPad->SetTicks();  
    
    if (fHPSD->Integral() > 0)
      gPad->SetLogy();
    
    fHPSD->Draw(opt);
    fHPSD->GetXaxis()->SetTitle("P(f)");
    fHPSD->GetYaxis()->SetTitle("Counts");  
    
    if (fPSDExpFit)
      {
        fPSDExpFit->SetLineColor(IsOscillating() ? kRed : kGreen);          
        fPSDExpFit->Draw("same");
      }
    
    c->cd(9);
    
    array = fPSDLoGain;
  
    if (!fPSDXaxis)
      CreatePSDXaxis(array->GetSize());

    TGraph *gr4 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray());
    gr4->SetTitle("Power Spectrum Density");
    gr4->GetXaxis()->SetTitle("Freq. [Hz]");
    gr4->GetYaxis()->SetTitle("P(f)");  
    gr4->SetBit(kCanDelete);
    gr4->Draw("AL");
    
    c->Modified();
    c->Update();

    c->cd(10);
}



Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1) 
{
  gRandom->SetSeed();

  if (fHBlindPixelCharge->GetIntegral() != 0)
    {
      *fLog << err << "Histogram " << fHBlindPixelCharge->GetTitle() << " is already filled. " << endl;
      *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
      return kFALSE;
    }

  if (!InitFit(fBlindPixelChargefirst,fBlindPixelChargelast))
    return kFALSE;

  for (Int_t i=0;i<10000; i++) 
    fHBlindPixelCharge->Fill(fSinglePheFit->GetRandom());
  
  return kTRUE;
}

Bool_t MHCalibrationBlindPixel::InitFit(Axis_t min, Axis_t max)
{
  
  *fLog << inf << "min: " << min << endl;
  *fLog << "max: " << max << endl;

  //
  // 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      = fHBlindPixelCharge->Integral("width");
  const Double_t lambda_guess = 0.1;
  const Double_t maximum_bin  = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
  const Double_t norm         = entries/gkSq2Pi;

  //
  // Initialize the fit function
  //
  switch (fFitFunc)
    {
    case kEPoisson4:
      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6);
      break;
    case kEPoisson5:
      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6);
      break;
    case kEPoisson6:
      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6);
      break;
    case kEPolya:
      fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8);
      break;
    case kEMichele:
      break;

    default:
      *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
      return kFALSE;
      break;
    }
  
  const Double_t mu_0_guess = maximum_bin;
  const Double_t si_0_guess = 40.;
  const Double_t mu_1_guess = mu_0_guess + 100.;
  const Double_t si_1_guess = si_0_guess + si_0_guess;
  // Michele
  const Double_t lambda_1cat_guess = 0.5;
  const Double_t lambda_1dyn_guess = 0.5;
  const Double_t mu_1cat_guess = mu_0_guess + 50.;
  const Double_t mu_1dyn_guess = mu_0_guess + 20.;
  const Double_t si_1cat_guess = si_0_guess + si_0_guess;
  const Double_t si_1dyn_guess = si_0_guess;
  // Polya
  const Double_t excessPoisson_guess = 0.5;
  const Double_t delta1_guess     = 8.;
  const Double_t delta2_guess     = 5.;
  const Double_t electronicAmp_guess  = fgBlindPixelElectronicAmp;
  const Double_t electronicAmp_limit  = fgBlindPixelElectronicAmpError;

  *fLog << inf << "pedestal: " << fMeanPedestal << endl;
  *fLog << "sigma: "    << fSigmaPedestal << endl;

  //
  // Initialize boundaries and start parameters
  //
  switch (fFitFunc)
    {
      
    case kEPoisson4:
      if ((fMeanPedestal) && (fSigmaPedestal))
        fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm);
      else
        fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);

      fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");

      fSinglePheFit->SetParLimits(0,0.,0.5);

      if ((fMeanPedestal) && (fSigmaPedestal))
        fSinglePheFit->SetParLimits(1,
                                    fMeanPedestal-5.*fMeanPedestalErr,
                                    fMeanPedestal+5.*fMeanPedestalErr);
      else
        fSinglePheFit->SetParLimits(1,-3.,0.);

      fSinglePheFit->SetParLimits(2,min,max);

      if ((fMeanPedestal) && (fSigmaPedestal))
        fSinglePheFit->SetParLimits(3,
                                    fSigmaPedestal-5.*fSigmaPedestalErr,
                                    fSigmaPedestal+5.*fSigmaPedestalErr);
      else
        fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);

      fSinglePheFit->SetParLimits(4,0.,(max-min));
      fSinglePheFit->SetParLimits(5,norm-(0.5*norm),norm+(0.5*norm));
      break;
    case kEPoisson5:
    case kEPoisson6:
      fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
      fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
      fSinglePheFit->SetParLimits(0,0.,1.);
      fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
      fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
      fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
      fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
      fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
      break;

    case kEPolya:
      if ((fMeanPedestal) && (fSigmaPedestal))
        fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
                                     delta1_guess,delta2_guess,
                                     electronicAmp_guess,
                                     fSigmaPedestal,
                                     norm, 
                                     fMeanPedestal);
      else
        fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
                                     delta1_guess,delta2_guess,
                                     electronicAmp_guess,
                                     si_0_guess,
                                     norm, mu_0_guess);
      fSinglePheFit->SetParNames("#lambda","b_{tot}",
                                 "#delta_{1}","#delta_{2}",
                                 "amp_{e}","#sigma_{0}",
                                 "Area", "#mu_{0}");
      fSinglePheFit->SetParLimits(0,0.,1.);
      fSinglePheFit->SetParLimits(1,0.,1.); 
      fSinglePheFit->SetParLimits(2,6.,12.);    
      fSinglePheFit->SetParLimits(3,3.,8.);    
      fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
                                    electronicAmp_guess+electronicAmp_limit);    
      if ((fMeanPedestal) && (fSigmaPedestal))
        fSinglePheFit->SetParLimits(5,
                                    fSigmaPedestal-3.*fSigmaPedestalErr,
                                    fSigmaPedestal+3.*fSigmaPedestalErr);
      else
        fSinglePheFit->SetParLimits(5,min,(max-min)/1.5);

      fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
      if ((fMeanPedestal) && (fSigmaPedestal))
        fSinglePheFit->SetParLimits(7,
                                    fMeanPedestal-3.*fMeanPedestalErr,
                                    fMeanPedestal+3.*fMeanPedestalErr);
      else
        fSinglePheFit->SetParLimits(7,-35.,15.);
      break;

    case kEMichele:


      break;

    default:
      *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
      return kFALSE;
      break;
    }

  return kTRUE;
}

void MHCalibrationBlindPixel::ExitFit(TF1 *f)
{
  

  //
  // Finalize
  //
  switch (fFitFunc)
    {
      
    case kEPoisson4:
    case kEPoisson5:
    case kEPoisson6:
    case kEPoisson7:
      fLambda = fSinglePheFit->GetParameter(0);
      fMu0    = fSinglePheFit->GetParameter(1);
      fMu1    = fSinglePheFit->GetParameter(2);
      fSigma0 = fSinglePheFit->GetParameter(3);
      fSigma1 = fSinglePheFit->GetParameter(4);
      
      fLambdaErr = fSinglePheFit->GetParError(0);
      fMu0Err    = fSinglePheFit->GetParError(1);
      fMu1Err    = fSinglePheFit->GetParError(2);
      fSigma0Err = fSinglePheFit->GetParError(3);
      fSigma1Err = fSinglePheFit->GetParError(4);
      break;
    case kEPolya:
      fLambda =  fSinglePheFit->GetParameter(0);
      fMu0    =  fSinglePheFit->GetParameter(7);
      fMu1    = 0.;
      fSigma0 =  fSinglePheFit->GetParameter(5);
      fSigma1 = 0.;

      fLambdaErr = fSinglePheFit->GetParError(0);
      fMu0Err    = fSinglePheFit->GetParError(7);
      fMu1Err    = 0.;
      fSigma0Err = fSinglePheFit->GetParError(5);
      fSigma1Err = 0.;
    default:
      break;
    }
  
}


Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt) 
{

  //
  // Get the fitting ranges
  //
  rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst;
  rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast;
  if (!InitFit(rmin,rmax))
    return kFALSE;

  fHBlindPixelCharge->Fit(fSinglePheFit,opt);

  ExitFit(fSinglePheFit);

  fProb      = fSinglePheFit->GetProb();
  fChisquare = fSinglePheFit->GetChisquare();
  fNdf       = fSinglePheFit->GetNDF();

  // Perform the cross-check fitting only the pedestal:
  fSinglePhePedFit = new TF1("GausPed","gaus",rmin,0.);
  fHBlindPixelCharge->Fit(fSinglePhePedFit,opt);

  const Stat_t entries      = fHBlindPixelCharge->Integral("width");

  Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2);
  fLambdaCheck     = TMath::Log(entries/pedarea);
  fLambdaCheckErr  = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0)
                     + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2);

  *fLog << inf << "Results of the Blind Pixel Fit: " << endl;
  *fLog << "Chisquare: " << fChisquare << endl;
  *fLog << "DoF: " << fNdf << endl;
  *fLog << "Probability: " << fProb << endl;

  //
  // The fit result is accepted under condition that: 
  // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
  // 2) at least 50 events are in the single Photo-electron peak
  //
  if (fProb < gkProbLimit) 
    {
      *fLog << warn << "WARNING - Fit Probability " << fProb
            << " is smaller than the allowed value: " << gkProbLimit << endl;
      CLRBIT(fFlags,kFitOK);
      return kFALSE;
    }

  Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
  
  if (contSinglePhe < 50.) 
    {
      *fLog << warn << "WARNING - Statistics is too low: Only " << contSinglePhe
            << " in the Single Photo-Electron peak " << endl;
      CLRBIT(fFlags,kFitOK);      
      return kFALSE;
    } 
  else
    *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
  
  SETBIT(fFlags,kFitOK);
    
  return kTRUE;
}

 
void MHCalibrationBlindPixel::CutAllEdges()
{

  Int_t nbins = 25;

  StripZeros(fHBlindPixelCharge,nbins);

  fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
  fBlindPixelChargelast  = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);

}

Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt) 
{
  
  rmin = (rmin != 0.) ? rmin : 4.;
  rmax = (rmax != 0.) ? rmax : 9.;

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

  fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax);  
  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);

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

  fHBlindPixelTime->Fit(fTimeGausFit,opt);

  fMeanTime     = fTimeGausFit->GetParameter(2);
  fSigmaTime    = fTimeGausFit->GetParameter(3);
  fMeanTimeErr  = fTimeGausFit->GetParError(2);
  fSigmaTimeErr = fTimeGausFit->GetParError(3);

  *fLog << inf << "Results of the Times Fit: " << endl;
  *fLog << "Chisquare: "   << fTimeGausFit->GetChisquare() << endl;
  *fLog << "Ndf: "         << fTimeGausFit->GetNDF() << endl;

  return kTRUE;

}

