/* ======================================================================== *\
!
! *
! * 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.       //
//                                                                          //
//  Input Containers:                                                       //
//   MRawEvtData                                                            //
//                                                                          //
//  Output Containers:                                                      //
//   MCalibrationCam                                                        //
//                                                                          //
//                                                                          //
//  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) //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MCalibrationCalc.h"
#include "MCalibrationConfig.h"
#include "MCalibrationFits.h"

#include "MCalibrationCam.h"
#include "MCalibrationPix.h"
#include "MCalibrationBlindPix.h"
#include "MCalibrationPINDiode.h"

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

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

#include "MParList.h"
#include "MH.h"

#include "MRawRunHeader.h"
#include "MRawEvtData.h"       // MRawEvtData::GetNumPixels
#include "MRawEvtPixelIter.h"

#include "MTime.h"
#include "TMath.h"

ClassImp(MCalibrationCalc);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. b is the number of slices before the maximum slice,
// a the number of slices behind the maximum slice which is taken as signal.
//
MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
    : 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);
}

// --------------------------------------------------------------------------
//
// 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)
{

    fHistOverFlow = 0;
    fEvents       = 0;
    fCosmics      = 0;

    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
    if (!fRawEvt)
    {
        *fLog << 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);            
      }



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

    fNumHiGainSamples =  fRunHeader->GetNumSamplesHiGain();
    fNumLoGainSamples =  fRunHeader->GetNumSamplesLoGain();

    //
    // FIXME: The next statement simply does not work: 
    //        fRawEvt->GetNumPixels() returns always 0
    // fCalibrations->InitSize(fRawEvt->GetNumPixels());
    //

    fCalibrations->InitSize(577);    

    for (Int_t i=0;i<577;i++)
      {
        MCalibrationPix &pix = (*fCalibrations)[i];
        pix.DefinePixId(i);
      }
    
    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()
{

    fEvents++;

    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 Int_t pixid = pixel.GetPixelId();

        Int_t sum = pixel.GetSumHiGainSamples();

        MPedestalPix    &ped = (*fPedestals)[pixid];
 
        Float_t pedes  = ped.GetPedestal();
        Float_t pedrms = ped.GetPedestalRms();
        
	if ((float)sum < (pedes*fNumHiGainSamples)+(1.5*pedrms) )
           cosmicpix++;
     }


    if (cosmicpix > 50.)
      {
	fCosmics++;
	return kTRUE;
      }

    pixel.Reset();

    while (pixel.Next())
      {

        UShort_t sat = 0;
        UShort_t lowgainoverflow = 0;

        const Int_t pixid = pixel.GetPixelId();

        Byte_t *ptr = pixel.GetHiGainSamples();
        Byte_t mid  = pixel.GetIdxMaxHiGainSample();
        UInt_t max  = pixel.GetMaxHiGainSample();

        Int_t sum = (max > gkSaturationLimit              // overflow ?
                     ? sat++, pixel.GetSumLoGainSamples() // take Low Gain
                     : pixel.GetSumHiGainSamples());      // no overflow

        if (sat)
          {

            ptr = pixel.GetLoGainSamples();
            max = pixel.GetMaxLoGainSample();
            mid = pixel.GetIdxMaxLoGainSample();

            // 
            //  FIXME: It seems the conversion HiGain LoGain is already
            //         performed in the data?!?
            //
            sum = (max > gkSaturationLimit                // overflow of LoGain ??? -> GimmeABreak!!!
                   ? lowgainoverflow++, gkLoGainOverFlow  // OUCH (Florian was maybe right)
                   : sum*gkConversionHiLo    );          // OUFF (Florian was wrong) !! 
                   //                   : sum  );

            if (lowgainoverflow)
		{ 
            	  *fLog << err << dbginf << "Warning: Saturation of LoGain reached in pixel: " << pixid << " " 
                	<< "   sum = " << sum << endl;
	            fHistOverFlow++;
		}
          }

        MPedestalPix    &ped = (*fPedestals)[pixid];
        MCalibrationPix &pix = (*fCalibrations)[pixid];
        
        Float_t pedes  = ped.GetPedestal();
        Float_t pedrms = ped.GetPedestalRms();
        
        //
        // FIXME: This is preliminary, we will change to pedestals per slice!!!
        // Assume pedestals per time slice ==> multiply with number of slices
        //
        pedes  *= (sat ? fNumLoGainSamples : fNumHiGainSamples );
        pedrms *= (sat ? fNumLoGainSamples : fNumHiGainSamples );

        Float_t rsum      = (float)sum - pedes;
        
        switch(pixid)
          {
            
          case gkCalibrationBlindPixelId:

	    // 
	    // FIXME: This works only when the blind pixel ID is much larger than 
	    //        the rest of the pixels (which is the case right now)
	    //
	    if (!blindpixel.FillCharge(sum)) 
	      *fLog << warn << 
		"Overflow or Underflow occurred filling Blind Pixel sum = " << sum << endl;

	    if (!blindpixel.FillTime((int)mid)) 
	      *fLog << warn << 
		"Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl;
	    
	    if (!blindpixel.FillRChargevsTime(rsum,fEvents))
	      *fLog << warn << 
		"Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
            
          case gkCalibrationPINDiodeId:
            if (!pindiode.FillCharge(sum)) 
              *fLog << warn << 
                "Overflow or Underflow occurred filling HCharge: means = " << sum << endl;
            if (!pindiode.FillTime((int)mid)) 
              *fLog << warn << 
                "Overflow or Underflow occurred filling HTime: time = " << (int)mid << endl;
            if (!pindiode.FillRChargevsTime(rsum,fEvents))
              *fLog << warn << 
                "Overflow or Underflow occurred filling HChargevsN: eventnr = " << fEvents << endl;

          default:

            if (!pix.FillCharge(sum))
              *fLog << warn << "Could not fill Charge of pixel: " << pixid 
                    << " signal = " << sum << endl;

            //
            // Fill the reduced charge into the control histo for better visibility
            //
            if (!pix.FillRChargevsTime(rsum,fEvents))
              *fLog << warn << "Could not fill red. Charge vs. EvtNr of pixel: " << pixid 
                    << " signal = " << rsum  << " event Nr: " << fEvents << endl;


            if (!pix.FillTime((int)mid)) 
            *fLog << warn << "Could not fill Time of pixel: " << pixid << " time = " << (int)mid << endl;


          } /* switch(pixid) */

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

    return kTRUE;
}

