#include <TStyle.h>
#include <TCanvas.h>
#include <TSystem.h>
#include <TF1.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TMath.h>
#include <TGraph.h>
#include <TFitResultPtr.h>
#include <TFitResult.h>
#include <TFile.h>

#include <cstdio>
#include <stdio.h>
#include <stdint.h>

#include "MH.h"
#include "MArrayI.h"
#include "MLog.h"
#include "MLogManip.h"
#include "MDirIter.h"
#include "MFillH.h"
#include "MEvtLoop.h"
#include "MCamEvent.h"
#include "MHCamEvent.h"
#include "MGeomApply.h"
#include "MTaskList.h"
#include "MParList.h"
#include "MContinue.h"
#include "MBinning.h"
#include "MDrsCalibApply.h"
#include "MDrsCalibration.h"
#include "MRawFitsRead.h"
#include "MBadPixelsCam.h"
#include "MStatusArray.h"
#include "MStatusDisplay.h"
#include "MTaskInteractive.h"
#include "MPedestalSubtractedEvt.h"
#include "MHCamera.h"
#include "MGeomCamFAMOUS.h"
#include "MRawRunHeader.h"
#include "MPedestalCam.h"
#include "MPedestalPix.h"
#include "MParameters.h"
#include "PixelMap.h"

using namespace std;
using namespace TMath;

MHCamera *hgain = 0;
MHCamera *hbase = 0;

PixelMap pmap;

// Structure to store the extracted value and the extracted time
struct Single
{
    float    fSignal;
    float    fTime;
};

//Storage class to kKeep a list of the extracted single
class MSingles : public MParContainer, public MCamEvent
{
    Int_t fIntegrationWindow;

    vector<vector<Single>> fData;

public:
    MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30)
    {
        fName = name ? name : "MSingles";
        fName = title ? title : "Storeage for vector of singles";
    }

    void InitSize(const UInt_t i)
    {
        fData.resize(i);
    }

    vector<Single> &operator[](UInt_t i) { return fData[i]; }
    vector<Single> &GetVector(UInt_t i)  { return fData[i]; }

    UInt_t GetNumPixels() const { return fData.size(); }

    void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; }
    Int_t GetIntegrationWindow() const { return fIntegrationWindow; }

    Bool_t GetPixelContent(Double_t &, Int_t , const MGeomCam &, Int_t) const
    {
        return kTRUE;
    }
    void   DrawPixelContent(Int_t) const { }

    ClassDef(MSingles, 1)
};

// Histogram class to extract the baseline
class MHBaseline : public MH
{
    TH2F fBaseline;

    MPedestalCam *fPedestal;

    // The baseline is only extracted where also the signal is extracted
    // FIXME: Make sure this is consistent with MExtractSingles
    UShort_t fSkipStart;
    UShort_t fSkipEnd;

public:
    MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10)
    {
        fName = "MHBaseline";

        // Setup the histogram
        fBaseline.SetName("Baseline");
        fBaseline.SetTitle("Median spectrum");
        fBaseline.SetXTitle("Pixel [idx]");
        fBaseline.SetYTitle("Median baseline [mV]");
        fBaseline.SetDirectory(NULL);
    }

    Bool_t ReInit(MParList *plist)
    {
        fPedestal = (MPedestalCam*)plist->FindCreateObj("MPedestalCam");
        if (!fPedestal)
            return kFALSE;

        const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
        if (!header)
        {
            *fLog << err << "MRawRunHeader not found... abort." << endl;
            return kFALSE;
        }

        // Bin width should be around 1 dac count which is about 0.5mV
        MBinning binsx, binsy;
        binsx.SetEdges(64, -0.5, 64-0.5);
        binsy.SetEdges(100, -20.5, 29.5);

        // Setup binnning
        MH::SetBinning(fBaseline, binsx, binsy);

        return kTRUE;
    }

    // Fill the samples into the histogram
    Int_t Fill(const MParContainer *par, const Stat_t)
    {
        const MPedestalSubtractedEvt *evt = dynamic_cast<const MPedestalSubtractedEvt*>(par);

        const Int_t n = evt->GetNumSamples()-fSkipStart-fSkipEnd;

        // Loop over all pixels
        for (int pix=0; pix<64; pix++)
        {
            // Get samples for each pixel
            const Float_t *ptr = evt->GetSamples(pix);

            // Average two consecutive samples
            for (int i=0; i<n; i+=2)
            {
                const Double_t val = 0.5*ptr[i+fSkipStart]+0.5*ptr[i+1+fSkipStart];
                fBaseline.Fill(pix, val);
            }
        }

        return kTRUE;
    }

    // Extract the baseline value from the distrbutions
    Bool_t Finalize()
    {
        if (!fPedestal)
            return kTRUE;

        fPedestal->InitSize(fBaseline.GetNbinsX());
        fPedestal->SetNumEvents(GetNumExecutions());

        Int_t    cnt = 0;
        Double_t avg = 0;
        Double_t rms = 0;

        // Loop over all 'pixels'
        for (int x=0; x<64; x++)
        {
            // Get the corresponding slice from the histogram
            TH1D *hist = fBaseline.ProjectionY("proj", x+1, x+1);

            // Get the maximum bin
            const Int_t bin  = hist->GetMaximumBin();

            // Calculate a parabola through this and the surrounding points
            // on logarithmic values (that's a gaussian)

            //
            // Quadratic interpolation
            //
            // calculate the parameters of a parabula such that
            //    y(i)  = a + b*x(i) +   c*x(i)^2
            //    y'(i) =     b      + 2*c*x(i)
            //
            //

            // -------------------------------------------------------------------------
            // a = y2;
            // b = (y3-y1)/2;
            // c = (y3+y1)/2 - y2;

            const Double_t v1 = hist->GetBinContent(bin-1);
            const Double_t v2 = hist->GetBinContent(bin);
            const Double_t v3 = hist->GetBinContent(bin+1);
            if (v1<=0 || v2<=0 || v3<=0)
                continue;

            const Double_t y1 = TMath::Log(v1);
            const Double_t y2 = TMath::Log(v2);
            const Double_t y3 = TMath::Log(v3);

            //Double_t a = y2;
            const Double_t b = (y3-y1)/2;
            const Double_t c = (y1+y3)/2 - y2;
            if (c>=0)
                continue;

            const Double_t w  = -1./(2*c);
            const Double_t dx =  b*w;

            if (dx<-1 || dx>1)
                continue;

            // y = exp( - (x-k)^2 / s^2 / 2 )
            //
            // -2*s^2 * log(y) = x^2 - 2*k*x + k^2
            //
            // a = (k/s0)^2/2
            // b = k/s^2
            // c = -1/(2s^2)      -->    s = sqrt(-1/2c)

            const Double_t binx = hist->GetBinCenter(bin);
            const Double_t binw = hist->GetBinWidth(bin);

            //const Double_t p = hist->GetBinCenter(hist->GetMaximumBin());
            const Double_t p = binx + dx*binw;

            avg += p;
            rms += p*p;

            cnt++;

            // Store baseline and sigma
            MPedestalPix &pix = (*fPedestal)[x];

            pix.SetPedestal(p);
            pix.SetPedestalRms(TMath::Sqrt(w)*binw);

            delete hist;
        }

        avg /= cnt;
        rms /= cnt;

        cout << "Baseline(" << cnt << "): " << avg << " +- " << sqrt(rms-avg*avg) << endl;

        return kTRUE;
    }

    // Draw histogram
    void Draw(Option_t *)
    {
        TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);

        AppendPad("");

        pad->SetBorderMode(0);
        pad->SetFrameBorderMode(0);
        fBaseline.Draw("colz");
    }


    ClassDef(MHBaseline, 1);
};

