/* ======================================================================== *\
!
! *
! * 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), fGeom(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 << "MRawEvtData not found... aborting." << endl;
      return kFALSE;
    }

    const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!runheader)
      *fLog << warn << "Warning - cannot check file type, MRawRunHeader not found." << endl;
    else
      if (runheader->IsMonteCarloRun())
          return kTRUE;
    
    fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
    if (!fCam)
        return kFALSE;

    fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
    if (!fPINDiode)
        return kFALSE;

    fCam->SetPINDiode(fPINDiode);

    fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
    if (!fBlindPixel)
        return kFALSE;

    fCam->SetBlindPixel(fBlindPixel);

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

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

    fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
    if (!fSignals)
    {
        *fLog << err << "MExtractedSignalCam not found... 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 << "MRawRunHeader not found... aborting." << endl;
      return kFALSE;
    }

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

    fCam->SetGeomCam(fGeom);

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

    UInt_t npixels = fGeom->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());
    }
    
    if (fExcludedPixelsFile.IsNull())
        return kTRUE;

    //
    // Look for file to exclude pixels from analysis
    //
    gSystem->ExpandPathName(fExcludedPixelsFile);

    //
    // Initialize reading the file
    //
    ifstream in(fExcludedPixelsFile);
    if (!in)
    {
        *fLog << warn << "Cannot open file '" << fExcludedPixelsFile << "'" << endl;
        return kTRUE;
    }

    *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile << "'" << endl;

    //
    // Read the file and count the number of entries
    //
    UInt_t pixel = 0;

    while (++fNumExcludedPixels)
    {
        in >> pixel;
        if (!in)
            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 << "'" << " is empty " << endl;
    else
        fCam->SetNumPixelsExcluded(fNumExcludedPixels);

    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;
  Float_t avinnerped  = 0;
  Float_t avinnerprms = 0;
  Float_t avinnernum  = 0;
  Float_t avouterped  = 0;
  Float_t avouterprms = 0;
  Float_t avouternum  = 0;

  for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
    {

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

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

      if (ped == -1.) 
	  pix.SetExcluded();
      
      //
      // Check if the pixel has been excluded from the fits
      //
      if (pix.IsExcluded())
        continue;

      if (fGeom->GetPixRatio(pixid) == 1.)
      {
	  avinnerped  += ped;
	  avinnerprms += prms;
	  avinnernum++;
      }
      else
      {
	  avouterped  += ped;
	  avouterprms += prms;
	  avouternum++;
      }
      //
      // 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;
  }

  MCalibrationChargePix *avinnerpix = fCam->GetAverageInnerPix();
  MCalibrationChargePix *avouterpix = fCam->GetAverageOuterPix();
  //
  // set the pedestans in the calibration camera
  //
  if (avinnerpix->IsHiGainSaturation())
  {
      avinnerpix->SetPedestal(avinnerped/avinnernum  * fNumLoGainSamples,
			 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumLoGainSamples),
			 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumLoGainSamples/avinnernum));
      avinnerpix->SetNumLoGainSamples((Float_t)fNumLoGainSamples);
      avinnerpix->ApplyLoGainConversion();
  }
  else
  {
      avinnerpix->SetPedestal(avinnerped/avinnernum  * fNumHiGainSamples,
			 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumHiGainSamples),
			 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumHiGainSamples/avinnernum));
  }

  if (avouterpix->IsHiGainSaturation())
  {
      avouterpix->SetPedestal(avouterped/avouternum  * fNumLoGainSamples,
			      avouterprms/avouternum * TMath::Sqrt((Float_t)fNumLoGainSamples),
			      avouterprms/avouternum * TMath::Sqrt((Float_t)fNumLoGainSamples/avouternum));
      avouterpix->SetNumLoGainSamples((Float_t)fNumLoGainSamples);
      avouterpix->ApplyLoGainConversion();
  }
  else
  {
      avouterpix->SetPedestal(avouterped/avouternum  * fNumHiGainSamples,
			      avouterprms/avouternum * TMath::Sqrt((Float_t)fNumHiGainSamples),
			      avouterprms/avouternum * TMath::Sqrt((Float_t)fNumHiGainSamples/avouternum));
  }
  
  if (!avinnerpix->CheckChargeValidity() || !avinnerpix->CheckTimeValidity())
      if (!avinnerpix->CalcReducedSigma())
	  avinnerpix->CalcFFactorMethod();

  if (!avouterpix->CheckChargeValidity() || !avouterpix->CheckTimeValidity())
      if (!avouterpix->CalcReducedSigma())
	  avouterpix->CalcFFactorMethod();


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