/* ======================================================================== *\
!
! *
! * 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   02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                               
// MHCalibrationChargeCam                                               
//
// Fills the extracted signals of MExtractedSignalCam into the MH-classes: 
//
// MHCalibrationChargeHiGainPix, MHCalibrationChargeLoGainPix for each pixel 
// and additionally for an average of the inner and the outer pixels, respectively.
// 
// By default, subsequently the hi-gain classes are treated unless 
// more than fNumHiGainSaturationLimit (default: 1%) of the overall FADC 
// slices show saturation. In that case, the low-gain classes are treated. 
// If more than fNumLoGainSaturationLimit (default: 1%) of the overall 
// low-gain FADC slices saturate, the pixel is declared not valid and no further 
// treatment is pursued. 
// 
// The filled histograms are fitted to a Gaussian and the mean and sigma with 
// its errors and the fit probability are extracted. If none of these values are 
// NaN's and if the probability is bigger than fProbLimit (default: 0.5%), the fit 
// is declared valid. Otherwise, the fit is repeated within ranges of the previous mean 
// +- 5 sigma. In case that this does not help, the histogram means and RMS's are taken directly,
// but the flag kFitValid is set to FALSE. Outliers of more than fPickUpLimit (default: 5) sigmas 
// from the mean are counted as PickUp events. 
//
// Additionally, the slice number with the highest value is stored and a corresponding 
// histogram is filled. This histogram serves only for a rough cross-check if the 
// signal does not lie at the edges of chose extraction window. 
//
// The class also fills arrays with the signal vs. event number, creates a fourier 
// spectrum out of it and investigates if the projected frequencies follow an exponential 
// distribution. In case that the probability of the exponential fit is less than 
// fProbLimit, the pixel is declared HiGainOscillating or LoGainOscillating, respectively.
// 
// The results are written into MCalibrationChargeCam.
//                                                               
/////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationChargeCam.h"

#include <TVirtualPad.h>
#include <TCanvas.h>
#include <TPad.h>

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

#include "MParList.h"

#include "MHCalibrationChargeHiGainPix.h"
#include "MHCalibrationChargeLoGainPix.h"
#include "MHCalibrationChargePix.h"

#include "MCalibrationChargeCam.h"
#include "MCalibrationChargePix.h"

#include "MGeomCam.h"
#include "MGeomPix.h"

#include "MRawEvtData.h"
#include "MRawEvtPixelIter.h"

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

ClassImp(MHCalibrationChargeCam);

using namespace std;

const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.01;
const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.01;
//
//
//
MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title)
{
    fName  = name  ? name  : "MHCalibrationChargeCam";
    fTitle = title ? title : "Class to fill the calibration histograms ";

    fHiGainArray = new TObjArray;
    fHiGainArray->SetOwner();
    
    fLoGainArray = new TObjArray;
    fLoGainArray->SetOwner();

    fAverageHiGainInnerPix = new MHCalibrationChargeHiGainPix("AverageHiGainInnerPix","Average HiGain FADC sums of inner pixels");
    fAverageLoGainInnerPix = new MHCalibrationChargeLoGainPix("AverageLoGainInnerPix","Average LoGain FADC sums of inner pixels");
    fAverageHiGainOuterPix = new MHCalibrationChargeHiGainPix("AverageHiGainOuterPix","Average HiGain FADC sums of outer pixels");
    fAverageLoGainOuterPix = new MHCalibrationChargeLoGainPix("AverageLoGainOuterPix","Average LoGain FADC sums of outer pixels");

/*
    fAverageHiGainInnerPix->GetHGausHist()->SetName("HCalibrationChargeAverageInnerHiGainPix");
    fAverageHiGainOuterPix->GetHGausHist()->SetName("HCalibrationChargeAverageOuterHiGainPix");
    fAverageLoGainInnerPix->GetHGausHist()->SetName("HCalibrationChargeAverageInnerLoGainPix");
    fAverageLoGainOuterPix->GetHGausHist()->SetName("HCalibrationChargeAverageOuterLoGainPix");
*/
    fAverageHiGainInnerPix->GetHGausHist()->SetTitle("Summed FADC slices average Inner pixels HiGain");
    fAverageLoGainInnerPix->GetHGausHist()->SetTitle("Summed FADC slices average Inner pixels LoGain");
    fAverageHiGainOuterPix->GetHGausHist()->SetTitle("Summed FADC slices average Outer pixels HiGain");
    fAverageLoGainOuterPix->GetHGausHist()->SetTitle("Summed FADC slices average Outer pixels LoGain");

    fAverageHiGainInnerPix->GetHAbsTime()->SetTitle("Absolute Arrival Time average Inner pixels HiGain");
    fAverageLoGainInnerPix->GetHAbsTime()->SetTitle("Absolute Arrival Time average Inner pixels LoGain");
    fAverageHiGainOuterPix->GetHAbsTime()->SetTitle("Absolute Arrival Time average Outer pixels HiGain");
    fAverageLoGainOuterPix->GetHAbsTime()->SetTitle("Absolute Arrival Time average Outer pixels LoGain");

    fAverageHiGainInnerPix->SetChargeLast(1000.);
    fAverageLoGainInnerPix->SetChargeLast(100.);
    fAverageHiGainOuterPix->SetChargeLast(1000.);
    fAverageLoGainOuterPix->SetChargeLast(100.);

    fAverageHiGainInnerPix->SetChargeNbins(4000);
    fAverageLoGainInnerPix->SetChargeNbins(4000);
    fAverageHiGainOuterPix->SetChargeNbins(4000);
    fAverageLoGainOuterPix->SetChargeNbins(4000);

    SetNumHiGainSaturationLimit();
    SetNumLoGainSaturationLimit();

    fNumInnerPixels = 0;
    fNumOuterPixels = 0; 
}

