/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 2014 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2014
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MTreatSaturation
//
//////////////////////////////////////////////////////////////////////////////
#include "MTreatSaturation.h"

#include <algorithm>

#include <TRandom.h>

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

#include "MParList.h"

#include "MRawRunHeader.h"
#include "MExtralgoSpline.h"

#include "MArrivalTimeCam.h"
#include "MArrivalTimePix.h"

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

#include "MPedestalSubtractedEvt.h"

ClassImp(MTreatSaturation);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MTreatSaturation::MTreatSaturation(const char *name, const char *title)
    : fRaw(0), fEvt(0)
{
    fName  = name  ? name  : "MTreatSaturation";
    fTitle = title ? title : "Base class for replacing saturation with pulses";
}

// --------------------------------------------------------------------------
//
Int_t MTreatSaturation::PreProcess(MParList *pList)
{
    fRaw = (MRawEvtData*)pList->FindObject("MRawEvtData");
    if (!fRaw)
    {
        *fLog << err << "MRawEvtData not found... aborting." << endl;
        return kFALSE;
    }

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

// --------------------------------------------------------------------------
//
Int_t MTreatSaturation::Process()
{
    const UInt_t fHiGainFirst =  50;
    const UInt_t fHiGainLast  = 225;

    const UInt_t npix = fEvt->GetNumPixels();

    const UInt_t rangehi = fHiGainLast - fHiGainFirst - 1;

    if (fDev1.GetSize()!=rangehi)
    {
        fDev1.Set(rangehi);
        fDev2.Set(rangehi);
    }

    for (UInt_t ch=0; ch<npix; ch++)
    {
        Float_t *ptr = fEvt->GetSamples(ch);

        const Float_t *pbeg = ptr+fHiGainFirst+1;
        const Float_t *pend = ptr+fHiGainLast-1;

        // Find maximum
        const Float_t *pmax = pbeg;

        for (const Float_t *p=pbeg; p<pend; p++)
            if (*p>*pmax)
                pmax = p;

        if (*pmax<1800)
            continue;

        MExtralgoSpline s(pbeg, rangehi, fDev1.GetArray(), fDev2.GetArray());
        s.SetExtractionType(MExtralgoSpline::kAmplitudeRel);

        const double leading   = s.SearchYdn(1800., pmax-pbeg);
        const double falling   = s.SearchYup(1800., pmax-pbeg);
        const double width     = falling-leading;
        const double amplitude = 1800/(0.898417 - 0.0187633*width + 0.000163919*width*width - 6.87369e-7*width*width*width + 1.13264e-9*width*width*width*width);
        const double time      = leading-(4.46944 - 0.277272*width + 0.0035479*width*width);

        double baseline = 0;
        for (const Float_t *p=ptr+10; p<ptr+50; p++)
            baseline += *ptr;
        baseline /= 40;

        const UInt_t beg = TMath::CeilNint(leading);
        const UInt_t end = TMath::FloorNint(falling);

        ptr += fHiGainFirst+1;
        for (UInt_t i=0; i<end-beg; i++)
        {
            cout << i << " " << ptr[beg-i] << " ";
            ptr[beg+i] = amplitude*(1-1/(1+exp((i-time)/2.14)))*exp((time-i)/38.8)+baseline;
            cout << ptr[beg+i] << endl;
        }

        cout << endl;
    }

    return kTRUE;
}

