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

//////////////////////////////////////////////////////////////////////////////
//
//   MCalibrationChargeCalc
//
//   Task to calculate the calibration conversion factors from the FADC
//   time slices. The integrated time slices have to be delivered by an 
//   MExtractedSignalCam. The pedestals by an MPedestalCam.
//
//   The output container MCalibrationCam holds one entry of type MCalibrationPix 
//   for every pixel. It is filled in the following way:
//
//   ProProcess: Search for MPedestalCam, MExtractedSignalCam and MExtractedSignalBlindPixel
//               Initialize MCalibrationCam
//               Initialize pulser light wavelength
//               
//   ReInit:     MCalibrationCam::InitSize(NumPixels) is called which allocates
//               memory in a TClonesArray of type MCalibrationPix
//               Initialize number of used FADC slices
//               Optionally exclude pixels from calibration               
//
//   Process:    Every MCalibrationPix holds a histogram class,
//               MHCalibrationPixel which itself hold histograms of type:
//               HCharge(npix) (distribution of summed FADC time slice
//               entries)
//               HTime(npix) (distribution of position of maximum)
//               HChargevsN(npix) (distribution of charges vs. event number.
//
//  PostProcess:  All histograms HCharge(npix) are fitted to a Gaussian
//                All histograms HTime(npix) are fitted to a Gaussian
//                The histogram HBlindPixelCharge (blind pixel) is fitted to
//                a single PhE fit
//
//                The histograms of the PIN Diode are fitted to Gaussians
//
//                Fits can be excluded via the commands:
//                MalibrationCam::SkipBlindPixelFits()  (skip all blind
//                pixel fits)
//
//                Hi-Gain vs. Lo-Gain Calibration (very memory-intensive)
//                can be skipped with the command:
//                MalibrationCam::SkipHiLoGainCalibration()
//
//  Input Containers:
//   MRawEvtData
//   MPedestalCam
//   MExtractedSignalCam
//   MExtractedSignalBlindPixel
//
//  Output Containers:
//   MCalibrationCam
//
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargeCalc.h"

// FXIME: Usage of fstream is a preliminary workaround!
#include <fstream>

#include <TSystem.h>
#include <TH1.h>

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MRawRunHeader.h"
#include "MRawEvtPixelIter.h"

#include "MPedestalCam.h"
#include "MPedestalPix.h"

#include "MCalibrationChargeCam.h"
#include "MCalibrationChargePINDiode.h"
#include "MCalibrationPix.h"

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

#include "MCalibrationBlindPix.h"

ClassImp(MCalibrationChargeCalc);

using namespace std;

const UInt_t MCalibrationChargeCalc::fgBlindPixelIdx   = 559;
const UInt_t MCalibrationChargeCalc::fgPINDiodeIdx     = 9999;
const UInt_t MCalibrationChargeCalc::fgBlindPixelSinglePheCut = 400;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
    : fPedestals(NULL), fCam(NULL), 
      fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL),
      fSignals(NULL), fBlindPixel(NULL), fPINDiode(NULL)
{

    fName  = name  ? name  : "MCalibrationChargeCalc";
    fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";

    AddToBranchList("MRawEvtData.fHiGainPixId");
    AddToBranchList("MRawEvtData.fLoGainPixId");
    AddToBranchList("MRawEvtData.fHiGainFadcSamples");
    AddToBranchList("MRawEvtData.fLoGainFadcSamples");

    Clear();
    SetBlindPixelIdx();
    SetPINDiodeIdx();
}

void MCalibrationChargeCalc::Clear(const Option_t *o)
{
  
    SETBIT(fFlags, kUseBlindPixelFit);
    SETBIT(fFlags, kUseQualityChecks);
    SETBIT(fFlags, kHiLoGainCalibration);

    CLRBIT(fFlags, kHiGainOverFlow);
    CLRBIT(fFlags, kLoGainOverFlow);

    fNumBlindPixelSinglePhe = 0;
    fNumBlindPixelPedestal  = 0;  

    fNumHiGainSamples  = 0;
    fNumLoGainSamples  = 0;
    fConversionHiLo    = 0;
    fNumExcludedPixels = 0;

}


// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MPedestalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MHCalibrationBlindPixel
//  - MCalibrationCam
//
// The following output containers are only searched, but not created
//
//  - MTime
//
Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
{

    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
    if (!fRawEvt)
    {
      *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
      return kFALSE;
    }

    const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!runheader)
      *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
    else
      if (runheader->GetRunType() == kRTMonteCarlo)
        {
          return kTRUE;
        }
    
    fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
    if (!fCam)
      {
        *fLog << err << dbginf << "MCalibrationChargeCam could not be created ... aborting." << endl;        
        return kFALSE;
      }

    fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
    if (!fPINDiode)
      {
        *fLog << err << dbginf << "MCalibrationChargePINDiode could not be created ... aborting." << endl;        
        return kFALSE;
      }

    fCam->SetPINDiode(fPINDiode);

    fEvtTime      = (MTime*)pList->FindObject("MTime");

    fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
    if (!fPedestals)
      {
        *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
        return kFALSE;
      }


    fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
    if (!fSignals)
      {
        *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
        return kFALSE;
      }

    fBlindPixel = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel");
    if (!fBlindPixel)
      {
        *fLog << err << dbginf << "Cannot find MExtractedSignalBlindPixel ... aborting" << endl;
        return kFALSE;
      }
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// The ReInit searches for the following input containers:
//  - MRawRunHeader
//
Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
{
 
    fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!fRunHeader)
    {
      *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
      return kFALSE;
    }


    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!cam)
    {
      *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
      return kFALSE;
    }

    fCam->SetGeomCam(cam);

    fNumHiGainSamples =  fSignals->GetNumUsedHiGainFADCSlices();
    fNumLoGainSamples =  fSignals->GetNumUsedLoGainFADCSlices();
    fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);

    UInt_t npixels = cam->GetNumPixels();

    for (UInt_t i=0; i<npixels; i++)
      {
        
        MCalibrationPix &pix = (*fCam)[i];
        pix.DefinePixId(i);

        pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
                                    fSignals->GetLastUsedSliceHiGain());
        pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
                                    fSignals->GetLastUsedSliceLoGain());
        
        if (!TESTBIT(fFlags,kUseQualityChecks))
          pix.SetExcludeQualityCheck();

        // Exclude the blind pixel and the PIN Diode from normal pixel calibration:
        if (i == fBlindPixelIdx)
          pix.SetExcluded();

        if (i == fPINDiodeIdx)
          pix.SetExcluded();

     }
    
    //
    // Look for file to exclude pixels from analysis
    //
    if (!fExcludedPixelsFile.IsNull())
      {
        
        fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
        
        //
        // Initialize reading the file
        //
        ifstream in(fExcludedPixelsFile.Data(),ios::in);

        if (in)
          {
            *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
            //
            // Read the file and count the number of entries
            //
            UInt_t pixel = 0;
            
            while (++fNumExcludedPixels)
              {
                
                in >> pixel;

                if (!in.good())
                  break;
                //
                // Check for out of range
                //
                if (pixel > npixels)
                  {
                    *fLog << warn << "WARNING: To be excluded pixel: " << pixel 
                          << " is out of range " << endl;
                    continue;
                  }
                //
                // Exclude pixel
                //
                MCalibrationPix &pix = (*fCam)[pixel];
                pix.SetExcluded();
                
                *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
                
              }
            
            if (--fNumExcludedPixels == 0)
              *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data() 
                    << "'" << " is empty " << endl;
            else
              fCam->SetNumPixelsExcluded(fNumExcludedPixels);
            
          }
        else
          *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
      }
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Calculate the integral of the FADC time slices and store them as a new
// pixel in the MCerPhotEvt container.
//
Int_t MCalibrationChargeCalc::Process()
{

  
  MRawEvtPixelIter pixel(fRawEvt);
  
  //
  // Create a loop to fill the calibration histograms
  // Search for: a signal in MExtractedSignalCam 
  // Search for: a signal in MExtractedSignalBlindPixel
  // Fill histograms with:
  //             charge
  //             charge vs. event nr.
  //
  //
  // Initialize pointers to blind pixel and individual pixels
  //
  // FIXME: filling the bind pixel histograms in this class is preliminary 
  //        and will be replaced soon to fill them with MFillH
  //
  MCalibrationBlindPix       &blindpixel = *(fCam->GetBlindPixel());
  MExtractedSignalBlindPixel &blindsig   = (*fBlindPixel);

  const UInt_t signal = blindsig.GetExtractedSignal();
  
  if (!blindpixel.FillCharge(signal)) 
    *fLog << warn << 
      "Overflow or Underflow occurred filling Blind Pixel sum = " << signal << endl;

  blindpixel.FillGraphs(signal,0);
  
  TH1I *hist;
  
  if (signal > fBlindPixelSinglePheCut)
    {
      hist = (blindpixel.GetHist())->GetHSinglePheFADCSlices();
      fNumBlindPixelSinglePhe++;
    }
  else
    {
      hist = (blindpixel.GetHist())->GetHPedestalFADCSlices();
      fNumBlindPixelPedestal++;
    }

  pixel.Jump(fBlindPixelIdx);
  
  const Byte_t *ptr = pixel.GetHiGainSamples();

  for (Int_t i=1;i<16;i++)
    hist->Fill(i,*ptr++);
  
  ptr = pixel.GetLoGainSamples();
  for (Int_t i=16;i<31;i++)
    hist->Fill(i,*ptr++);

  pixel.Reset();
  
  while (pixel.Next())
    {
      
      const UInt_t pixid = pixel.GetPixelId();
      
      MCalibrationPix            &pix   =  (*fCam)[pixid];
      MExtractedSignalPix        &sig   =  (*fSignals)     [pixid];
      
      const Float_t sumhi  = sig.GetExtractedSignalHiGain();
      const Float_t sumlo  = sig.GetExtractedSignalLoGain();

      Float_t abstime = 0.;

      if (sig.IsLoGainUsed())
        abstime = (Float_t)pixel.GetIdxMaxLoGainSample();
      else
        abstime = (Float_t)pixel.GetIdxMaxHiGainSample();

      if (pix.IsExcluded())
        continue;
      
      pix.FillGraphs(sumhi,sumlo);

      if (sig.IsLoGainUsed())
        {
          
          if (!pix.FillChargeLoGain(sumlo))
            *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid 
                  << " signal = " << sumlo << endl;
          
          if (!pix.FillAbsTimeLoGain(abstime)) 
            *fLog << warn << "Could not fill Lo Gain Abs. Time of pixel: " 
                  << pixid << " time = " << abstime << endl;
        }
      else
        {
          if (!pix.FillChargeHiGain(sumhi))
            *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid 
                  << " signal = " << sumhi << endl;
          
          if (!pix.FillAbsTimeHiGain(abstime))
            *fLog << warn << "Could not fill Hi Gain Abs. Time of pixel: " 
                  << pixid << " time = " << abstime << endl;
        }
      
    } /* while (pixel.Next()) */

  return kTRUE;
}

