/* ======================================================================== *\
!
! *
! * 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
//
//////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationBlindPixel.h"
#include "MHCalibrationConfig.h"

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

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

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

ClassImp(MHCalibrationBlindPixel);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title)
    :   fBlindPixelChargeNbins(1000),
        fBlindPixelTimeNbins(32),
        fBlindPixelChargevsNbins(10000),
        fBlindPixelTimeFirst(-0.25),
        fBlindPixelTimeLast(15.75),
        fHBlindPixelPSD(NULL), 
        fSinglePheFit(NULL),  
        fTimeGausFit(NULL),
        fSinglePhePedFit(NULL),
        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 = -200.;
    fBlindPixelChargelast  = 800.;

    fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices",
				  fBlindPixelChargeNbins,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",
                                fBlindPixelTimeNbins,fBlindPixelTimeFirst,fBlindPixelTimeLast);
    fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
    fHBlindPixelTime->SetYTitle("Nr. of events");
    fHBlindPixelTime->Sumw2();
    fHBlindPixelTime->SetDirectory(NULL);

    fHBlindPixelChargevsN = new TH1I("HBlindPixelChargevsN","Sum of Charges vs. Event Number",
                                     fBlindPixelChargevsNbins,-0.5,(Axis_t)fBlindPixelChargevsNbins-0.5);
    fHBlindPixelChargevsN->SetXTitle("Event Nr.");
    fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices");
    fHBlindPixelChargevsN->SetDirectory(NULL);

    Clear();
}

MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
{

  if (fFitLegend)
    delete fFitLegend;

  delete fHBlindPixelCharge;
  delete fHBlindPixelTime;
  delete fHBlindPixelChargevsN;

  if (fHBlindPixelPSD)
    delete fHBlindPixelPSD;
  if (fSinglePheFit)  
    delete fSinglePheFit;
  if (fTimeGausFit)
    delete fTimeGausFit;
  if(fSinglePhePedFit)
    delete fSinglePhePedFit;

}

void MHCalibrationBlindPixel::Clear(Option_t *o)
{
  
  fBlindPixelChargefirst = -200.;
  fBlindPixelChargelast  = 800.;

  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 = kEPoisson4;

  if (fFitLegend)
    delete fFitLegend;
  if (fHBlindPixelPSD)
    delete fHBlindPixelPSD;
  if (fSinglePheFit)  
    delete fSinglePheFit;
  if (fTimeGausFit)
    delete fTimeGausFit;
  if(fSinglePhePedFit)
    delete fSinglePhePedFit;

  return;
}

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

  fHBlindPixelCharge->Reset();
  fHBlindPixelTime->Reset();
  fHBlindPixelChargevsN->Reset();
  
}

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

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

Bool_t MHCalibrationBlindPixel::FillBlindPixelChargevsN(Stat_t rq, Int_t t)
{
    return fHBlindPixelChargevsN->Fill(t,rq) > -1;
}


// -------------------------------------------------------------------------
//
// Draw a legend with the fit results
//
void MHCalibrationBlindPixel::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 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 (fFitOK)
    {
      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,3);

    gROOT->SetSelectedPad(NULL);

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

    fHBlindPixelCharge->Draw(opt);
    
    if (fFitOK)
      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->Modified();
    c->Update();
    
    c->cd(3);
    gPad->SetLogy(1);
    gPad->SetBorderMode(0);
    fHBlindPixelTime->Draw(opt);
    
    c->cd(4);

    fHBlindPixelChargevsN->Draw(opt);

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

    c->cd(5);
    
    //    fHBlindPixelPSD = (TH1F*)MFFT::PowerSpectrumDensity(fHBlindPixelChargevsN);
    //    TH1F *hist = MFFT::PowerSpectrumDensity((*this)->fHBlindPixelChargevsN);

    //    fHBlindPixelPSD->Draw(opt);
    c->Modified();
    c->Update();

}



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)
{
  
  //
  // 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.5;
  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 = 20.;
  const Double_t mu_1_guess = mu_0_guess + 50.;
  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  = 0.0025;

  //
  // 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.,1.);

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

      fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));

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

      fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
      fSinglePheFit->SetParLimits(5,norm-0.5,norm+0.5);
      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:
      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,0.,0.005);    
      fSinglePheFit->SetParLimits(5,min,(max-min)/1.5);
      fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
      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 << inf << "Chisquare: " << fChisquare << endl;
  *fLog << inf << "DoF: " << fNdf << endl;
  *fLog << inf << "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 << err << "ERROR: Fit Probability " << fProb 
            << " is smaller than the allowed value: " << gkProbLimit << endl;
      fFitOK = kFALSE;
      return kFALSE;
    }

  Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
  
  if (contSinglePhe < 50.) 
    {
      *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe 
            << " in the Single Photo-Electron peak " << endl;
      fFitOK = kFALSE;
      return kFALSE;
    } 
  else
    *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
  
  fFitOK = kTRUE;
    
  return kTRUE;
}

 
void MHCalibrationBlindPixel::CutAllEdges()
{

  Int_t nbins = 30;

  CutEdges(fHBlindPixelCharge,nbins);

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

  CutEdges(fHBlindPixelChargevsN,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 << inf << "Chisquare: "   << fTimeGausFit->GetChisquare() << endl;
  *fLog << inf << "Ndf: "         << fTimeGausFit->GetNDF() << endl;

  return kTRUE;

}
