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

//////////////////////////////////////////////////////////////////////////////
//
//   MCalibrationCalc
//
//   This is a task which calculates the number of photons from the FADC
//   time slices. At the moment it integrates simply the FADC values.
//
//
//  The class MCalibrationCam hold one entry of type MCalibrationPix for
//  every pixel. It is filled in the following way:
//  PreParocess: MalibrationCam::InitSize(577) is called which allocates
//               memory in an TClonesArray of type MCalibrationPix and
//               all pointers to NULL.
//
//  Process:     The NULL pointer is tested on every pixel via
//               MalibrationCam::IsPixelUsed(npix).
//
//               In case, IsPixelUsed returns NULL,
//               MalibrationCam::AddPixel(npix) is invoked which creates a
//               new MCalibrationPix(npix) in the npix's entry
//               of the TClonesArray.
//
//               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 histogram HBlindPixelTime (blind pixel) is fitted to a
//               Gaussian
//               The histograms of the PIN Diode are fitted to Gaussians
//
//               Fits can be excluded via the commands:
//               MalibrationCam::SetSkipTimeFits()   (skip all time fits)
//               MalibrationCam::SetSkipBlindPixelFits()  (skip all blind
//               pixel fits)
//               MalibrationCam::SetSkipPinDiodeFits()  (skip all PIN Diode
//               fits)
//
//  Input Containers:
//   MRawEvtData
//
//  Output Containers:
//   MCalibrationCam
//
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrationCalc.h"

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

// FXIME: This has to be removed!!!!
#include "MCalibrationConfig.h"

#include <TSystem.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 "MCalibrationCam.h"
#include "MCalibrationPix.h"

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

#include "MCalibrationBlindPix.h"
#include "MCalibrationPINDiode.h"

ClassImp(MCalibrationCalc);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
    : fPedestals(NULL), fCalibrations(NULL), fSignals(NULL),
      fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL),
      fEvents(0), fHistOverFlow(0), fCosmics(0),
      fNumHiGainSamples(0), fNumLoGainSamples(0), fConversionHiLo(0.),
      fNumExcludedPixels(0),
      fColor(kEBlue)
{

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

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

    SETBIT(fFlags, kUseTimeFits);
    SETBIT(fFlags, kUseBlindPixelFit);
    SETBIT(fFlags, kUsePinDiodeFit);

}

MCalibrationBlindPix *MCalibrationCalc::GetBlindPixel() const
{
    return fCalibrations->GetBlindPixel();
}

MCalibrationPINDiode *MCalibrationCalc::GetPINDiode() const
{
    return fCalibrations->GetPINDiode();
}