Int_t MCalibrationChargeCalc::PostProcess()
{
  *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;

  //
  // Cut edges to make fits and viewing of the hists easier  
  //
  fCam->CutEdges();

  // 
  // Fit the blind pixel
  //
  if (TESTBIT(fFlags,kUseBlindPixelFit))
  {
      //
      // Get pointer to blind pixel
      //
      MCalibrationBlindPix &blindpixel = *(fCam->GetBlindPixel());
      
      *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;

      //
      // retrieve the histogram containers
      //
      MHCalibrationBlindPixel *hist = blindpixel.GetHist();
          
      //
      // retrieve mean and sigma of the blind pixel pedestal, 
      // so that we can use it for the fit
      //
      // FIXME: This part has to be migrated into the future MHCalibrationChargeBlindPix class
      // 
      const UInt_t nentries    = fPedestals->GetTotalEntries();
      const UInt_t nslices     = 12;
      const Float_t sqrslice   = TMath::Sqrt((Float_t)nslices);

      MPedestalPix &pedpix  = (*fPedestals)[fBlindPixelIdx];
      if (&pedpix)
        {
          const Float_t pedestal    = pedpix.GetPedestal()*nslices;
          const Float_t pederr      = pedpix.GetPedestalRms()*nslices/nentries;
          const Float_t pedsigma    = pedpix.GetPedestalRms()*sqrslice;
          const Float_t pedsigmaerr = pederr/2.;
          
          hist->SetMeanPedestal(pedestal);
          hist->SetMeanPedestalErr(pederr);
          hist->SetSigmaPedestal(pedsigma);
          hist->SetSigmaPedestalErr(pedsigmaerr);
        }
      
      if (!blindpixel.FitCharge())
        {
          *fLog << warn << "Could not fit the blind pixel! " << endl;
          *fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl;
          fCam->SetBlindPixelMethodValid(kFALSE);
        }
      else
        fCam->SetBlindPixelMethodValid(kTRUE);
      
      if (blindpixel.CheckOscillations())
        fCam->SetBlindPixelMethodValid(kFALSE);

      TH1I *sphehist = hist->GetHSinglePheFADCSlices();
      TH1I *pedhist  = hist->GetHPedestalFADCSlices();

      if (fNumBlindPixelSinglePhe > 1)
        sphehist->Scale(1./fNumBlindPixelSinglePhe);
      if (fNumBlindPixelPedestal > 1)
        pedhist->Scale(1./fNumBlindPixelPedestal);
      
      blindpixel.DrawClone();
  }
  else 
      *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;

  *fLog << inf << "total: " << GetNumExecutions() << " sphe: " << fNumBlindPixelSinglePhe << " ped: " << fNumBlindPixelPedestal << endl;

  
  *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;

  //
  // loop over the pedestal events and check if we have calibration
  //
  for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
    {

      MCalibrationPix &pix = (*fCam)[pixid];

      //
      // Check if the pixel has been excluded from the fits
      //
      if (pix.IsExcluded())
        continue;

      //
      // get the pedestals
      //
      const Float_t ped    = (*fPedestals)[pixid].GetPedestal();
      const Float_t prms   = (*fPedestals)[pixid].GetPedestalRms();

      //
      // set them in the calibration camera
      //
      pix.SetPedestal(ped,prms,(Float_t)fNumHiGainSamples,(Float_t)fNumLoGainSamples);

      //
      // perform the Gauss fits to the charges
      //
      pix.FitCharge();

      //
      // check also for oscillations
      // 
      pix.CheckOscillations();

      //
      // calculate the F-Factor method
      //
      pix.CalcFFactorMethod();
    }

  if (TESTBIT(fFlags,kUseBlindPixelFit) && fCam->IsBlindPixelMethodValid())
  {
      if (!fCam->CalcFluxInsidePlexiglass())
      {
          *fLog << warn << "Could not calculate the number of photons from the blind pixel " << endl;
          *fLog << "You can try to calibrate using the MCalibrationChargeCalc::SkipBlindPixelFit()" << endl;
        fCam->SetBlindPixelMethodValid(kFALSE);          
      }
  }
  else
    *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl;

  if (!fPINDiode->CalcFluxOutsidePlexiglass())
  {
      *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
      fCam->SetPINDiodeMethodValid(kFALSE);
  }
  else
  {
      fCam->SetPINDiodeMethodValid(kTRUE);
  }

  fCam->SetReadyToSave();
  
  return kTRUE;
}