/* ======================================================================== *\
!
! *
! * 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(Int_t pix, const char *name, const char *title)
    : fFitOK(kFALSE), fPixId(pix), fTGausFit(NULL), fQGausFit(NULL), fFitLegend(NULL)
{

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

    TString qname  = "HQ";
    TString qtitle = "Distribution of Summed FADC Slices Pixel ";
    qname  += pix;
    qtitle += pix;

    // Create a large number of bins, later we will rebin
    fQfirst = -0.5;
    fQlast  = gkStartQlast - 0.5;
    fQnbins = gkStartPixelBinNr;

    fHQ = new TH1I( qname.Data(),qtitle.Data(),
                   fQnbins,fQfirst,fQlast);
    fHQ->SetXTitle("Sum FADC Slices");
    fHQ->SetYTitle("Nr. of events");
    fHQ->Sumw2();

    TString tname  = "HT";
    TString ttitle = "Distribution of Mean Arrival Times Pixel ";
    tname  += pix;
    ttitle += pix;

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

    fHT = new TH1I(tname.Data(),ttitle.Data(),
                  nbins,tfirst,tlast);
    fHT->SetXTitle("Mean Arrival Times [FADC slice nr]");
    fHT->SetYTitle("Nr. of events");
    fHT->Sumw2();

    TString qvsnname  = "HQvsN";
    TString qvsntitle = "Sum of Charges vs. Event Number Pixel ";
    qvsnname  += pix;
    qvsntitle += pix;

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

    fHQvsN = new TH1I(qvsnname.Data(),qvsntitle.Data(),
                     nbins,nfirst,nlast);
    fHQvsN->SetXTitle("Event Nr.");
    fHQvsN->SetYTitle("Sum of FADC slices");

    fQChisquare = -1.;
    fQProb      = -1.;
    fQNdf       = -1;
    
    fTChisquare = -1.;
    fTProb      = -1.;
    fTNdf       = -1;

}

MHCalibrationPixel::~MHCalibrationPixel()
{

  delete fHQ;
  delete fHT;
  delete fHQvsN;

  if (fQGausFit)
    delete fQGausFit;
  if (fTGausFit)
    delete fTGausFit;
  if (fFitLegend)
    delete fFitLegend;

}

void MHCalibrationPixel::Reset()
{
  
  for (Int_t i = fHQ->FindBin(fQfirst); i <= fHQ->FindBin(fQlast); i++)
      fHQ->SetBinContent(i, 1.e-20);

  for (Int_t i = 0; i < 16; i++)
      fHT->SetBinContent(i, 1.e-20);
  
  fQlast     = gkStartQlast;

  fHQ->GetXaxis()->SetRangeUser(0.,fQlast);

  return;
}


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

  fHQ->Reset();
  fHT->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);

  char line1[32];
  sprintf(line1,"Mean: Q_{#mu} = %2.2f #pm %2.2f",fQMean,fQMeanErr);
  fFitLegend->AddText(line1);

  char line4[32];
  sprintf(line4,"Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fQSigma,fQSigmaErr);
  fFitLegend->AddText(line4);

  char line7[32];
  sprintf(line7,"#chi^{2} / N_{dof}: %4.2f / %3i",fQChisquare,fQNdf);
  fFitLegend->AddText(line7);

  char line8[32];
  sprintf(line8,"Probability: %4.2f ",fQProb);
  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();

    fHQ->DrawCopy(opt);
    
    if (fQGausFit)
      {
        if (fFitOK)
          fQGausFit->SetLineColor(kGreen);          
        else
          fQGausFit->SetLineColor(kRed);

        fQGausFit->DrawCopy("same");
        c->Modified();
        c->Update();
      }

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

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

    gPad->SetLogy(1);
    fHT->DrawCopy(opt);

    if (fHT->GetFunction("GausTime"))
      {
        TF1 *tfit = fHT->GetFunction("GausTime");
        if (tfit->GetProb() < 0.01)
          tfit->SetLineColor(kRed);
        else
          tfit->SetLineColor(kGreen);

        tfit->DrawCopy("same");
        c->Modified();
        c->Update();
      }
    
    c->Modified();
    c->Update();

    c->cd(4);
    fHQvsN->DrawCopy(opt);
}



Bool_t MHCalibrationPixel::FitT(Axis_t rmin, Axis_t rmax, Option_t *option)
{
  
  if (fTGausFit)
    return kFALSE;

  rmin = (rmin != 0.) ? rmin : -0.5;
  rmax = (rmax != 0.) ? rmax : 15.5;

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

  fTGausFit = new TF1("GausTime","gaus",rmin,rmax);  

  if (!fTGausFit) 
    {
    *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
    return kFALSE;
    }

  fTGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
  fTGausFit->SetParNames("Area","#mu","#sigma");
  fTGausFit->SetParLimits(0,0.,entries);
  fTGausFit->SetParLimits(1,rmin,rmax);
  fTGausFit->SetParLimits(2,0.,rmax-rmin);

  fHT->Fit("GausTime",option);

  fTChisquare = fTGausFit->GetChisquare();
  fTNdf       = fTGausFit->GetNDF();
  fTProb      = fTGausFit->GetProb();
  fTMean      = fTGausFit->GetParameter(1);
  fTSigma     = fTGausFit->GetParameter(2);

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

}

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

  if (fQGausFit)
    return kFALSE;

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

  //
  // 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  = fHQ->GetEntries();
  const Double_t ar_guess = entries/gkSq2Pi;
  const Double_t mu_guess = fHQ->GetBinCenter(fHQ->GetMaximumBin());
  const Double_t si_guess = mu_guess/500.;

  fQGausFit = new TF1("QGausFit","gaus",rmin,rmax);

  if (!fQGausFit) 
    {
    *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
    return kFALSE;
    }
  
  fQGausFit->SetParameters(ar_guess,mu_guess,si_guess);
  fQGausFit->SetParNames("Area","#mu","#sigma");
  fQGausFit->SetParLimits(0,0.,entries);
  fQGausFit->SetParLimits(1,rmin,rmax);
  fQGausFit->SetParLimits(2,0.,rmax-rmin);

  fHQ->Fit("QGausFit",option);

  fQChisquare = fQGausFit->GetChisquare();
  fQNdf       = fQGausFit->GetNDF();
  fQProb      = fQGausFit->GetProb();
  fQMean      = fQGausFit->GetParameter(1);
  fQMeanErr   = fQGausFit->GetParError(1);
  fQSigma     = fQGausFit->GetParameter(2);
  fQSigmaErr  = fQGausFit->GetParError(2);

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

 
void MHCalibrationPixel::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 = 100;

  CutEdges(fHQ,nbins);

  fQfirst = fHQ->GetBinLowEdge(fHQ->GetXaxis()->GetFirst());
  fQlast  = fHQ->GetBinLowEdge(fHQ->GetXaxis()->GetLast())+fHQ->GetBinWidth(0);
  fQnbins = nbins;

  CutEdges(fHQvsN,0);

}

void MHCalibrationPixel::PrintQFitResult()
{
  
  *fLog << "Results of the Summed Charges Fit: "                 << endl;
  *fLog << "Chisquare: "        << fQChisquare                   << endl;
  *fLog << "DoF: "              << fQNdf                         << endl;
  *fLog << "Probability: "      << fQProb                        << endl;
  *fLog                                                          << endl;
  
}

void MHCalibrationPixel::PrintTFitResult()
{

  *fLog << "Results of the Arrival Time Slices Fit: "             << endl;
  *fLog << "Chisquare: "   << fTChisquare                         << endl;
  *fLog << "Ndf: "         << fTNdf                               << endl;
  *fLog << "Probability: " << fTProb                              << endl;
  *fLog                                                           << endl;

}