// --------------------------------------------------------------------------
//
// 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
//  - MTime
//
Int_t MCalibrationCalc::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;
        }
    
    fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
    if (!fCalibrations)
      {
        *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;        
        return kFALSE;
      }


    switch (fColor)
      {
      case kEBlue:
        fCalibrations->SetColor(MCalibrationCam::kECBlue);
	break;        
      case kEGreen:
        fCalibrations->SetColor(MCalibrationCam::kECGreen);      
	break;
      case kEUV:
        fCalibrations->SetColor(MCalibrationCam::kECUV);            
	break;
      case kECT1:
        fCalibrations->SetColor(MCalibrationCam::kECCT1);            
	break;
      default:
        fCalibrations->SetColor(MCalibrationCam::kECCT1); 
      }

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


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

    UInt_t npixels = cam->GetNumPixels();

    fCalibrations->InitSize(npixels);    
    
    for (UInt_t i=0; i<npixels; i++)
      {
        
        MCalibrationPix &pix = (*fCalibrations)[i];
        pix.DefinePixId(i);
        MHCalibrationPixel *hist = pix.GetHist();
        
        hist->SetTimeFitRangesHiGain(fSignals->GetFirstUsedSliceHiGain(),
                                     fSignals->GetLastUsedSliceHiGain());
        hist->SetTimeFitRangesLoGain(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
                //
                MCalibrationPix &pix = (*fCalibrations)[pixel];
                pix.SetExcluded();
                
                *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
                
              }
            
            if (--fNumExcludedPixels == 0)
              *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data() 
                    << "'" << " is empty " << endl;
            else
              fCalibrations->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 MCalibrationCalc::Process()
{

    Int_t cosmicpix = 0;

    MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
    MCalibrationPINDiode &pindiode   = *(fCalibrations->GetPINDiode());

    MRawEvtPixelIter pixel(fRawEvt);

    //
    // Create a first loop to sort out the cosmics ...
    // 
    // This is a very primitive check for the number of cosmicpixs
    // The cut will be applied in the fit, but for the blind pixel,
    // we need to remove this event
    //
    // FIXME: In the future need a much more sophisticated one!!!
    //

    while (pixel.Next())
      {
	
	const UInt_t pixid = pixel.GetPixelId();
	
	MExtractedSignalPix &sig =  (*fSignals)[pixid];
        MPedestalPix        &ped =  (*fPedestals)[pixid];
        Float_t pedrms           = ped.GetPedestalRms()*fSqrtHiGainSamples;
	Float_t sumhi            = sig.GetExtractedSignalHiGain();
        
        //
        // We consider a pixel as presumably due to cosmics 
        // if its sum of FADC slices is lower than 3 pedestal RMS
        //
	if (sumhi < 3.*pedrms )  
	  cosmicpix++;
     }

    //
    // If the camera contains more than 230 
    // (this is the number of outer pixels plus about 50 inner ones)
    // presumed pixels due to cosmics, then the event is discarted. 
    // This procedure is more or less equivalent to keeping only events 
    // with at least 350 pixels with high signals.
    //
    if (cosmicpix > 230.)
      {
	fCosmics++;
	return kCONTINUE;
      }

    pixel.Reset();
    fEvents++;


    //
    // Create a second loop to do fill the calibration histograms
    // 

    while (pixel.Next())
      {

	const UInt_t pixid = pixel.GetPixelId();
	
	MCalibrationPix &pix = (*fCalibrations)[pixid];

        if (pix.IsExcluded())
          continue;

	MExtractedSignalPix &sig =  (*fSignals)[pixid];
        
	Float_t sumhi  = sig.GetExtractedSignalHiGain();
	Float_t sumlo  = sig.GetExtractedSignalLoGain();
	Float_t mtime  = sig.GetMeanArrivalTime();


        switch(pixid)
          {
            
          case gkCalibrationBlindPixelId:

	    if (!blindpixel.FillCharge(sumhi)) 
	      *fLog << warn << 
		"Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl;

	    if (!blindpixel.FillTime((int)mtime)) 
	      *fLog << warn << 
		"Overflow or Underflow occurred filling Blind Pixel time = " << mtime << endl;
	    
	    if (!blindpixel.FillRChargevsTime(sumhi,fEvents))
	      *fLog << warn << 
		"Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
            break;
            
          case gkCalibrationPINDiodeId:
            if (!pindiode.FillCharge(sumhi)) 
              *fLog << warn << 
                "Overflow or Underflow occurred filling PINDiode: sum = " << sumhi << endl;
            if (!pindiode.FillTime((int)mtime)) 
              *fLog << warn << 
                "Overflow or Underflow occurred filling PINDiode: time = " << mtime << endl;
            if (!pindiode.FillRChargevsTime(sumhi,fEvents))
              *fLog << warn << 
                "Overflow or Underflow occurred filling PINDiode: eventnr = " << fEvents << endl;
	    break;

          default:

	    pix.SetChargesInGraph(sumhi,sumlo);

	    if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
	      *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: " 
		    << pixid << " signal = " << sumlo  << " event Nr: " << fEvents << endl;
		
	    if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
	      *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: " 
		    << pixid << " signal = " << sumhi  << " event Nr: " << fEvents << endl;

	    if (sig.IsLoGainUsed())
	      {
		
		if (!pix.FillChargeLoGain(sumlo))
		  *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid 
			<< " signal = " << sumlo << endl;

		if (!pix.FillTimeLoGain((int)mtime)) 
		  *fLog << warn << "Could not fill Lo Gain Time of pixel: " 
			<< pixid << " time = " << mtime << endl;
		
	      }
	    else
	      {
		if (!pix.FillChargeHiGain(sumhi))
		  *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid 
			<< " signal = " << sumhi << endl;
		
		if (!pix.FillTimeHiGain((int)mtime)) 
		  *fLog << warn << "Could not fill Hi Gain Time of pixel: " 
			<< pixid << " time = " << mtime << endl;
		
	      }
	    break;
	    
          } /* switch(pixid) */

      } /* while (pixel.Next()) */

    return kTRUE;
}

Int_t MCalibrationCalc::PostProcess()
{

  *fLog << inf << endl;

  if (fEvents == 0)
    {
      *fLog << err << GetDescriptor() 
            << ": This run contains only cosmics or pedestals, " 
            << "cannot find events with more than 350 illuminated pixels. " << endl;
      return kFALSE;
    }
  
  if (fEvents < fCosmics)
    *fLog << warn << GetDescriptor() 
          << ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl;


  *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;

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

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

      if (!blindpixel.FitCharge())
        {
          *fLog << err << dbginf << "Could not fit the blind pixel! " << endl;
          blindpixel.DrawClone();
          return kFALSE;
        }

      //blindpixel.DrawClone();
    }
  else 
    *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << 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 = (*fCalibrations)[pixid];

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

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


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

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

      //
      // Perform the Gauss fits to the arrival times
      //
      if (TESTBIT(fFlags,kUseTimeFits))
        pix.FitTime();
      
    }

  if (TESTBIT(fFlags,kUseBlindPixelFit))
    {
      
      if (!fCalibrations->CalcNumPhotInsidePlexiglass())
        {
          *fLog << err 
                << "Could not calculate the number of photons from the blind pixel " << endl;
          *fLog << err 
                << "You can try to calibrate using the MCalibrationCalc::SkipBlindPixelFit()" << endl;
          return kFALSE;
        }
    }
  else
    *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;

  fCalibrations->SetReadyToSave();
  
  if (GetNumExecutions()==0)
    return kTRUE;
  
  *fLog << inf << endl;
  *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;

  return kTRUE;
}