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

#include "MH.h"
#include "MLog.h"
#include "MLogManip.h"
#include "MDirIter.h"
#include "MFillH.h"
#include "MEvtLoop.h"
#include "MCamEvent.h"
#include "MGeomApply.h"
#include "MTaskList.h"
#include "MParList.h"
#include "MContinue.h"
#include "MBinning.h"
#include "MHDrsCalib.h"
#include "MDrsCalibApply.h"
#include "MRawFitsRead.h"
#include "MBadPixelsCam.h"
#include "MStatusDisplay.h"
#include "MTaskInteractive.h"
#include "MPedestalSubtractedEvt.h"
#include "MHCamera.h"
#include "MGeomCamFACT.h"
#include "MRawRunHeader.h"

using namespace std;

MPedestalSubtractedEvt *fEvent = 0;
MLog *fLog = &gLog;

struct Single
{
    float    fSignal;
    float    fTime;
};

class MSingles : public MParContainer, public MCamEvent
{
    Int_t fIntegrationWindow;

    vector<vector<Single>> fData;

public:
    MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0)
    {
        fName = name ? name : "MSingles";
    }

    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 &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
    {
        return kTRUE;
    }
    void   DrawPixelContent(Int_t num) const { }

    ClassDef(MSingles, 1)
};

class MH2F : public TH2F
{
public:
    Int_t Fill(Double_t x,Double_t y)
    {
        Int_t binx, biny, bin;
        fEntries++;
        binx = TMath::Nint(x+1);
        biny = fYaxis.FindBin(y);
        bin  = biny*(fXaxis.GetNbins()+2) + binx;
        AddBinContent(bin);
        if (fSumw2.fN) ++fSumw2.fArray[bin];
        if (biny == 0 || biny > fYaxis.GetNbins()) {
            if (!fgStatOverflows) return -1;
        }
        ++fTsumw;
        ++fTsumw2;
        fTsumwy  += y;
        fTsumwy2 += y*y;
        return bin;
    }
    ClassDef(MH2F, 1);
};

class MProfile2D : public TProfile2D
{
public:
    Int_t Fill(Double_t x, Double_t y, Double_t z)
    {
        Int_t bin,binx,biny;
        fEntries++;
        binx =TMath::Nint(x+1);
        biny =fYaxis.FindBin(y);
        bin = GetBin(binx, biny);
        fArray[bin] += z;
        fSumw2.fArray[bin] += z*z;
        fBinEntries.fArray[bin] += 1;
        if (fBinSumw2.fN)  fBinSumw2.fArray[bin] += 1;
        if (biny == 0 || biny > fYaxis.GetNbins()) {
            if (!fgStatOverflows) return -1;
        }
        ++fTsumw;
        ++fTsumw2;
        fTsumwy  += y;
        fTsumwy2 += y*y;
        fTsumwz  += z;
        fTsumwz2 += z*z;
        return bin;
    }
    ClassDef(MProfile2D, 1);
};



class MHSingles : public MH
{
    MH2F       fSignal;
    MH2F       fTime;
    MProfile2D fPulse;

    UInt_t fNumEntries;

    MSingles               *fSingles;   //!
    MPedestalSubtractedEvt *fEvent;     //!

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

        fSignal.SetName("Signal");
        fSignal.SetTitle("Title");
        fSignal.SetXTitle("X title");
        fSignal.SetYTitle("Y title");
        fSignal.SetDirectory(NULL);

        fTime.SetName("Time");
        fTime.SetTitle("Title");
        fTime.SetXTitle("X title");
        fTime.SetYTitle("Y title");
        fTime.SetDirectory(NULL);