// Histogram class for the signal and time distribution as
// well as the pulse shape
class MHSingles : public MH
{
    TH2F fSignalAll;
    TH2F fSignalOut;
    TH2F fSignal1;
    TH2F fSignal2;
    TH2F fSignal3;
    TH2F fSignal4;
    TH2F fSignal5;
    TH2F fSignal6;
    TH2F fSignal7;
    TH2F fSignal8;

    TProfile fProfAll;
    TProfile fProf1;
    TProfile fProf2;
    TProfile fProf3;
    TProfile fProf4;
    TProfile fProf5;
    TProfile fProf6;
    TProfile fProf7;
    TProfile fProf8;

    UInt_t fNumEntries;

    MSingles               *fSingles;   //!
    MPedestalSubtractedEvt *fEvent;     //!
    MBadPixelsCam          *fBadPix;    //!
    MPedestalCam           *fPedestal;  //!

public:
    MHSingles() : fNumEntries(0), fSingles(0), fEvent(0), fPedestal(0)
    {
        fName = "MHSingles";

        // Setup histograms
        fSignalAll.SetName("SignalAll");
        fSignalAll.SetTitle("");
        fSignalAll.SetXTitle("Time [ns]");
        fSignalAll.SetYTitle("Amplitude [a.u.]");
        fSignalAll.GetXaxis()->CenterTitle();
        fSignalAll.GetXaxis()->SetTitleSize(0.045);
        fSignalAll.GetYaxis()->SetTitleSize(0.045);
        fSignalAll.SetDirectory(NULL);
        fSignalAll.SetStats(kFALSE);

        fSignalOut.SetName("SignalOut");
        fSignalOut.SetTitle("");
        fSignalOut.SetXTitle("Time [ns]");
        fSignalOut.SetYTitle("Amplitude [a.u.]");
        fSignalOut.GetXaxis()->CenterTitle();
        fSignalOut.GetXaxis()->SetTitleSize(0.045);
        fSignalOut.GetYaxis()->SetTitleSize(0.045);
        fSignalOut.SetDirectory(NULL);
        fSignalOut.SetStats(kFALSE);

        fSignal1.SetName("Signal1");
        fSignal1.SetTitle("");
        fSignal1.SetXTitle("Time [ns]");
        fSignal1.SetYTitle("Amplitude [a.u.]");
        fSignal1.GetXaxis()->CenterTitle();
        fSignal1.GetXaxis()->SetTitleSize(0.045);
        fSignal1.GetYaxis()->SetTitleSize(0.045);
        fSignal1.SetDirectory(NULL);
        fSignal1.SetStats(kFALSE);

        fSignal2.SetName("Signal2");
        fSignal2.SetTitle("");
        fSignal2.SetXTitle("Time [ns]");
        fSignal2.SetYTitle("Amplitude [a.u.]");
        fSignal2.GetXaxis()->CenterTitle();
        fSignal2.GetXaxis()->SetTitleSize(0.045);
        fSignal2.GetYaxis()->SetTitleSize(0.045);
        fSignal2.SetDirectory(NULL);
        fSignal2.SetStats(kFALSE);

        fSignal3.SetName("Signal3");
        fSignal3.SetTitle("");
        fSignal3.SetXTitle("Time [ns]");
        fSignal3.SetYTitle("Amplitude [a.u.]");
        fSignal3.GetXaxis()->CenterTitle();
        fSignal3.GetXaxis()->SetTitleSize(0.045);
        fSignal3.GetYaxis()->SetTitleSize(0.045);
        fSignal3.SetDirectory(NULL);
        fSignal3.SetStats(kFALSE);

        fSignal4.SetName("Signal4");
        fSignal4.SetTitle("");
        fSignal4.SetXTitle("Time [ns]");
        fSignal4.SetYTitle("Amplitude [a.u.]");
        fSignal4.GetXaxis()->CenterTitle();
        fSignal4.GetXaxis()->SetTitleSize(0.045);
        fSignal4.GetYaxis()->SetTitleSize(0.045);
        fSignal4.SetDirectory(NULL);
        fSignal4.SetStats(kFALSE);

        fSignal5.SetName("Signal5");
        fSignal5.SetTitle("");
        fSignal5.SetXTitle("Time [ns]");
        fSignal5.SetYTitle("Amplitude [a.u.]");
        fSignal5.GetXaxis()->CenterTitle();
        fSignal5.GetXaxis()->SetTitleSize(0.045);
        fSignal5.GetYaxis()->SetTitleSize(0.045);
        fSignal5.SetDirectory(NULL);
        fSignal5.SetStats(kFALSE);

        fSignal6.SetName("Signal6");
        fSignal6.SetTitle("");
        fSignal6.SetXTitle("Time [ns]");
        fSignal6.SetYTitle("Amplitude [a.u.]");
        fSignal6.GetXaxis()->CenterTitle();
        fSignal6.GetXaxis()->SetTitleSize(0.045);
        fSignal6.GetYaxis()->SetTitleSize(0.045);
        fSignal6.SetDirectory(NULL);
        fSignal6.SetStats(kFALSE);

        fSignal7.SetName("Signal7");
        fSignal7.SetTitle("");
        fSignal7.SetXTitle("Time [ns]");
        fSignal7.SetYTitle("Amplitude [a.u.]");
        fSignal7.GetXaxis()->CenterTitle();
        fSignal7.GetXaxis()->SetTitleSize(0.045);
        fSignal7.GetYaxis()->SetTitleSize(0.045);
        fSignal7.SetDirectory(NULL);
        fSignal7.SetStats(kFALSE);

        fSignal8.SetName("Signal8");
        fSignal8.SetTitle("");
        fSignal8.SetXTitle("Time [ns]");
        fSignal8.SetYTitle("Amplitude [a.u.]");
        fSignal8.GetXaxis()->CenterTitle();
        fSignal8.GetXaxis()->SetTitleSize(0.045);
        fSignal8.GetYaxis()->SetTitleSize(0.045);
        fSignal8.SetDirectory(NULL);
        fSignal8.SetStats(kFALSE);

        fProfAll.SetName("ProfAll");
        fProf1.SetName("Prof1");
        fProf2.SetName("Prof2");
        fProf3.SetName("Prof3");
        fProf4.SetName("Prof4");
        fProf5.SetName("Prof5");
        fProf6.SetName("Prof6");
        fProf7.SetName("Prof7");
        fProf8.SetName("Prof8");

        fProfAll.SetErrorOption("s");
        fProf1.SetErrorOption("s");
        fProf2.SetErrorOption("s");
        fProf3.SetErrorOption("s");
        fProf4.SetErrorOption("s");
        fProf5.SetErrorOption("s");
        fProf6.SetErrorOption("s");
        fProf7.SetErrorOption("s");
        fProf8.SetErrorOption("s");

        fProfAll.GetYaxis()->SetTitleOffset(0.8);
        fProf1.GetYaxis()->SetTitleOffset(0.8);
        fProf2.GetYaxis()->SetTitleOffset(0.8);
        fProf3.GetYaxis()->SetTitleOffset(0.8);
        fProf4.GetYaxis()->SetTitleOffset(0.8);
        fProf5.GetYaxis()->SetTitleOffset(0.8);
        fProf6.GetYaxis()->SetTitleOffset(0.8);
        fProf7.GetYaxis()->SetTitleOffset(0.8);
        fProf8.GetYaxis()->SetTitleOffset(0.8);
    }

