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

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//   MExtractSignal                                                         //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MExtractSignal.h"

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

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

#include "MGeomCam.h"

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

#include "MParList.h"

#include "MRawEvtData.h" 
#include "MRawEvtPixelIter.h"

#include "TMath.h"

ClassImp(MExtractSignal);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
MExtractSignal::MExtractSignal(const char *name, const char *title)
  : fSaturationLimit(254)
{

    fName  = name  ? name  : "MExtractSignal";
    fTitle = title ? title : "Task to extract the signal from the FADC slices";

    AddToBranchList("MRawEvtData.*");

    SetRange();
}

void MExtractSignal::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
{

    fNumHiGainSamples = hilast-hifirst+1;
    fNumLoGainSamples = lolast-lofirst+1;

    fHiGainFirst = hifirst;
    fLoGainFirst = lofirst;

    fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
    fSqrtLoGainSamples = TMath::Sqrt((Float_t)fNumLoGainSamples);
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MPedestalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MExtractedSignalCam
//
Int_t MExtractSignal::PreProcess(MParList *pList)
{
    fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData"));
    if (!fRawEvt)
    {
      *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl;
        return kFALSE;
    }


    fSignals = (MExtractedSignalCam*)pList->FindCreateObj(AddSerialNumber("MExtractedSignalCam"));
    if (!fSignals)
        return kFALSE;

    fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainFirst+fNumHiGainSamples-1,
                                fLoGainFirst, fLoGainFirst+fNumLoGainSamples-1);

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


// --------------------------------------------------------------------------
//
// The ReInit searches for the following input containers:
//  - MRawRunHeader
//
Bool_t MExtractSignal::ReInit(MParList *pList )
{
    MGeomCam *cam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!cam)
    {
        *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
        return kFALSE;
    }

    // FIXME: This solution may change in the future, MExtractedSignal
    //        must be similar to MCerPhotEvt not to MPedestalCam
    //        (Have to think about the mean size of both solutions)
    fSignals->InitSize(cam->GetNumPixels());    

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Calculate the integral of the FADC time slices and store them as a new
// pixel in the MExtractedSignalCam container.
//
Int_t MExtractSignal::Process()
{
    MRawEvtPixelIter pixel(fRawEvt);
    fSignals->Clear();

    UInt_t satlo=0;

    while (pixel.Next())
    {
        UShort_t satHi = 0;
        UShort_t satLo = 0;

        Byte_t *ptr   = pixel.GetHiGainSamples();
        Byte_t *first = ptr + fHiGainFirst;
        Byte_t *last  = ptr + fHiGainFirst + fNumHiGainSamples;

	UInt_t  sumHi   = 0;
	UInt_t  sumLo   = 0;

        const Int_t pixid = pixel.GetPixelId();
        const MPedestalPix  &ped = (*fPedestals)[pixid]; 
	MExtractedSignalPix &pix = (*fSignals)[pixid];

        const Float_t pedes   = ped.GetPedestal();
        const Float_t pedrms  = ped.GetPedestalRms();

	Byte_t maxhi  = 0;
	Byte_t midhi  = 0;

	for (ptr=first;ptr<last;ptr++)
	  {
	    
	    if (*ptr > maxhi)
	      {
		maxhi = *ptr;
		midhi = ptr-first;
	      }

	    sumHi += *ptr;
	    
	    if (*ptr >= fSaturationLimit)
	      satHi++;

	  }

	ptr   = pixel.GetLoGainSamples();
        first = ptr + fLoGainFirst;
        last  = ptr + fLoGainFirst + fNumLoGainSamples;
	Byte_t maxlo  = 0;
	Byte_t midlo  = 0;

	for (ptr=first;ptr<last;ptr++)
	  {

	    if (*ptr > maxlo)
	      {
		maxlo = *ptr;
		midlo = ptr-first;
	      }

	    sumLo += *ptr;

            if (*ptr >= fSaturationLimit)
                satLo++;
	  }

        if (satLo)
            satlo++;

	pix.SetExtractedSignal((Float_t)sumHi - pedes*(Float_t)fNumHiGainSamples,
			       pedrms*fSqrtHiGainSamples,
			       ((Float_t)sumLo - pedes*(Float_t)fNumLoGainSamples),
			       pedrms*fSqrtLoGainSamples);

	pix.SetGainSaturation(satHi, satHi, satLo);

        // FIXME: Arrival time has to be stored in MArrivalTime!
	if (satHi)
	  pix.SetArrivalTimes((Float_t)(midlo+fLoGainFirst),0.,0.);
	else
	  pix.SetArrivalTimes((Float_t)(midhi+fHiGainFirst),0.,0.);

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

    if (satlo)
        *fLog << warn << "WARNING - Lo Gain saturated in " << satlo << " pixels." << endl;


    fSignals->SetReadyToSave();

    return kTRUE;
}