        fPulse.SetName("Pulse");
        fPulse.SetTitle("Title");
        fPulse.SetXTitle("X title");
        fPulse.SetYTitle("Y title");
        fPulse.SetDirectory(NULL);
    }

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

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

        const Int_t w = fSingles->GetIntegrationWindow();

        MBinning binsx, binsy;
        binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
        binsy.SetEdges(2*150, -2*7.5*w, 2*(60-7.5)*w);

        MH::SetBinning(fSignal, binsx, binsy);

        const UShort_t roi = header->GetNumSamples();

        // Correct for DRS timing!!
        MBinning binst(roi, -0.5, roi-0.5);
        MH::SetBinning(fTime, binsx, binst);

        MBinning binspy(2*roi, -roi-0.5, roi-0.5);
        MH::SetBinning(fPulse, binsx, binspy);

        return kTRUE;
    }

    Int_t Fill(const MParContainer *par, const Stat_t weight=1)
    {
        const Float_t *ptr = fEvent->GetSamples(0);
        const Short_t  roi = fEvent->GetNumSamples();

        for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
        {
            const vector<Single> &result = fSingles->GetVector(i);

            for (unsigned int j=0; j<result.size(); j++)
            {
                fSignal.Fill(i, result[j].fSignal);
                fTime.Fill  (i, result[j].fTime);

                if (!ptr)
                    continue;

                const Short_t tm = floor(result[j].fTime);

                for (int s=0; s<roi; s++)
                    fPulse.Fill(i, s-tm, ptr[s]);
            }

            ptr += roi;
        }

        fNumEntries++;

        return kTRUE;
    }

    TH2 *GetSignal() { return &fSignal; }
    TH2 *GetTime()   { return &fTime; }
    TH2 *GetPulse()  { return &fPulse; }

    UInt_t GetNumEntries() const { return fNumEntries; }

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

        AppendPad("");

        pad->Divide(2,2);

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

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

        pad->cd(3);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        fPulse.Draw("colz");
    }

    ClassDef(MHSingles, 1)
};

MStatusDisplay *d = new MStatusDisplay;

MSingles *fSingles;

TH1F *fluct1 = new TH1F("", "", 200, -200, 200);
TH1F *fluct2 = new TH1F("", "", 100, -50, 50);

TGraph weights("weights.txt");

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

    d->AddTab("Fluct");
    fluct2->Draw();
    fluct1->Draw("same");

    fSingles->SetIntegrationWindow(20);

    //if (weights.GetN()==0)
    //    return kFALSE;

    return kTRUE;
}

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

    const size_t navg             = 10;
    const float  threshold        = 5;
    const UInt_t start_skip       = 20;
    const UInt_t end_skip         = 10;
    const UInt_t integration_size = 10*3;
    const UInt_t max_search_window = 30;

    vector<float> val(roi-navg);//result of first sliding average
    for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
    {
        vector<Single> &result = fSingles->GetVector(pix);
        result.clear();

        const Float_t *ptr = fEvent->GetSamples(pix);

        //Sliding window
        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;
        }

        //peak finding via threshold
        for (UInt_t i=start_skip; i<roi-navg-end_skip-30; 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+5; k<i+max_search_window; k++)
            {
                if (val[k] > val[k_max])
                    k_max = k;
            }

            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;
                }
            }
            if (k_half_max > roi-navg-end_skip-max_search_window - integration_size)
                continue;
            if (k_half_max == 0)
                continue;
            if (k_max - k_half_max > 14)
                continue;

            fluct2->Fill(k_max - k_half_max);

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

            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];
            result.push_back(single);

            i += 5+integration_size;
        }
    }
    return kTRUE;
}

Int_t PostProcess()
{
    return kTRUE;
}

Double_t fcn(Double_t *x, Double_t *par);

