/* ======================================================================== *\
!
! *
! * 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 <TMath.h>
#include <TPad.h>

#include <TMinuit.h>
#include <TFitter.h>

#include <TF1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TPaveText.h>
#include <TRandom.h>

#include <TText.h>

#include "MBinning.h"
#include "MParList.h"

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

ClassImp(MHCalibrationBlindPixel);

using namespace std;
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title)
    :     fFitLegend(NULL), fFitOK(kFALSE), 
          fLambda(0.), fMu0(0.), fMu1(0.), fSigma0(0.), fSigma1(0.),
          fLambdaErr(0.), fMu0Err(0.), fMu1Err(0.), fSigma0Err(0.), fSigma1Err(0.),    
          fChisquare(0.), fProb(0.),  fNdf(0),
          fMeanTime(0.),  fMeanTimeErr(0.), fSigmaTime(0.), fSigmaTimeErr(0.),
          fLambdaCheck(0.), fLambdaCheckErr(0.),
          fFitFunc(kEPoisson4)
{

    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 = -1000.;
    fBlindPixelChargelast  = gkStartBlindPixelBinNr;
    fBlindPixelChargenbins = gkStartBlindPixelBinNr+(int)fBlindPixelChargefirst;

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

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

    fHBlindPixelTime = new TH1I("HBlindPixelTime","Distribution of Mean Arrival Times",nbins,tfirst,tlast);
    fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
    fHBlindPixelTime->SetYTitle("Nr. of events");
    fHBlindPixelTime->Sumw2();
    fHBlindPixelTime->SetDirectory(NULL);

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

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

}

MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
{

  if (fFitLegend)
    delete fFitLegend;

  delete fHBlindPixelCharge;
  delete fHBlindPixelTime;
  delete fHBlindPixelChargevsN;
  
}



void MHCalibrationBlindPixel::ResetBin(Int_t i)
{
    fHBlindPixelCharge->SetBinContent (i, 1.e-20);
    fHBlindPixelTime->SetBinContent(i, 1.e-20);
}


// -------------------------------------------------------------------------
//
// 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 *opt) const
{

  gROOT->SetSelectedPad(NULL);

  TObject *newobj = Clone();

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

  //  ((MHCalibrationBlindPixel*)newobj)->Draw(opt);
  newobj->Draw(opt);
  return newobj;
}
  


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

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

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

    c->Divide(2,2);

    gROOT->SetSelectedPad(NULL);

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

    fHBlindPixelCharge->DrawCopy(opt);
    
    if (fFitOK)
      fSinglePheFit.SetLineColor(kGreen);          
    else
      fSinglePheFit.SetLineColor(kRed);
    
    fSinglePheFit.DrawCopy("same");
    c->Modified();
    c->Update();

    fSinglePhePedFit.SetLineColor(kBlue);
    fSinglePhePedFit.DrawCopy("same");

    c->cd(2);
    DrawLegend();
    c->Modified();
    c->Update();
    
    c->cd(3);
    gPad->SetLogy(1);
    gPad->SetBorderMode(0);
    fHBlindPixelTime->DrawCopy(opt);
    
    if (fHBlindPixelTime->GetFunction("GausTime"))
      {
        TF1 *tfit = fHBlindPixelTime->GetFunction("GausTime");
        if (tfit->GetProb() < 0.01)
          tfit->SetLineColor(kRed);
        else
          tfit->SetLineColor(kGreen);
        
        tfit->DrawCopy("same");
        c->Modified();
        c->Update();
      }
    
    c->cd(4);

    fHBlindPixelChargevsN->DrawCopy(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;
    }

  TF1 simulateSinglePhe;
  
  InitFit(simulateSinglePhe,fBlindPixelChargefirst,fBlindPixelChargelast);
  
  for (Int_t i=0;i<10000; i++) 
    fHBlindPixelCharge->Fill(simulateSinglePhe.GetRandom());
  
  return kTRUE;
}

void MHCalibrationBlindPixel::InitFit(TF1& f, 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 mu_0_guess = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
  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;
  const Double_t norm       = entries/gkSq2Pi;

  //
  // Initialize
  //
  switch (fFitFunc)
    {
      
    case kEPoisson4:
      f = TF1("SinglePheFit",&fPoissonKto4,min,max,6);
      f.SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
      f.SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
      f.SetParLimits(0,0.,1.);
      f.SetParLimits(1,min,(max-min)/1.5);
      f.SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
      f.SetParLimits(3,1.0,(max-min)/2.0);
      f.SetParLimits(4,1.0,(max-min)/2.5);
      f.SetParLimits(5,norm-0.1,norm+0.1);
      break;
      
    case kEPoisson5:
      f = TF1("SinglePheFit",&fPoissonKto5,min,max,6);
      f.SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
      f.SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
      f.SetParLimits(0,0.,1.);
      f.SetParLimits(1,min,(max-min)/1.5);
      f.SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
      f.SetParLimits(3,1.0,(max-min)/2.0);
      f.SetParLimits(4,1.0,(max-min)/2.5);
      f.SetParLimits(5,norm-0.1,norm+0.1);
      break;

    case kEPoisson6:
      f = TF1("SinglePheFit",&fPoissonKto6,min,max,6);
      f.SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
      f.SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
      f.SetParLimits(0,0.,1.);
      f.SetParLimits(1,min,(max-min)/1.5);
      f.SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
      f.SetParLimits(3,1.0,(max-min)/2.0);
      f.SetParLimits(4,1.0,(max-min)/2.5);
      f.SetParLimits(5,norm-0.1,norm+0.1);
      break;

    case kEPoisson7:
      break;

    case kEPolya:
      break;

    case kEMichele:
      break;
      
    }
}

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

  InitFit(fSinglePheFit,rmin,rmax);

  fHBlindPixelCharge->Fit(&fSinglePheFit,opt);
  
  ExitFit(fSinglePheFit);

  fProb      = fSinglePheFit.GetProb();
  fChisquare = fSinglePheFit.GetChisquare();
  fNdf       = fSinglePheFit.GetNDF();

  // Perform the cross-check fitting only the pedestal:
  fSinglePhePedFit = TF1("GausPed","gaus",rmin,fMu0);
  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 100 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;
    }

  if (fProb < 0.01) 
    *fLog << warn << "WARNING: Fit Probability " << fProb << " is smaller than 1% " << endl;

  Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
  
  if (contSinglePhe < 100.) 
    {
      *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);
  fBlindPixelChargenbins = nbins;

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

}