Int_t MCalibrationCalc::PostProcess()
{

  *fLog << inf << endl;
  *fLog << GetDescriptor() << " Cut Histogram Edges" << endl;
  //
  // Cut edges to make fits and viewing of the hists easier  
  //
  fCalibrations->CutEdges();

  //
  // Get pointer to blind pixel
  //
  MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());

  *fLog << GetDescriptor() << " Fitting the Blind Pixel" << endl;

  // 
  // Fit the blind pixel
  //
  if (TESTBIT(fFlags,kUseBlindPixelFit))
    {
      if (blindpixel.FitCharge())
        if (!fCalibrations->CalcNrPhotInnerPixel())
          *fLog << err << dbginf << "Could not calculate Number of photons from the blind pixel " << endl;
      else
        *fLog << err << dbginf << "Could not fit the blind pixel " << endl;
      
      if (!blindpixel.FitTime())
        *fLog << warn << "Could not the Times of the blind pixel " << endl;

      blindpixel.Draw();

    }

  *fLog << 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];

      const Float_t ped    = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
      const Float_t prms   = (*fPedestals)[pixid].GetPedestalRms() * TMath::Sqrt((float)fNumHiGainSamples);

      pix.SetPedestal(ped,prms);

      if (TESTBIT(fFlags,kUseTimeFits))
        pix.FitTime();
      
      pix.FitCharge();
    }

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

  return kTRUE;
}