int singlepe()
{
    // ======================================================

    const char *drsfile = "/fact/raw/2012/06/01/20120601_003.drs.fits.gz";

    MDirIter iter;
    iter.AddDirectory("/fact/raw/2012/06/01", "20120601_005.fits.gz");

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

    // 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 ? "/scratch_nfs/FACTmap111030.txt" : NULL;

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

    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("Singles");
    gLog << "Extract singles with '" << drsfile << "'" << endl;
    gLog << endl;

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

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

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

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

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

    //MStatusDisplay *d = new MStatusDisplay;

    MBadPixelsCam badpixels;
    badpixels.InitSize(1440);
    badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);

    //  Twin pixel
    //     113
    //     115
    //     354
    //     423
    //    1195
    //    1393

    //MDrsCalibrationTime timecam;

    MH::SetPalette();

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

    gLog << endl;
    gLog.Separator("Processing external light pulser run");

    MTaskList tlist1;

    MSingles singles;
    MRawRunHeader header;

    MParList plist1;
    plist1.AddToList(&tlist1);
    plist1.AddToList(&drscalib300);
    plist1.AddToList(&badpixels);
    plist1.AddToList(&singles);
    plist1.AddToList(&header);

    MEvtLoop loop1("DetermineCalConst");
    loop1.SetDisplay(d);
    loop1.SetParList(&plist1);

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

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

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

    MGeomApply apply1;

    MDrsCalibApply drsapply1;

    MTaskInteractive mytask;
    mytask.SetPreProcess(PreProcess);
    mytask.SetProcess(Process);
    mytask.SetPostProcess(PostProcess);

    MFillH fill("MHSingles");

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

    tlist1.AddToList(&read1);