    Bool_t SetupFill(const MParList *plist)
    {
        fSingles = (MSingles*)plist->FindObject("MSingles");
        if (!fSingles)
        {
            *fLog << err << "MSingles not found... abort." << endl;
            return kFALSE;
        }

        fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
        if (!fEvent)
        {
            *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl;
            return kFALSE;
        }

        fBadPix = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
        if (!fBadPix)
        {
            *fLog << err << "MBadPixelsCam not found... abort." << endl;
            return kFALSE;
        }

        fPedestal = (MPedestalCam*)plist->FindObject("MPedestalCam");
        if (!fPedestal)
        {
            *fLog << err << "MPedestalCam not found... abort." << endl;
            return kFALSE;
        }

        fNumEntries = 0;

        return kTRUE;
    }

    Bool_t ReInit(MParList *plist)
    {
        const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
        if (!header)
        {
            *fLog << err << "MRawRunHeader not found... abort." << endl;
            return kFALSE;
        }

        // Setup binning
        //const Int_t w = fSingles->GetIntegrationWindow();

        MBinning binsx, binsy;
        binsx.SetEdges(1024,  -512, 512);
        binsy.SetEdges( 800,  -20,   20);

        MH::SetBinning(fSignalAll, binsx, binsy);
        MH::SetBinning(fSignalOut, binsx, binsy);
        MH::SetBinning(fSignal1,   binsx, binsy);
        MH::SetBinning(fSignal2,   binsx, binsy);
        MH::SetBinning(fSignal3,   binsx, binsy);
        MH::SetBinning(fSignal4,   binsx, binsy);
        MH::SetBinning(fSignal5,   binsx, binsy);
        MH::SetBinning(fSignal6,   binsx, binsy);
        MH::SetBinning(fSignal7,   binsx, binsy);
        MH::SetBinning(fSignal8,   binsx, binsy);
        MH::SetBinning(fProfAll,   binsx);
        MH::SetBinning(fProf1,     binsx);
        MH::SetBinning(fProf2,     binsx);
        MH::SetBinning(fProf3,     binsx);
        MH::SetBinning(fProf4,     binsx);
        MH::SetBinning(fProf5,     binsx);
        MH::SetBinning(fProf6,     binsx);
        MH::SetBinning(fProf7,     binsx);
        MH::SetBinning(fProf8,     binsx);

        fSignal1.GetXaxis()->SetRangeUser(-25, 300);
        fSignal1.GetYaxis()->SetRangeUser(-1.5, 2.5);

        fSignal2.GetXaxis()->SetRangeUser(-25, 300);
        fSignal2.GetYaxis()->SetRangeUser(-1.5, 3.5);

        fSignal3.GetXaxis()->SetRangeUser(-25, 300);
        fSignal3.GetYaxis()->SetRangeUser(-1.5, 4.5);

        fSignal4.GetXaxis()->SetRangeUser(-25, 300);
        fSignal4.GetYaxis()->SetRangeUser(-1.5, 5.5);

        fSignal5.GetXaxis()->SetRangeUser(-25, 300);
        fSignal5.GetYaxis()->SetRangeUser(-1.5, 6.5);

        fSignal6.GetXaxis()->SetRangeUser(-25, 300);
        fSignal6.GetYaxis()->SetRangeUser(-1.5, 7.5);

        fSignal7.GetXaxis()->SetRangeUser(-25, 300);
        fSignal7.GetYaxis()->SetRangeUser(-1.5, 8.5);

        fSignal8.GetXaxis()->SetRangeUser(-25, 300);
        fSignal8.GetYaxis()->SetRangeUser(-1.5, 9.5);

        return kTRUE;
    }