// --------------------------------------------------------------------------
//
// Delete the TClonesArray of MHCalibrationChargePix containers
// Delete the MHCalibrationPINDiode and the MHCalibrationBlindPix
//
// Delete the histograms if they exist
//
MHCalibrationChargeCam::~MHCalibrationChargeCam()
{
  delete fHiGainArray;
  delete fLoGainArray;

  delete fAverageHiGainInnerPix;
  delete fAverageLoGainInnerPix;
  delete fAverageHiGainOuterPix;
  delete fAverageLoGainOuterPix;
}

// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
MHCalibrationChargeHiGainPix &MHCalibrationChargeCam::operator[](UInt_t i)
{
  return *static_cast<MHCalibrationChargeHiGainPix*>(fHiGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
const MHCalibrationChargeHiGainPix &MHCalibrationChargeCam::operator[](UInt_t i) const
{
  return *static_cast<MHCalibrationChargeHiGainPix*>(fHiGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
MHCalibrationChargeLoGainPix &MHCalibrationChargeCam::operator()(UInt_t i)
{
  return *static_cast<MHCalibrationChargeLoGainPix*>(fLoGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
const MHCalibrationChargeLoGainPix &MHCalibrationChargeCam::operator()(UInt_t i) const
{
  return *static_cast<MHCalibrationChargeLoGainPix*>(fLoGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Our own clone function is necessary since root 3.01/06 or Mars 0.4
// I don't know the reason
//
TObject *MHCalibrationChargeCam::Clone(const char *) const
{

    const Int_t nhi = fHiGainArray->GetSize();
    const Int_t nlo = fLoGainArray->GetSize();
    //
    // FIXME, this might be done faster and more elegant, by direct copy.
    //
    MHCalibrationChargeCam *cam = new MHCalibrationChargeCam;
//    MHCalibrationChargeCam cam;
    Copy(*cam);

/*
    cam->fHiGainArray->Expand(nhi);
    cam->fLoGainArray->Expand(nlo);

    for (int i=0; i<nhi; i++)
    {
        delete (*cam->fHiGainArray)[i];
        (*cam->fHiGainArray)[i] = (*fHiGainArray)[i]->Clone();
    }
    for (int i=0; i<nlo; i++)
    {
        delete (*cam->fLoGainArray)[i];
        (*cam->fLoGainArray)[i] = (*fLoGainArray)[i]->Clone();
    }

    delete cam->fAverageHiGainInnerPix;
    delete cam->fAverageLoGainInnerPix;
    delete cam->fAverageHiGainOuterPix;
    delete cam->fAverageLoGainOuterPix;

    cam->GetAverageHiGainInnerPix() = *(MHCalibrationChargeHiGainPix*)fAverageHiGainInnerPix->Clone();
    cam->GetAverageLoGainInnerPix() = *(MHCalibrationChargeLoGainPix*)fAverageLoGainInnerPix->Clone();

    cam->GetAverageHiGainOuterPix() = *(MHCalibrationChargeHiGainPix*)fAverageHiGainOuterPix->Clone();
    cam->GetAverageLoGainOuterPix() = *(MHCalibrationChargeLoGainPix*)fAverageLoGainOuterPix->Clone();

    fAverageHiGainInnerPix->Copy(cam->GetAverageHiGainInnerPix());
    fAverageLoGainInnerPix->Copy(cam->GetAverageLoGainInnerPix());
    fAverageHiGainOuterPix->Copy(cam->GetAverageHiGainOuterPix());
    fAverageLoGainOuterPix->Copy(cam->GetAverageLoGainOuterPix());
*/

    return cam;
}

// --------------------------------------------------------------------------
//
// To setup the object we get the number of pixels from a MGeomCam object
// in the Parameter list.
//
Bool_t MHCalibrationChargeCam::SetupFill(const MParList *pList)
{
  
  fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");

  if (!fRawEvt)
  {
      gLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
      return kFALSE;
  }

  fGeom = (MGeomCam*)pList->FindObject("MGeomCam");

  if (!fGeom)
  {
      gLog << err << dbginf << "MGeomCam not found... aborting." << endl;
      return kFALSE;
  }

  fHiGainArray->Delete();
  fLoGainArray->Delete();

  fAverageHiGainInnerPix->Init();
  fAverageLoGainInnerPix->Init();  
  fAverageHiGainOuterPix->Init();
  fAverageLoGainOuterPix->Init();  

  return kTRUE;
}

Bool_t MHCalibrationChargeCam::ReInit(MParList *pList)
{

  fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
  if (!fCam)
  {
      gLog << err << GetDescriptor() << ": ERROR: Could not find MCalibrationChargeCam ... aborting " << endl;
      return kFALSE;
  }

  MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
  if (!signal)
    {
      gLog << err << "No argument in MExtractedSignalCam::ReInit... abort." << endl;
      return kFALSE;
    }

  const Int_t n = signal->GetSize();

  if (fHiGainArray->GetEntries()==0)
  {
      fHiGainArray->Expand(n);
  
      for (Int_t i=0; i<n; i++)
      {
	  (*fHiGainArray)[i] = new MHCalibrationChargeHiGainPix;
	  (*this)[i].Init();
	  (*this)[i].ChangeHistId(i);
	  
      }
  }

  if (fLoGainArray->GetEntries()==0)
  {
      fLoGainArray->Expand(n);
      
      for (Int_t i=0; i<n; i++)
      {
	  (*fLoGainArray)[i] = new MHCalibrationChargeLoGainPix;
	  (*this)(i).Init();
	  (*this)(i).ChangeHistId(i);
      }
  }

  return kTRUE;
}
  

// --------------------------------------------------------------------------
Bool_t MHCalibrationChargeCam::Fill(const MParContainer *par, const Stat_t w)
{

  MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
  if (!signal)
    {
      gLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
      return kFALSE;
    }
  
  Int_t n = signal->GetSize();

  if (fHiGainArray->GetEntries() != n)
    {
      gLog << err << "ERROR - Size mismatch... abort." << endl;
      return kFALSE;
    }

  if (fLoGainArray->GetEntries() != n)
    {
      gLog << err << "ERROR - Size mismatch... abort." << endl;
      return kFALSE;
    }

  //
  // First the extracted signal
  // 
  Float_t sumhiinnertot = 0.;
  Float_t sumloinnertot = 0.;
  Float_t sumhioutertot = 0.;
  Float_t sumlooutertot = 0.;
  Int_t   sumhiinnersat = 0;
  Int_t   sumloinnersat = 0;
  Int_t   sumhioutersat = 0;
  Int_t   sumlooutersat = 0;

  for (int i=0; i<n; i++)
    {

      const MExtractedSignalPix &pix = (*signal)[i];
      
      const Float_t sumhi  = pix.GetExtractedSignalHiGain();
      const Float_t sumlo  = pix.GetExtractedSignalLoGain();
      
      (*this)[i].FillHistAndArray(sumhi);
      (*this)(i).FillHistAndArray(sumlo);

      const Int_t sathi = (Int_t)pix.GetNumHiGainSaturated();
      const Int_t satlo = (Int_t)pix.GetNumLoGainSaturated();

      (*this)[i].SetSaturated(sathi); 
      (*this)(i).SetSaturated(satlo); 

      if (fGeom->GetPixRatio(i) == 1.)
      {
	  sumhiinnertot += sumhi;
	  sumloinnertot += sumlo;
	  sumhiinnersat += sathi;
	  sumloinnersat += satlo;
	  fNumInnerPixels++;
      }
      else
      {
	  sumhioutertot += sumhi;
	  sumlooutertot += sumlo;
	  sumhioutersat += sathi;
	  sumlooutersat += satlo;
	  fNumOuterPixels++;
      }

    }

  fAverageHiGainInnerPix->FillHistAndArray(sumhiinnertot/fNumInnerPixels); 
  fAverageLoGainInnerPix->FillHistAndArray(sumloinnertot/fNumInnerPixels); 
  fAverageHiGainOuterPix->FillHistAndArray(sumhioutertot/fNumOuterPixels); 
  fAverageLoGainOuterPix->FillHistAndArray(sumlooutertot/fNumOuterPixels); 

  fAverageHiGainInnerPix->SetSaturated(sumhiinnersat/fNumInnerPixels); 
  fAverageLoGainInnerPix->SetSaturated(sumloinnersat/fNumInnerPixels); 
  fAverageHiGainOuterPix->SetSaturated(sumhioutersat/fNumOuterPixels); 
  fAverageLoGainOuterPix->SetSaturated(sumlooutersat/fNumOuterPixels); 

  //
  // Now the times
  //
  sumhiinnertot = sumloinnertot = sumhioutertot = sumlooutertot = 0.;
  MRawEvtPixelIter pixel(fRawEvt);
  while (pixel.Next())
    {
      
      const UInt_t pixid = pixel.GetPixelId();
      
      const Float_t timehi = (Float_t)pixel.GetIdxMaxHiGainSample();
      const Float_t timelo = (Float_t)pixel.GetIdxMaxLoGainSample();

      (*this)[pixid].FillAbsTime(timehi);
      (*this)(pixid).FillAbsTime(timelo+15.);

      if (fGeom->GetPixRatio(pixel.GetPixelId()) == 1.)
      {
	  sumhiinnertot += timehi;
	  sumloinnertot += timelo;
      }
      else
      {
	  sumhioutertot += timehi;
	  sumlooutertot += timelo;
      }
    }

  fAverageHiGainInnerPix-> FillAbsTime(sumhiinnertot/fNumInnerPixels); 
  fAverageLoGainInnerPix-> FillAbsTime(sumloinnertot/fNumInnerPixels); 
  fAverageHiGainOuterPix-> FillAbsTime(sumhioutertot/fNumOuterPixels); 
  fAverageLoGainOuterPix-> FillAbsTime(sumlooutertot/fNumOuterPixels); 

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// 1) Return if the charge distribution is already succesfully fitted  
//    or if the histogram is empty
// 2) Fit the histograms with a Gaussian
// 3) In case of failure set the bit kFitted to false and take histogram means and RMS
// 4) Check for pickup noise
// 5) Check the fourier spectrum 
// 5) Retrieve the results and store them in this class
//
Bool_t MHCalibrationChargeCam::Finalize()
{
    
    for (int i=0; i<fHiGainArray->GetSize(); i++)
    {
     
	MHCalibrationChargeHiGainPix &histhi = (*this)[i];
	MCalibrationChargePix        &pix    = (*fCam)[i];

	FinalizeHiGainHists(histhi,pix);
    }

    for (int i=0; i<fLoGainArray->GetSize(); i++)
    {
     
	MHCalibrationChargeLoGainPix &histlo = (*this)(i);
	MCalibrationChargePix        &pix    = (*fCam)[i];

	FinalizeLoGainHists(histlo,pix);
    }

    FinalizeHiGainHists(*fAverageHiGainInnerPix,*fCam->GetAverageInnerPix());
    FinalizeLoGainHists(*fAverageLoGainInnerPix,*fCam->GetAverageInnerPix());
    FinalizeHiGainHists(*fAverageHiGainOuterPix,*fCam->GetAverageOuterPix());
    FinalizeLoGainHists(*fAverageLoGainOuterPix,*fCam->GetAverageOuterPix());

    fCam->GetAverageInnerPix()->SetSigmaCharge   (fCam->GetAverageInnerPix()->GetSigmaCharge()   *TMath::Sqrt((Float_t)fNumInnerPixels));
    fCam->GetAverageOuterPix()->SetSigmaCharge   (fCam->GetAverageOuterPix()->GetSigmaCharge()   *TMath::Sqrt((Float_t)fNumOuterPixels));
    fCam->GetAverageInnerPix()->SetSigmaChargeErr(fCam->GetAverageInnerPix()->GetSigmaChargeErr()*TMath::Sqrt((Float_t)fNumInnerPixels));
    fCam->GetAverageOuterPix()->SetSigmaChargeErr(fCam->GetAverageOuterPix()->GetSigmaChargeErr()*TMath::Sqrt((Float_t)fNumOuterPixels));

    return kTRUE;
}

void MHCalibrationChargeCam::FinalizeHiGainHists(MHCalibrationChargeHiGainPix &hist, MCalibrationChargePix &pix)
{

    if (hist.IsEmpty())
	return;
    
    if (hist.GetSaturated() > fNumHiGainSaturationLimit*hist.GetHGausHist()->GetEntries())
    {
	pix.SetHiGainSaturation();
	return;
    }
    
    //
    // 2) Fit the Hi Gain histograms with a Gaussian
    //
    if (hist.FitGaus())
	pix.SetHiGainFitted();
    //
    // 3) In case of failure set the bit kFitted to false and take histogram means and RMS
    //
    else if (hist.RepeatFit())
	pix.SetHiGainFitted();
    else
	hist.BypassFit();
    
    //
    // 4) Check for pickup
    //
    hist.CountPickup();
    
    //
    // 5) Check for oscillations
    // 
    hist.CreateFourierSpectrum();
    
    //
    // 6) Retrieve the results and store them in this class
    //
    pix.SetHiGainMeanCharge(     hist.GetMean()     );
    pix.SetHiGainMeanChargeErr(  hist.GetMeanErr()  );
    pix.SetHiGainSigmaCharge(    hist.GetSigma()    );
    pix.SetHiGainSigmaChargeErr( hist.GetSigmaErr() );
    pix.SetHiGainChargeProb    ( hist.GetProb()     );
    
    pix.SetAbsTimeMean         ( hist.GetAbsTimeMean());
    pix.SetAbsTimeRms          ( hist.GetAbsTimeRms() );
    
    pix.SetHiGainOscillating   ( !hist.IsFourierSpectrumOK() );
    pix.SetHiGainNumPickup     (  hist.GetPickup()           );
    
}


void MHCalibrationChargeCam::FinalizeLoGainHists(MHCalibrationChargeLoGainPix &hist, MCalibrationChargePix &pix)
{


    if (hist.IsEmpty())
            return;

    if (hist.GetSaturated() > fNumLoGainSaturationLimit*hist.GetHGausHist()->GetEntries())
    {
	pix.SetLoGainSaturation();
	return;
    }
    //
    // 2) Fit the Lo Gain histograms with a Gaussian
    //
    if (hist.FitGaus())
	pix.SetLoGainFitted();
    //
    // 3) In case of failure set the bit kFitted to false and take histogram means and RMS
    //
    else if (hist.RepeatFit())
	pix.SetLoGainFitted();
    else
	hist.BypassFit();
    
    //
    // 4) Check for pickup
    //
    hist.CountPickup();
    
    //
    // 5) Check for oscillations
    // 
    hist.CreateFourierSpectrum();
    
    //
    // 6) Retrieve the results and store them in this class
    //
    pix.SetLoGainMeanCharge(     hist.GetMean()     );
    pix.SetLoGainMeanChargeErr(  hist.GetMeanErr()  );
    pix.SetLoGainSigmaCharge(    hist.GetSigma()    );
    pix.SetLoGainSigmaChargeErr( hist.GetSigmaErr() );
    pix.SetLoGainChargeProb    ( hist.GetProb()     );
    
    if (pix.IsHiGainSaturation())
    {
	pix.SetAbsTimeMean     ( hist.GetAbsTimeMean());
	pix.SetAbsTimeRms      ( hist.GetAbsTimeRms() );
    }	    
    
    pix.SetLoGainOscillating   ( !hist.IsFourierSpectrumOK() );
    pix.SetLoGainNumPickup     (  hist.GetPickup()           );
    
}    

Bool_t MHCalibrationChargeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// What MHCamera needs in order to draw an individual pixel in the camera
//
void MHCalibrationChargeCam::DrawPixelContent(Int_t idx) const
{

    if (idx == -1)
	fAverageHiGainInnerPix->DrawClone();
    if (idx == -2)
	fAverageHiGainOuterPix->DrawClone();
    if (idx == -3)
	fAverageLoGainInnerPix->DrawClone();
    if (idx == -4)
	fAverageLoGainOuterPix->DrawClone();
    (*this)[idx].DrawClone();

}


void MHCalibrationChargeCam::Draw(const Option_t *opt)
{

  TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);  
  pad->SetBorderMode(0);

  pad->Divide(2,2);

  pad->cd(1);

  if (!fAverageHiGainInnerPix->IsEmpty())
    gPad->SetLogy();
  fAverageHiGainInnerPix->Draw(opt);

  pad->cd(2);

  if (!fAverageLoGainInnerPix->IsEmpty())
    gPad->SetLogy();
  fAverageLoGainInnerPix->Draw(opt);

  pad->cd(3);

  if (!fAverageHiGainOuterPix->IsEmpty())
    gPad->SetLogy();
  fAverageHiGainOuterPix->Draw(opt);

  pad->cd(4);

  if (!fAverageLoGainOuterPix->IsEmpty())
    gPad->SetLogy();
  fAverageLoGainOuterPix->Draw(opt);

}