/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   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 MCalibrationChargePix 
//   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 MCalibrationChargePix
//               Initialize number of used FADC slices
//               Optionally exclude pixels from calibration               
//
//   Process:    Every MCalibrationChargePix 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 "MCalibrationChargePix.h"
#include "MCalibrationChargePINDiode.h"
#include "MCalibrationChargeBlindPix.h"

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


ClassImp(MCalibrationChargeCalc);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
    : fPedestals(NULL), fCam(NULL), 
      fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL),
      fSignals(NULL), fPINDiode(NULL), fBlindPixel(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();
}

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

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

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

    fCam->SetBlindPixel(fBlindPixel);

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

    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++)
      {
        
        MCalibrationChargePix &pix = (*fCam)[i];
        pix.DefinePixId(i);

        pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
                                    fSignals->GetLastUsedSliceHiGain());
        pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
                                    fSignals->GetLastUsedSliceLoGain());
        
     }
    
    //
    // 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
                //
                MCalibrationChargePix &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()
{
  return kTRUE;
}

Int_t MCalibrationChargeCalc::PostProcess()
{

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

      MCalibrationChargePix &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();
      const Float_t num    = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());

      //
      // set them in the calibration camera
      //
      if (pix.IsHiGainSaturation())
      {
	  pix.SetPedestal(ped  * fNumLoGainSamples,
			  prms * TMath::Sqrt((Float_t)fNumLoGainSamples),
			  prms * fNumLoGainSamples / num);
	  pix.SetNumLoGainSamples((Float_t)fNumLoGainSamples);
          pix.ApplyLoGainConversion();
      }
      else
      {
	  pix.SetPedestal(ped  * fNumHiGainSamples,
			  prms * TMath::Sqrt((Float_t)fNumHiGainSamples),
			  prms * fNumHiGainSamples / num);
      }

      if (!pix.CheckChargeValidity() || !pix.CheckTimeValidity())
	  continue;

      nvalid++;

      if (!pix.CalcReducedSigma())
	  continue;

      pix.CalcFFactorMethod();

    }

  //
  // The Michele check ...
  //
  if (nvalid == 0)
  {
      *fLog << err << GetDescriptor() << ": Dear Michele! All pixels have non-valid calibration. " 
	    << "Did you forget to fill the histograms (filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
      return kFALSE;
  }

  if (!fBlindPixel->CheckChargeFitValidity())
  {
      *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
      fCam->SetBlindPixelMethodValid(kFALSE);
  }
  else
  { 
      if (!fBlindPixel->CalcFluxInsidePlexiglass())
      {
	  *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
	  fCam->SetBlindPixelMethodValid(kFALSE);
      }
      else
      {
	  fCam->SetBlindPixelMethodValid(kTRUE);
	  fCam->ApplyBlindPixelCalibration();
      }
  }

  if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
  {
      *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
      fCam->SetPINDiodeMethodValid(kFALSE);
  }
  else
  { 
      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->ApplyPINDiodeCalibration();
      }
  }
  fCam->SetReadyToSave();
  
  return kTRUE;
}