    // Fill singles into histograms
    Int_t Fill(const MParContainer *, const Stat_t)
    {
        // Get pointer to samples to fill samples
        const Float_t *ptr = fEvent->GetSamples(0);
        const Short_t  roi = fEvent->GetNumSamples();

        if (!ptr)
            return kTRUE;

        const Int_t nx = fSignalAll.GetNbinsX();
        const Int_t ny = fSignalAll.GetNbinsY();

        const double x0 = fSignalAll.GetXaxis()->GetBinCenter(1);
        const double y0 = fSignalAll.GetYaxis()->GetBinCenter(1);

        const double wx = fSignalAll.GetXaxis()->GetBinCenter(nx)-x0;
        const double wy = fSignalAll.GetYaxis()->GetBinCenter(ny)-y0;

        // Loop over all pixels
        for (unsigned int i=0; i<64; i++, ptr+=roi)
        {
            int idx = i;//pmap.hw(i).index;

            if ((*fBadPix)[i].IsUnsuitable())
                continue;

            if (!hgain->IsUsed(idx))
                continue;

            if (idx!=0)
                continue;

            const double gain  = hgain->GetBinContent(idx+1)/30;
            const double base  = hbase->GetBinContent(idx+1);

            const double ped = (*fPedestal)[i].GetPedestal();

            // loop over all singles
            const vector<Single> &result = fSingles->GetVector(i);
            if (result.size()!=1)
                continue;

            // Fill pulse shape
            const double tm  = result[0].fTime;
            const double sig = result[0].fSignal/30 - base;

            bool ok = true;//tm>3 && (ptr[int(tm)-3]-ped-base)/gain > -5;

            double mean = 0;

            for (int s=20; s<roi-20; s++)
            {
                double x = (s - tm)/2; // [ns]
                double y = (ptr[s]-ped-base)/gain;

                if (x>-100 && x<-10)
                    mean += y/90;

                if (y<-3)
                    ok = false;
                if (x>5 && x<40 && y<=x/40*-0.5)
                    ok = false;
                if (x>-150 && x<-5 && y>1.5)
                    ok = false;
                if (y>2.5)
                    ok = false;
            }

            // Nothing at or bleow zero from 5 to 40machen

            bool ok1 = sig>gain*0.5 && sig<gain*1.5;
            bool ok2 = sig>gain*1.5 && sig<gain*2.5;
            bool ok3 = sig>gain*2.5 && sig<gain*3.5;
            bool ok4 = sig>gain*3.5 && sig<gain*4.5;
            bool ok5 = sig>gain*4.5 && sig<gain*5.5;
            bool ok6 = sig>gain*5.5 && sig<gain*6.5;
            bool ok7 = sig>gain*6.5 && sig<gain*7.5;
            bool ok8 = sig>gain*7.5 && sig<gain*8.5;


            for (int s=20; s<roi-20; s++)
            {
                double x = (s - tm)/2; // [ns]
                double y = (ptr[s]-ped-base)/gain -  mean;

                int ix = TMath::Nint((x-x0)/wx*nx);
                int iy = TMath::Nint((y-y0)/wy*ny);

                if (ix<0 || iy<0 || ix>=nx || iy>=ny)
                    continue;

                int bin = (iy+1)*(nx+2)+(ix+1);

                if (ok)
                {
                    fSignalAll.AddBinContent(bin);
                    fProfAll.Fill(x, y);

                    if (ok1)
                    {
                        fSignal1.AddBinContent(bin);
                        fProf1.Fill(x, y);
                    }
                    if (ok2)
                    {
                        fSignal2.AddBinContent(bin);
                        fProf2.Fill(x, y);
                    }
                    if (ok3)
                    {
                        fSignal3.AddBinContent(bin);
                        fProf3.Fill(x, y);
                    }
                    if (ok4)
                    {
                        fSignal4.AddBinContent(bin);
                        fProf4.Fill(x, y);
                    }
                    if (ok5)
                    {
                        fSignal5.AddBinContent(bin);
                        fProf5.Fill(x, y);
                    }
                    if (ok6)
                    {
                        fSignal6.AddBinContent(bin);
                        fProf6.Fill(x, y);
                    }
                    if (ok7)
                    {
                        fSignal7.AddBinContent(bin);
                        fProf7.Fill(x, y);
                    }
                    if (ok8)
                    {
                        fSignal8.AddBinContent(bin);
                        fProf8.Fill(x, y);
                    }
                }
                else
                    fSignalOut.AddBinContent(bin);
            }

            if (ok)
            {
                fSignalAll.SetEntries(fSignalAll.GetEntries()+1);

                if (ok1)
                    fSignal1.SetEntries(fSignal1.GetEntries()+1);
                if (ok2)
                    fSignal2.SetEntries(fSignal2.GetEntries()+1);
                if (ok3)
                    fSignal3.SetEntries(fSignal3.GetEntries()+1);
                if (ok4)
                    fSignal4.SetEntries(fSignal4.GetEntries()+1);
                if (ok5)
                    fSignal5.SetEntries(fSignal5.GetEntries()+1);
                if (ok6)
                    fSignal6.SetEntries(fSignal6.GetEntries()+1);
                if (ok7)
                    fSignal7.SetEntries(fSignal7.GetEntries()+1);
                if (ok8)
                    fSignal8.SetEntries(fSignal8.GetEntries()+1);
            }
            else
                fSignalOut.SetEntries(fSignalOut.GetEntries()+1);
        }

        fNumEntries++;

        return kTRUE;
    }

