/* ======================================================================== *\
!
! *
! * 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 "MCalibrationFits.h"

#include <TStyle.h>
#include <TMath.h>

#include <TFitter.h>

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

#include "MParList.h"

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

ClassImp(MHCalibrationPixel);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationPixel::MHCalibrationPixel(const char *name, const char *title)
      : fPixId(-1),
        fChargeGausFit(NULL), 
	fTimeGausFit(NULL), 
	fFitLegend(NULL),
	fLowerFitRange(0.),
        fFitOK(kFALSE),
        fChargeChisquare(-1.),
        fChargeProb(-1.),
        fChargeNdf(-1),
        fTimeChisquare(-1.),
        fTimeProb(-1.),
        fTimeNdf(-1)
{ 

    fName  = name  ? name  : "MHCalibrationPixel";
    fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits";

    // Create a large number of bins, later we will rebin
    fChargeFirst = -0.5;
    fChargeLast  = 10000. - 0.5;
    fChargeNbins = 20000;

    fHCharge = new TH1I("HCharge","Distribution of Summed FADC Slices Pixel ",
                   fChargeNbins,fChargeFirst,fChargeLast);
    fHCharge->SetXTitle("Sum FADC Slices");
    fHCharge->SetYTitle("Nr. of events");
    fHCharge->Sumw2();

    fHCharge->SetDirectory(NULL);

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

    fHTime = new TH1I("HTime","Distribution of Mean Arrival Times Pixel ",
                  ntbins,tfirst,tlast);
    fHTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
    fHTime->SetYTitle("Nr. of events");
    fHTime->Sumw2();

    fHTime->SetDirectory(NULL);

    TString qvsntitle = "Sum of Charges vs. Event Number Pixel ";

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

    fHChargevsN = new TH1I("HChargevsN",qvsntitle.Data(),
                     nqbins,nfirst,nlast);
    fHChargevsN->SetXTitle("Event Nr.");
    fHChargevsN->SetYTitle("Sum of FADC slices");

    fHChargevsN->SetDirectory(NULL);

}

MHCalibrationPixel::~MHCalibrationPixel()
{

  delete fHCharge;
  delete fHTime;
  delete fHChargevsN;

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

}


void MHCalibrationPixel::ChangeHistId(Int_t id)
{

  fPixId = id;
  
  TString nameQ = TString(fHCharge->GetName());
  nameQ += id;
  fHCharge->SetName(nameQ.Data());

  TString nameT = TString(fHTime->GetName());
  nameT += id;
  fHTime->SetName(nameT.Data());

  TString nameQvsN  = TString(fHChargevsN->GetName());
  nameQvsN += id;
  fHChargevsN->SetName(nameQvsN.Data());

  TString titleQ = TString(fHCharge->GetTitle());
  titleQ += id;
  fHCharge->SetTitle(titleQ.Data());

  TString titleT = TString(fHTime->GetTitle());
  titleT += id;
  fHTime->SetTitle(titleT.Data());

  TString titleQvsN  = TString(fHChargevsN->GetTitle());
  titleQvsN += id;
  fHChargevsN->SetTitle(titleQvsN.Data());
}


void MHCalibrationPixel::Reset()
{
  
  for (Int_t i = fHCharge->FindBin(fChargeFirst); 
       i <= fHCharge->FindBin(fChargeLast); i++)
    fHCharge->SetBinContent(i, 1.e-20);

  for (Int_t i = 0; i < 16; i++)
      fHTime->SetBinContent(i, 1.e-20);
  
  fChargeLast     = 9999.5;

  fHCharge->GetXaxis()->SetRangeUser(0.,fChargeLast);

  return;
}


// -------------------------------------------------------------------------
//
// Set the binnings and prepare the filling of the histograms
//
Bool_t MHCalibrationPixel::SetupFill(const MParList *plist)
{

  fHCharge->Reset();
  fHTime->Reset();

  return kTRUE;
}



// -------------------------------------------------------------------------
//
// Draw a legend with the fit results
//
void MHCalibrationPixel::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 Gauss Fit:");
  fFitLegend->SetTextSize(0.05);

  const TString line1 = 
    Form("Mean: Q_{#mu} = %2.2f #pm %2.2f",fChargeMean,fChargeMeanErr);

  fFitLegend->AddText(line1);


  const TString line4 =
    Form("Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fChargeSigma,fChargeSigmaErr);

  fFitLegend->AddText(line4);


  const TString line7 =
    Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChargeChisquare,fChargeNdf);

  fFitLegend->AddText(line7);


  const TString line8 =
    Form("Probability: %4.3f ",fChargeProb);

  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 MHCalibrationPixel::Draw(Option_t *opt) 
{

    gStyle->SetOptFit(0);
    gStyle->SetOptStat(1111111);
    
    TCanvas *c = MakeDefCanvas(this,600,900); 

    gROOT->SetSelectedPad(NULL);

    c->Divide(2,2);

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

    fHCharge->Draw(opt);
    
    if (fChargeGausFit)
      {
        if (fFitOK)
          fChargeGausFit->SetLineColor(kGreen);          
        else
          fChargeGausFit->SetLineColor(kRed);

        fChargeGausFit->Draw("same");
        c->Modified();
        c->Update();
      }

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

    c->cd(3);
    gStyle->SetOptStat(1111111);

    gPad->SetLogy(1);
    fHTime->Draw(opt);

    if (fTimeGausFit)
      {
        if (fTimeChisquare > 1.)
          fTimeGausFit->SetLineColor(kRed);
        else
          fTimeGausFit->SetLineColor(kGreen);

        fTimeGausFit->Draw("same");
        c->Modified();
        c->Update();
      }
    
    c->Modified();
    c->Update();

    c->cd(4);
    fHChargevsN->Draw(opt);
}



Bool_t MHCalibrationPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *option)
{

  if (fTimeGausFit)
    return kFALSE;

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

  const Stat_t   entries     = fHTime->GetEntries();
  const Double_t mu_guess    = fHTime->GetBinCenter(fHTime->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 << err << dbginf << "Could not create fit function for Gauss 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);

  fHTime->Fit(fTimeGausFit,option);

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

  fHTime->Fit(fTimeGausFit,option);

  fTimeChisquare = fTimeGausFit->GetChisquare();
  fTimeNdf       = fTimeGausFit->GetNDF();
  fTimeProb      = fTimeGausFit->GetProb();

  fTimeMean      = fTimeGausFit->GetParameter(1);
  fTimeSigma     = fTimeGausFit->GetParameter(2);

  if (fTimeChisquare > 1.) 
    {
      *fLog << warn << "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 = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirst;
  Axis_t rmax = 0.;

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

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

  fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLast);

  if (!fChargeGausFit) 
    {
    *fLog << err << dbginf << "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,fChargeLast);
  fChargeGausFit->SetParLimits(2,0.,fChargeLast-rmin);
  fChargeGausFit->SetRange(rmin,fChargeLast);

  fHCharge->Fit(fChargeGausFit,option);
  
  Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2);
  
  rmin = (rtry < rmin ? rmin : rtry);
  rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2);
  fChargeGausFit->SetRange(rmin,rmax);  

  fHCharge->Fit(fChargeGausFit,option);

  //  rmin = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2);
  //  rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2);
  //  fChargeGausFit->SetRange(rmin,rmax);  

  // fHCharge->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 Probability is greater than gkProbLimit (default 0.01 == 99%)
  //
  if (fChargeProb < gkProbLimit) 
    {
      *fLog << warn << "Prob: " << fChargeProb << " is smaller than the allowed value: " << gkProbLimit << endl;
      fFitOK = kFALSE;
      return kFALSE;
    }
  
  
  fFitOK = kTRUE;
    
  return kTRUE;
}

 
void MHCalibrationPixel::CutAllEdges()
{

  Int_t nbins = 50;

  CutEdges(fHCharge,nbins);

  fChargeFirst = fHCharge->GetBinLowEdge(fHCharge->GetXaxis()->GetFirst());
  fChargeLast  = fHCharge->GetBinLowEdge(fHCharge->GetXaxis()->GetLast())+fHCharge->GetBinWidth(0);
  fChargeNbins = nbins;

  CutEdges(fHChargevsN,0);

}

void MHCalibrationPixel::PrintChargeFitResult()
{
  
  *fLog << "Results of the Summed Charges Fit: "                      << endl;
  *fLog << "Chisquare: "        << fChargeChisquare                   << endl;
  *fLog << "DoF: "              << fChargeNdf                         << endl;
  *fLog << "Probability: "      << fChargeProb                        << endl;
  *fLog                                                               << endl;
  
}

void MHCalibrationPixel::PrintTimeFitResult()
{

  *fLog << "Results of the Time Slices Fit: "                        << endl;
  *fLog << "Chisquare: "   << fTimeChisquare                         << endl;
  *fLog << "Ndf: "         << fTimeNdf                               << endl;
  *fLog << "Probability: " << fTimeProb                              << endl;
  *fLog                                                              << endl;

}