//    tlist1.AddToList(&cont1);
    tlist1.AddToList(&apply1);
    tlist1.AddToList(&drsapply1);
    tlist1.AddToList(&mytask);
    tlist1.AddToList(&fill);

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

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

    gStyle->SetOptFit(kTRUE);

    MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
    if (!h)
        return 0;

    TH2 *hsignal = h->GetSignal();
    TH2 *htime   = h->GetTime();
    TH2 *hpulse  = h->GetPulse();

    d->AddTab("Time");
    TH1 *ht = htime->ProjectionY();
    ht->Scale(1./1440);
    ht->Draw();

    d->AddTab("Pulse");
    TH1 *hp = hpulse->ProjectionY();
    hp->Scale(1./1440);
    hp->Draw();


    // ------------------ Spectrum Fit ---------------
    // AFTER this the Code for the spektrum fit follows
    // access the spektrumhistogram by the pointer *hist

    MGeomCamFACT fact;
    MHCamera hcam(fact);

    //built an array of Amplitude spectra
    TH1F hAmplitude("Rate",      "Average number of dark counts per event",
                    100,  0,  10);
    TH1F hGain     ("Gain",      "Gain distribution",
                     50,  200,   350);
    TH1F hGainRMS  ("RelSigma",  "Distribution of rel. Sigma",
                    100,   0,   30);
    TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability",
                    100,   0,    25);
    TH1F hOffset   ("Offset",    "Baseline offset distribution",
                    50, -7.5,  2.5);
    TH1F hFitProb  ("FitProb",   "FitProb distribution",
                     50,   0,  100);
    TH1F hSum1     ("Sum1",      "Sum",
                    100,    0,  2000);
    TH1F hSum2     ("Sum2",      "Sum",
                    100,    0,   10);

    // define fit function for Amplitudespectrum
    TF1 func("spektrum", fcn, 0, 1000, 5);
    func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
    func.SetLineColor(kBlue);

    // define fit function for Amplitudespectrum
    TF1 func2("sum_spektrum", fcn, 0, 2000, 5);
    func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
    func2.SetLineColor(kBlue);

    // define fit function for Amplitudespectrum
    TF1 func3("gain_sum_spektrum", fcn, 0, 10, 5);
    func3.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
    func3.SetLineColor(kBlue);



    // Map for which pixel shall be plotted and which not
    TArrayC usePixel(1440);
    memset(usePixel.GetArray(), 1, 1440);

    // List of Pixel that should be ignored in camera view
    usePixel[424]  = 0;
    usePixel[923]  = 0;
    usePixel[1208] = 0;
    usePixel[583]  = 0;
    usePixel[830]  = 0;
    usePixel[1399] = 0;
    usePixel[113]  = 0;
    usePixel[115]  = 0;
    usePixel[354]  = 0;
    usePixel[423]  = 0;
    usePixel[1195] = 0;
    usePixel[1393] = 0;

    int count_ok = 0;

    // Begin of Loop over Pixels
    for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
    {
        if (usePixel[pixel]==0)
            continue;

        //Projectipon of a certain Pixel out of Ampl.Spectrum
        TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
        hist->SetDirectory(0);

        const Int_t    maxbin   = hist->GetMaximumBin();
        const Double_t binwidth = hist->GetBinWidth(maxbin);
        const Double_t ampl     = hist->GetBinContent(maxbin);

        double fwhm = 0;
        for (int i=1; i<maxbin; i++)
            if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2)
            {
                fwhm = (2*i+1)*binwidth;
                break;
            }

        if (fwhm==0)
        {
            cout << "Could not determine start value for sigma." << endl;
            continue;
        }

        // Amplitude, Gain, Sigma, Crosstalk Prob, Shift
        // par[1] + par[4] is the position of the first maximum
        Double_t par[5] =
        {
            ampl, hist->GetBinCenter(maxbin), fwhm*1.4, 0.1, -10
        };

        func.SetParameters(par);

        const double min = par[1]-fwhm*1.4;
        const double max = par[1]*5;

        //Fit Pixels spectrum
        const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", min, max);

        const bool ok = int(rc)==0;

        const double fit_prob   = rc->Prob();
        const float  fGain      = func.GetParameter(1);
        const float  fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0);
        const float  f2GainRMS  = func.GetParameter(2);
        const float  fCrosstlk  = func.GetParameter(3);
        const float  fOffset    = func.GetParameter(4)/singles.GetIntegrationWindow();

        hAmplitude.Fill(fAmplitude);
        hGain.Fill(fGain);
        //hGainRMS.Fill(f2GainRMS);
        hCrosstlk.Fill(fCrosstlk*100);
        hOffset.Fill(fOffset);
        hFitProb.Fill(100*fit_prob);
        hGainRMS.Fill(100*f2GainRMS/fGain);

        hcam.SetBinContent(pixel+1, fGain);

        usePixel[pixel] = ok;

        if (!ok)
        {
            cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
            continue;
        }

        for (int b=1; b<=hist->GetNbinsX(); b++)
        {
            hSum1.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
            hSum2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
        }

        if (count_ok==0)
        {
            TCanvas &c = d->AddTab(Form("Pix%d", pixel));
            c.cd();
            hist->SetName(Form("Pix%d", pixel));
            hist->DrawCopy()->SetDirectory(0);
            func.DrawCopy("SAME");
        }



        count_ok++;

        delete hist;
    }

    //------------------fill histograms-----------------------

    hcam.SetUsed(usePixel);

    TCanvas &canv = d->AddTab("Gain");
    canv.cd();
    hcam.SetNameTitle( "gain","Gain distribution over whole camera");
    hcam.DrawCopy();

    gROOT->SetSelectedPad(0);
    TCanvas &c11 = d->AddTab("SumHist");
    gPad->SetLogy();
    hSum1.DrawCopy();

    //------------------------fit sum spectrum-------------------------------
    const Int_t    maxbin   = hSum1.GetMaximumBin();
    const Double_t binwidth = hSum1.GetBinWidth(maxbin);
    const Double_t ampl     = hSum1.GetBinContent(maxbin);

    double fwhm2 = 0;

    //Calculate full width half Maximum
    for (int i=1; i<maxbin; i++)
        if (hSum1.GetBinContent(maxbin-i)+hSum1.GetBinContent(maxbin+i)<ampl*2)
        {
            fwhm2 = (2*i+1)*binwidth;
            break;
        }
    if (fwhm2==0)
    {
        cout << "Could not determine start value for sigma." << endl;
    }

    //Set fit parameters
    Double_t par[5] =
    {
        ampl, hSum1.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10
    };
    func2.SetParameters(par);

    //calculate fit range
    const double min = par[1]-fwhm2*1.4;
    const double max = par[1]*5;

    //Fit and draw spectrum
    hSum1.Fit(&func2, "NOQS", "", min, max);
    func2.DrawCopy("SAME");

    //-------------------END OF: fit sum spectrum-----------------------------

    // Draw gain corrected sum specetrum
    gROOT->SetSelectedPad(0);
    TCanvas &c12 = d->AddTab("GainHist");
    gPad->SetLogy();
    hSum2.DrawCopy();
    const Double_t fMaxPos  = hSum2.GetBinCenter(maxbin);

    //--------fit gausses to peaks of gain corrected sum specetrum -----------
    UInt_t nfunc    = 5;

    // create array for gaus funtions to fit to peaks of spectrum
    TF1 **gaus      = new TF1*[nfunc];

    // create histo for height of Maximum amplitudes vs. pe
    TH1D hMaxHeights("MaxHeights",      "peak heights",
                      20,  0,   20);

    // fit gauss functions to spectrum
    for (UInt_t nr = 0; nr<nfunc; nr++)
    {
        char fname[20];
        sprintf(fname,"gaus%d",nr);

        //create gauss functions
        gaus[nr] = new TF1( fname,"gaus",
                            fMaxPos*(nr+1)-0.6,     fMaxPos*(nr+1)+0.6);
        //fit and draw gauss functions
        hSum2.Fit(  gaus[nr],   "N0QS",     "",
                    fMaxPos*(nr+1)-0.6,     fMaxPos*(nr+1)+0.6);
        gaus[nr]->DrawCopy("SAME");

        //fill histo for height of Maximum amplitudes vs. pe
        hMaxHeights.SetBinContent(nr, gaus[nr]->GetParameter(0));
        delete gaus[nr];
    }
    delete[] gaus;

    //------------------fit gain corrected sum spectrum-----------------------

    const Int_t    maxbin2   = hSum2.GetMaximumBin();
    const Double_t binwidth2 = hSum2.GetBinWidth(maxbin);
    const Double_t ampl2     = hSum2.GetBinContent(maxbin);

    //Calculate full width half Maximum
    fwhm2 = 0;
    for (int i=1; i<maxbin; i++)
        if (hSum2.GetBinContent(maxbin2-i)+hSum2.GetBinContent(maxbin2+i)<ampl2*2)
        {
            fwhm2 = (2*i+1)*binwidth2;
            break;
        }
    if (fwhm2==0)
    {
        cout << "Could not determine start value for sigma." << endl;
    }

    //Set fit parameters
    Double_t par2[5] =
    {
        ampl2, hSum2.GetBinCenter(maxbin2), fwhm2*1.4, 0.1, 0
    };

    func3.SetParameters(par2);

    const double min2 = par2[1]-fwhm2*1.4;
    const double max2 = par2[1]*5;

    hSum2.Fit(&func3, "N0QS", "", min2, max2);
    func3.DrawCopy("SAME");
    //-------------------END OF: fit sum spectrum-----------------------------

    TCanvas &c15 = d->AddTab("PeakHeights");
    gPad->SetLogy();
    hMaxHeights.DrawCopy();

    TCanvas &c2 = d->AddTab("Result");
    c2.Divide(3,2);
    c2.cd(1);
    gPad->SetLogy();
    hAmplitude.DrawCopy();

    c2.cd(2);
    gPad->SetLogy();
    hGain.DrawCopy();

    c2.cd(3);
    gPad->SetLogy();
    hGainRMS.DrawCopy();

    c2.cd(4);
    gPad->SetLogy();
    hCrosstlk.DrawCopy();

    c2.cd(5);
    gPad->SetLogy();
    hOffset.DrawCopy();

    c2.cd(6);
    gPad->SetLogy();
    hFitProb.DrawCopy();

    /*
    TCanvas &c3 = d->AddTab("Separation");
    gPad->SetLogy();
    hSep.DrawCopy();
     */
    return 0;
}

Double_t fcn(Double_t *xx, Double_t *par)
{
    Double_t x     = xx[0];

    Double_t ampl  = par[0];
    Double_t gain  = par[1];
    Double_t sigma = par[2];
    Double_t cross = par[3];
    Double_t shift = par[4];

    Double_t xTalk = 1;

    Double_t y = 0;
    for (int order = 1; order < 14; order++)
    {
        Double_t arg   = (x - order*gain - shift)/sigma;

        Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);

        y += xTalk*gauss;

        xTalk *= cross;
    }

    return ampl*y;
}