    // Getter for histograms
    TH2 *GetSignalAll() { return &fSignalAll; }
    TH2 *GetSignalOut() { return &fSignalOut; }
    TH2 *GetSignal1() { return &fSignal1; }
    TH2 *GetSignal2() { return &fSignal2; }
    TH2 *GetSignal3() { return &fSignal3; }
    TH2 *GetSignal4() { return &fSignal4; }
    TH2 *GetSignal5() { return &fSignal5; }
    TH2 *GetSignal6() { return &fSignal6; }
    TH2 *GetSignal7() { return &fSignal7; }
    TH2 *GetSignal8() { return &fSignal8; }
    TH1 *GetProf1() { return &fProf1; }
    TH1 *GetProf2() { return &fProf2; }
    TH1 *GetProf3() { return &fProf3; }
    TH1 *GetProf4() { return &fProf4; }
    TH1 *GetProf5() { return &fProf5; }
    TH1 *GetProf6() { return &fProf6; }
    TH1 *GetProf7() { return &fProf7; }
    TH1 *GetProf8() { return &fProf8; }

    UInt_t GetNumEntries() const { return fNumEntries; }

    void Draw(Option_t *)
    {
        TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);

        AppendPad("");

        pad->Divide(2,2);

        pad->cd(1);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fSignalAll.Draw("colz");
        fProfAll.Draw("same");

        pad->cd(2);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fSignalOut.Draw("colz");

        pad->cd(3);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fSignal1.Draw("colz");
        fProf1.Draw("same");

        pad->cd(4);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fSignal8.Draw("colz");
        //fProf8.Draw("same");
    }

    void DrawCopy(Option_t *)
    {
        TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);

        AppendPad("");

        pad->Divide(2,2);

        pad->cd(1);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fSignalAll.DrawCopy("colz");

        pad->cd(2);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fSignalOut.DrawCopy("colz");

        pad->cd(3);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fSignal1.DrawCopy("colz");

        pad->cd(4);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fSignal8.DrawCopy("colz");
    }

    ClassDef(MHSingles, 2)
};

// Task to extract the singles
class MExtractSingles : public MTask
{
    MSingles *fSingles;
    MPedestalCam *fPedestal;
    MPedestalSubtractedEvt *fEvent;

    // On time for each pixel in samples
    MArrayI fExtractionRange;

    // Number of samples for sliding average
    size_t fNumAverage;

