/* ======================================================================== *\
!
! *
! * 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 necessary fits to extract the mean number of photons   //
//              out of the derived light flux                               //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationBlindPixel.h"
#include "MHCalibrationConfig.h"
#include "MCalibrationFits.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 "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)
  : fSinglePheFit(NULL), fTimeGausFit(NULL), fSinglePhePedFit(NULL),
    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.)
{

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

    fErrBlindPixelChargefirst = 0.;
    fErrBlindPixelChargelast  = gkStartBlindPixelBinNr;
    fErrBlindPixelChargenbins = gkStartBlindPixelBinNr;

    fHBlindPixelErrCharge = new TH1F("HBlindPixelErrCharge","Distribution of Variances of Summed FADC Slices",
                        fErrBlindPixelChargenbins,fErrBlindPixelChargefirst,fErrBlindPixelChargelast);
    fHBlindPixelErrCharge->SetXTitle("Variance Summed FADC Slices");
    fHBlindPixelErrCharge->SetYTitle("Nr. of events");
    fHBlindPixelErrCharge->Sumw2();

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

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

    fgSinglePheFitFunc = &gfKto8;
    fgSinglePheFitNPar = 5;
}

MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
{

  delete fHBlindPixelCharge;
  delete fHBlindPixelTime;
  delete fHBlindPixelErrCharge;
  
  if (fSinglePheFit)
    delete fSinglePheFit;
  if (fSinglePhePedFit)
    delete fSinglePhePedFit;
  if (fTimeGausFit)
    delete fTimeGausFit;
}



void MHCalibrationBlindPixel::ResetBin(Int_t i)
{
    fHBlindPixelCharge->SetBinContent (i, 1.e-20);
    fHBlindPixelErrCharge->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());
  fFitLegend->AddText(line1);

  const TString line6 =
  Form("Mean #lambda (check) = %2.2f #pm %2.2f",GetLambdaCheck(),GetLambdaCheckErr());
  fFitLegend->AddText(line6);

  const TString line2 = 
  Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err());
  fFitLegend->AddText(line2);

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

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

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

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

  const TString line8 =
  Form("Probability: %4.2f ",GetProb());
  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 MHCalibrationBlindPixel::Draw(Option_t *opt) 
{

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

    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 (fSinglePheFit)
      {
        if (fFitOK)
	  fSinglePheFit->SetLineColor(kGreen);          
        else
          fSinglePheFit->SetLineColor(kRed);

        fSinglePheFit->DrawCopy("same");
        c->Modified();
        c->Update();

	if (fSinglePhePedFit)
	  {
	    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->GetEntries() != 0)
    {
      *fLog << err << "Histogram " << fHBlindPixelCharge->GetTitle() << " is already filled. " << endl;
      *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
      return kFALSE;
    }
  
  TF1 *simulateSinglePhe = new TF1("simulateSinglePhe",fgSinglePheFitFunc,
                                   fBlindPixelChargefirst,fBlindPixelChargelast,fgSinglePheFitNPar);
  
  simulateSinglePhe->SetParameters(lambda,mu0,mu1,sigma0,sigma1);
  simulateSinglePhe->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1");
  simulateSinglePhe->SetNpx(fBlindPixelChargenbins);  

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


void MHCalibrationBlindPixel::ChangeFitFunc(BlindPixelFitFunc fitfunc, Int_t par)
{
  
  fgSinglePheFitFunc = fitfunc;
  fgSinglePheFitNPar = par;

}



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

  if (fSinglePheFit)
    return kFALSE;


  //
  // Get the fitting ranges
  //
  rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst;
  rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast;

  //
  // 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->GetEntries();
  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 + 100.;
  const Double_t si_1_guess = si_0_guess + si_0_guess;

  fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar+1);
  fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,entries);
  fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","area");
  fSinglePheFit->SetParLimits(0,0.,3.);
  fSinglePheFit->SetParLimits(1,rmin,rmax);
  fSinglePheFit->SetParLimits(2,rmin,rmax);
  fSinglePheFit->SetParLimits(3,1.0,rmax-rmin);
  fSinglePheFit->SetParLimits(4,1.7,rmax-rmin);
  fSinglePheFit->SetParLimits(5,0.,1.5*entries);
  //
  // Normalize the histogram to facilitate faster fitting of the area
  // For speed reasons, FKto8 is normalized to Sqrt(2 pi).
  // Therefore also normalize the histogram to that value 
  //
  // ROOT gives us another nice example of user-unfriendly behavior:
  // Although the normalization of the function fSinglePhe and the 
  // Histogram fHBlindPixelCharge agree (!!), the fit does not normalize correctly INTERNALLY
  // in the fitting procedure !!!
  // 
  // This has to do with the fact that the internal function histogramming 
  // uses 100 bins and does not adapt to the binning of the fitted histogram, unlike PAW does
  // (very important if you use Sumw2, see e.g. ROOTTALK: Mon May 26 1997 - 09:56:03 MEST)
  // 
  // So, WE have to adapt to that internal flaw of ROOT:
  //
  //  const Int_t  npx     = fSinglePheFit->GetNpx();
  //  const Int_t  bins    = fHBlindPixelCharge->GetXaxis()->GetLast()-fHBlindPixelCharge->GetXaxis()->GetFirst();
  //  fHBlindPixelCharge->Scale(gkSq2Pi*(float)bins/npx/entries);

  // 
  // we need this, otherwise, ROOT does not calculate the area correctly
  // don't ask me why it does not behave correctly, it's one of the nasty
  // mysteries of ROOT which takes you a whole day to find out :-)
  //
  //  fSinglePheFit->SetNpx(fChargenbins);  

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

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

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

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

  Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2);
  cout << "Parameter0: " << fSinglePhePedFit->GetParameter(0) << endl;
  cout << "Parameter2: " << fSinglePhePedFit->GetParameter(2) << endl;
  cout << "Pedarea: " << pedarea << endl;
  cout << "entries: " << entries << endl;
  fLambdaCheck     = TMath::Log((double)entries/pedarea);
  fLambdaCheckErr  = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0)
                     + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2);

  *fLog << "Results of the Blind Pixel Fit: " << endl;
  *fLog << "Chisquare: " << fChisquare << endl;
  *fLog << "DoF: " << fNdf << endl;
  *fLog << "Probability: " << fProb << endl;
  *fLog << "Integral: " << fSinglePheFit->Integral(rmin,rmax);

  //
  // The fit result is accepted under condition
  // The fit result is accepted under condition
  // The Probability is greater than gkProbLimit (default 0.001 == 99.7%)
  //
  if (fProb < gkProbLimit) 
    {
      *fLog << warn << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl;
      fFitOK = kFALSE;
      return kFALSE;
    }
  
  
  fFitOK = kTRUE;
    
  return kTRUE;
}

 
void MHCalibrationBlindPixel::CutAllEdges()
{

  //
  // The number 100 is necessary because it is the internal binning
  // of ROOT functions. A call to SetNpx() does NOT help
  // If you find another solution which WORKS!!, please tell me!!
  //
  Int_t nbins = 50;

  *fLog << "New number of bins in HSinCharge: " << CutEdges(fHBlindPixelCharge,nbins) << endl;

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

  *fLog << "New number of bins in HErrCharge: " << CutEdges(fHBlindPixelErrCharge,30) << endl;
  fErrBlindPixelChargefirst = fHBlindPixelErrCharge->GetBinLowEdge(fHBlindPixelErrCharge->GetXaxis()->GetFirst());
  fErrBlindPixelChargelast  = fHBlindPixelErrCharge->GetBinLowEdge(fHBlindPixelErrCharge->GetXaxis()->GetLast())+fHBlindPixelErrCharge->GetBinWidth(0);
  fErrBlindPixelChargenbins = nbins;

  CutEdges(fHBlindPixelChargevsN,0);

}

Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt) 
{
  
  if (fTimeGausFit)
    return kFALSE;

  rmin = (rmin != 0.) ? rmin : 4.;
  rmax = (rmax != 0.) ? rmax : 9.;

  const Stat_t   entries     = fHBlindPixelTime->GetEntries();
  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);

  Float_t prob = fTimeGausFit->GetProb();

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

  if (prob < gkProbLimit) 
    {
      *fLog << warn << "Fit of the Arrival times failed ! " << endl;
      return kFALSE;
    }
  
  return kTRUE;

}