    // The range in which the singles have to fit in
    // FIXME: Make sure this is consistent with MHBaseline
    UInt_t fStartSkip;
    UInt_t fEndSkip;

    UInt_t fIntegrationSize;
    UInt_t fMaxSearchWindow;

    Int_t    fMaxDist;
    Double_t fThreshold;

public:
    MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0),
        fExtractionRange(64),
        fNumAverage(10), fStartSkip(20), fEndSkip(10),
        fIntegrationSize(3*10), fMaxSearchWindow(30)
    {
    }

    void SetMaxDist(Int_t m) { fMaxDist = m; }
    void SetThreshold(Int_t th) { fThreshold = th; }

    UInt_t GetIntegrationSize() const { return fIntegrationSize; }

    const MArrayI &GetExtractionRange() const { return fExtractionRange; }


    Int_t PreProcess(MParList *plist)
    {
        fSingles = (MSingles*)plist->FindCreateObj("MSingles");
        if (!fSingles)
            return kFALSE;

        fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
        if (!fEvent)
        {
            *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl;
            return kFALSE;
        }

        fPedestal = (MPedestalCam*)plist->FindObject("MPedestalCam");
        if (!fPedestal)
        {
            *fLog << err << "MPedestalCam not found... abort." << endl;
            return kFALSE;
        }

        return kTRUE;
    }

    Int_t Process()
    {
        const UInt_t roi = fEvent->GetNumSamples();

        const size_t navg              = fNumAverage;
        const float  threshold         = fThreshold; 
        const UInt_t start_skip        = fStartSkip;
        const UInt_t end_skip          = fEndSkip;
        const UInt_t integration_size  = fIntegrationSize;
        const UInt_t max_search_window = fMaxSearchWindow;
        const UInt_t max_dist          = fMaxDist;

        vector<float> val(roi-navg);//result of first sliding average
        for (int pix=0; pix<64; pix++)
        {
            // Total number of samples as 'on time'
            fExtractionRange[pix] += roi-start_skip-navg-end_skip-integration_size;

            // Clear singles for this pixel
            vector<Single> &result = fSingles->GetVector(pix);
            result.clear();

            // Get pointer to samples
            const Float_t *ptr = fEvent->GetSamples(pix);

            // Get Baseline
            const Double_t ped = (*fPedestal)[pix].GetPedestal();

            // Subtract baseline and do a sliding average over
            // the samples to remove noise (mainly the 200MHz noise)
            for (size_t i=0; i<roi-navg; i++)
            {
                val[i] = 0;
                for (size_t j=i; j<i+navg; j++)
                    val[i] += ptr[j];
                val[i] /= navg;
                val[i] -= ped;
            }

            // peak finding via threshold
            UInt_t imax = roi-navg-end_skip-integration_size;
            for (UInt_t i=start_skip; i<imax; i++)
            {
                //search for threshold crossings
                if (val[i+0]>threshold ||
                    val[i+4]>threshold ||

                    val[i+5]<threshold ||
                    val[i+9]<threshold)
                    continue;

                //search for maximum after threshold crossing
                UInt_t k_max = i+5;
                for (UInt_t k=i+6; k<i+max_search_window; k++)
                {
                    if (val[k] > val[k_max])
                        k_max = k;
                }

                // Check if the findings make sense
                if (k_max == i+5 || k_max == i + max_search_window-1)
                    continue;

                //search for half maximum before maximum
                UInt_t k_half_max = 0;
                for (UInt_t k=k_max; k>k_max-25; k--)
                {
                    if (val[k-1] < val[k_max]/2 &&
                        val[k] >= val[k_max]/2 )
                    {
                        k_half_max = k;
                        break;
                    }
                }
                // Check if the finding makes sense
                if (k_half_max+integration_size > roi-navg-end_skip)
                    continue;
                if (k_half_max == 0)
                    continue;
                if (k_max - k_half_max > max_dist)
                    continue;

                // FIXME: Evaluate arrival time more precisely!!!
                // FIXME: Get a better integration window

                // Compile "single"
                Single single;
                single.fSignal = 0;
                single.fTime   = k_half_max;

                // Crossing of the threshold is at 5
                for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
                    single.fSignal += ptr[j];

                single.fSignal -= integration_size*ped;

                // Add single to list
                result.push_back(single);

                // skip falling edge
                i += 5+integration_size;

                // Remove skipped samples from 'on time'
                fExtractionRange[pix] -= i>imax ? 5+integration_size-(i-imax) : 5+integration_size;
            }
        }
        return kTRUE;
    }
};

int extract_pulse(
                    Int_t        max_dist   = 14,
                    Double_t     threshold  =  5,
                    const char  *runfile    = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_",
                    int          firstRunID =  6,
                    int          lastRunID  =  6,
                    const char  *drsfile    = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_003.drs.fits",
                    const char  *outpath    = "."
                   )
{
    // ======================================================
    if (!pmap.Read("../hawc/FAMOUSmap171218.txt"))
    {
        gLog << err << "../hawc/FAMOUSmap171218.txt not found." << endl;
        return 1;
    }

    // ======================================================


    MDirIter iter;

    // built output filename and path
    TString outfile = "20191027_006_006-pulse.root";
    TString infile  = "20190227_006_006-fit.root";

    //if (!gSystem->AccessPathName(outfile))
    //    return 0;

    // add input files to directory iterator
    for (Int_t fileNr = firstRunID; fileNr <= lastRunID; fileNr++)
    {
        TString filename = runfile;
        filename += Form("%03d.fits", fileNr);
        iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".fz"));
        iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".gz"));
        iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename));
    }

    // ======================================================

    TFile file(infile);
    MStatusArray arr;
    if (arr.Read()<=0)
    {
        gLog << "ERROR - Cannot open gain file '" << infile << "'" << endl;
        return 1;
    }

    hgain = (MHCamera*)arr.FindObjectInCanvas("Gain_copy",     "MHCamera", "Cams1");
    hbase = (MHCamera*)arr.FindObjectInCanvas("Baseline_copy", "MHCamera", "Cams1");
    if (!hgain || !hbase)
    {
        gLog << "ERROR - Cannot find Gain/Baseline" << endl;
        return 1;
    }

    // ======================================================

    // true:  Display correctly mapped pixels in the camera displays
    //        but the value-vs-index plot is in software/spiral indices
    // false: Display pixels in hardware/linear indices,
    //        but the order is the camera display is distorted
    bool usemap = true;

    // map file to use (get that from La Palma!)
    const char *map = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;

    // ------------------------------------------------------

    Long_t max0 = 0;
    Long_t max1 = 0;

    // ======================================================

    if (map && gSystem->AccessPathName(map, kFileExists))
    {
        gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
        return 1;
    }
    if (gSystem->AccessPathName(drsfile, kFileExists))
    {
        gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
        return 2;
    }

    // --------------------------------------------------------------------------------

    gLog.Separator("Pulses");
    gLog << all << "Extract pulses with '" << drsfile << "'" << endl;
    gLog << endl;

    // ------------------------------------------------------

    MDrsCalibration drscalib300;
    if (!drscalib300.ReadFits(drsfile))
        return 3;

    // ------------------------------------------------------

    iter.Sort();
    iter.Print();

    TString title;
    title =  iter.Next();
    title += "; ";
    title += max1;
    iter.Reset();

    // ======================================================

    MStatusDisplay *d = new MStatusDisplay;

    MBadPixelsCam badpixels;
    badpixels.InitSize(64);
    badpixels[7].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
    badpixels[17].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);

    MH::SetPalette();

    MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");

    MGeomApply apply;

    MDrsCalibApply drsapply;
    drsapply.SetMaxNumPrevEvents(0); // Switched off step-correction -> crashes due to the requirement of N*9 ch

    MRawFitsRead read;
    read.LoadMap(map);
    read.AddFiles(iter);

    // ======================================================

    gLog << endl;
    gLog.Separator("Calculating baseline");

    MTaskList tlist0;

    MRawRunHeader header;
    MPedestalCam  pedcam;

    MParList plist;
    plist.AddToList(&tlist0);
    plist.AddToList(&drscalib300);
    plist.AddToList(&header);
    plist.AddToList(&pedcam);

    // ------------------ Setup the tasks ---------------

    MFillH fill0("MHBaseline", "MPedestalSubtractedEvt");

    drsapply.SetRemoveSpikes(4);

    // ------------------ Setup eventloop and run analysis ---------------

    tlist0.AddToList(&read);
    tlist0.AddToList(&apply);
    tlist0.AddToList(&drsapply);
    tlist0.AddToList(&fill0);

    // ------------------ Setup and run eventloop ---------------

    MEvtLoop loop0(title);
    loop0.SetDisplay(d);
    loop0.SetParList(&plist);

    if (!loop0.Eventloop(max0))
        return 4;

    if (!loop0.GetDisplay())
        return 5;

    // ----------------------------------------------------------------

    MGeomCamFAMOUS fact;
    MHCamera hped(fact);
    hped.SetCamContent(pedcam);
    hped.SetCamError(pedcam, 4);
    hped.SetAllUsed();
    MHCamEvent display;
    display.SetHist(hped);

    d->AddTab("Baseline");
    display.Clone()->Draw();

    // ======================================================

    gLog << endl;
    gLog.Separator("Extracting singles");

    MTaskList tlist1;

    MSingles singles;

    plist.Replace(&tlist1);
    plist.AddToList(&badpixels);
    plist.AddToList(&singles);

    // ------------------ Setup the tasks ---------------

    MExtractSingles extract;
    extract.SetMaxDist(max_dist);
    extract.SetThreshold(threshold);

    MFillH fill("MHSingles");

    // ------------------ Setup eventloop and run analysis ---------------

    tlist1.AddToList(&read);
    tlist1.AddToList(&apply);
    tlist1.AddToList(&drsapply);
    tlist1.AddToList(&extract);
    tlist1.AddToList(&fill);

    // ------------------ Setup and run eventloop ---------------

    MEvtLoop loop1(title);
    loop1.SetDisplay(d);
    loop1.SetParList(&plist);

    if (!loop1.Eventloop(max1))
        return 6;

    if (!loop1.GetDisplay())
        return 7;

    if (fill.GetNumExecutions()==0)
        return 8;

    // =============================================================

    MHSingles* h = (MHSingles*) plist.FindObject("MHSingles");
    if (!h)
        return 9;

    gStyle->SetOptStat(10);

    TH1 *h1 = h->GetSignal1();
    TH1 *p1 = h->GetProf1();

    TH1 *h2 = h->GetSignal2();
    TH1 *p2 = h->GetProf2();

    TH1 *h3 = h->GetSignal3();
    TH1 *p3 = h->GetProf3();

    TH1 *h4 = h->GetSignal4();
    TH1 *p4 = h->GetProf4();

    TH1 *h5 = h->GetSignal5();
    TH1 *p5 = h->GetProf5();

    TH1 *h6 = h->GetSignal6();
    TH1 *p6 = h->GetProf6();

    TH1 *h7 = h->GetSignal7();
    TH1 *p7 = h->GetProf7();

    TH1 *h8 = h->GetSignal8();
    TH1 *p8 = h->GetProf8();

    TF1 pulse("pulse", "[0]*(1-exp(-x*x/[1]))*exp(-x/[2])", 0, 300);
    pulse.SetLineColor(kBlack);
    pulse.SetParameters(0, 1./0.094, 1./0.05);

    // =====

    h1->SetMinimum(h1->GetEntries()/1000);
    h2->SetMinimum(h2->GetEntries()/1000);
    h3->SetMinimum(h3->GetEntries()/1000);
    h4->SetMinimum(h4->GetEntries()/1000);
    h5->SetMinimum(h5->GetEntries()/1000);
    h6->SetMinimum(h6->GetEntries()/1000);
    h7->SetMinimum(h7->GetEntries()/1000);
    h8->SetMinimum(h8->GetEntries()/1000);

    // =====

    d->AddTab("1");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    h1->DrawCopy("colz");
    p1->DrawCopy("same");

    pulse.SetParameter(0, 1*1.3);
    p1->Fit(&pulse, "N0", "", 1.5, 25);
    TF1 *f1 = pulse.DrawCopy("same");

    // =====

    d->AddTab("2");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    h2->DrawCopy("colz");
    p2->DrawCopy("same");

    pulse.SetParameter(0, 2*1.3);
    p2->Fit(&pulse, "N0", "", 1.5, 25);
    TF1 *f2 = pulse.DrawCopy("same");

    // =====

    d->AddTab("3");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    h3->DrawCopy("colz");
    p3->DrawCopy("same");

    pulse.SetParameter(0, 3*1.3);
    p3->Fit(&pulse, "N0", "", 1.5, 25);
    TF1 *f3 = pulse.DrawCopy("same");

    // =====

    d->AddTab("4");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    h4->DrawCopy("colz");
    p4->DrawCopy("same");

    pulse.SetParameter(0, 4*1.3);
    p4->Fit(&pulse, "N0", "", 1.5, 25);
    TF1 *f4 = pulse.DrawCopy("same");

    // =====

    d->AddTab("5");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    h5->DrawCopy("colz");
    p5->DrawCopy("same");

    pulse.SetParameter(0, 5*1.3);
    p5->Fit(&pulse, "N0", "", 1.5, 25);
    TF1 *f5 = pulse.DrawCopy("same");

    // =====

    d->AddTab("6");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    h6->DrawCopy("colz");
    p6->DrawCopy("same");

    pulse.SetParameter(0, 6*1.3);
    p6->Fit(&pulse, "N0", "", 1.5, 25);
    TF1 *f6 = pulse.DrawCopy("same");

    // =====

    d->AddTab("7");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    h7->DrawCopy("colz");
    p7->DrawCopy("same");

    pulse.SetParameter(0, 7*1.3);
    p7->Fit(&pulse, "N0", "", 1.5, 25);
    TF1 *f7 = pulse.DrawCopy("same");

    // =====

    d->AddTab("8");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    h8->DrawCopy("colz");
    p8->DrawCopy("same");

    pulse.SetParameter(0, 6*1.3);
    p8->Fit(&pulse, "N0", "", 1.5, 25);
    TF1 *f8 = pulse.DrawCopy("same");

    // =====

    d->AddTab("Pulse");
    gPad->SetGrid();
    gPad->SetLeftMargin(0.07);
    gPad->SetTopMargin(0.01);
    gPad->SetRightMargin(0.15);
    p8->SetStats(kFALSE);
    p8->GetXaxis()->SetRangeUser(-25, 300);
    p8->GetXaxis()->CenterTitle();
    p8->GetXaxis()->SetTitleSize(0.045);
    p8->GetYaxis()->SetTitleSize(0.045);
    p8->GetYaxis()->SetTitleOffset(0.8);
    p8->SetYTitle("Amplitude");
    p8->SetXTitle("Time [ns]");
    p8->DrawCopy();
    p7->DrawCopy("same");
    p6->DrawCopy("same");
    p5->DrawCopy("same");
    p4->DrawCopy("same");
    p3->DrawCopy("same");
    p2->DrawCopy("same");
    p1->DrawCopy("same");

    f1->Draw("same");
    f2->Draw("same");
    f3->Draw("same");
    f4->Draw("same");
    f5->Draw("same");
    f6->Draw("same");
    f7->Draw("same");
    f8->Draw("same");

    d->SaveAs("paper_pulse.root");


    /*
    d->SaveAs(outfile);

    TFile f(outfile, "UPDATE");

    MParameterI par("NumEvents");
    par.SetVal(fill.GetNumExecutions());
    par.Write();

    MParameterI win("IntegrationWindow");
    win.SetVal(extract.GetIntegrationSize());
    win.Write();

    extract.GetExtractionRange().Write("ExtractionRange");

    if (firstRunID==lastRunID)
        header.Write();
    */
    return 0;
}