#include <TStyle.h>
#include <TCanvas.h>
#include <TSystem.h>
#include <TF1.h>
#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TMath.h>
#include <TFitResultPtr.h>
#include <TFitResult.h>
#include <TGProgressBar.h>        // TGHProgressBar
#include <TVirtualFitter.h>
#include <Getline.h>

#include <cstdio>
#include <stdio.h>
#include <stdint.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 "MStatusArray.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;
using namespace TMath;

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
{
    TH2F       fSignal;
    TH2F       fTime;
    TProfile2D fPulse;

    UInt_t fNumEntries;

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

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

        fSignal.SetName("Signal");
        fSignal.SetTitle("pulse integral spectrum");
        fSignal.SetXTitle("pixel [SoftID]");
        fSignal.SetYTitle("time [a.u.]");
        fSignal.SetDirectory(NULL);

        fTime.SetName("Time");
        fTime.SetTitle("arival time spectrum");
        fTime.SetXTitle("pixel [SoftID]");
        fTime.SetYTitle("time [a.u.]");
        fTime.SetDirectory(NULL);

        fPulse.SetName("Pulse");
        fPulse.SetTitle("average pulse shape spectrum");
        fPulse.SetXTitle("pixel [SoftID]");
        fPulse.SetYTitle("time [a.u.]");
        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");
    }

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

        AppendPad("");

        pad->Divide(2,2);

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

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

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

    ClassDef(MHSingles, 1)
};

MStatusDisplay *d = new MStatusDisplay;

MSingles *fSingles;

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()
{
    return kTRUE;
}

Int_t PostProcess()
{
    return kTRUE;
}

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

int recalcSinglepe(const char *file, const char *outfile="", TString option="")
{
    bool save = 0;

    // ======================================================
    //save options

    if ( option.Contains("S") ){
        save = 1;
        cout << "...will save" << endl;
    }

    //Maximum pe order to which the sepctrum will be fitted
    Int_t maxOrder = 5;

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

    cout << file << endl;
    TFile fi(file, "READ");
    if (!fi.IsOpen())
    {
        cout << "ERROR - Could not find file " << file << endl;
        return 2;
    }

    MStatusArray arr;
    if (arr.Read()<=0)
    {
        cout << "ERROR - could not read MStatusArray." << endl;
        return 2;
    }

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

    // Load histogramm for fluctuations of all Pixel
    TH1 *fluct = (TH1*)arr.FindObjectInCanvas("", "TH1F", "Fluct");
    if (!fluct)
    {
        cout << "WARNING - Reading of fluct failed." << endl;
        return 2;
    }

    d->AddTab("Fluct");
    fluct->DrawCopy();


    MHSingles* singles = (MHSingles*) arr.FindObjectInCanvas("MHSingles", "MHSingles", "MHSingles");
    if (!singles)
    {
        cout << "WARNING - Reading of singles failed." << endl;
        return 2;
    }

    d->AddTab("MHSingles");
    singles->DrawCopy("colz");

    TH1 *time = (TH1*)arr.FindObjectInCanvas("Time_py", "TH1D", "Time");
    if (!fluct)
    {
        cout << "WARNING - Reading of time failed." << endl;
        return 2;
    }

    d->AddTab("Time");
    time->DrawCopy();

    TH1 *pulse = (TH1*)arr.FindObjectInCanvas("Pulse_py", "TH1D", "Pulse");
    if (!fluct)
    {
        cout << "WARNING - Reading of pulse failed." << endl;
        return 2;
    }

    d->AddTab("Pulse");
    pulse->DrawCopy();

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

    // Load histogramm for Sum Spectrum of all Pixel
    TH1 *hSum = (TH1*)arr.FindObjectInCanvas("Sum1", "TH1F", "SumHist");
    if (!hSum)
    {
        cout << "WARNING - Reading of hsum failed." << endl;
        return 2;
    }

    // Load fit function for Sum Spectrum of all Pixel
    TF1 *func1 = (TF1*)arr.FindObjectInCanvas("spektrum", "TF1", "SumHist");
    if (!func1)
    {
        cout << "WARNING - Reading of func failed." << endl;
        return 2;
    }

    //Draw histogramm for Sum Spectrum of all Pixel
//    d->AddTab("SumHist");
//    gPad->SetLogy();
//    gPad->SetGrid();
//    hSum->DrawCopy("hist");
//    func1->DrawCopy("SAME");

    // Load TH2 histogramm for Pixel Spectra
    TH2F *hsignal = (TH2F*)arr.FindObjectInCanvas("Signal", "TH2F", "MHSingles");
    if (!hsignal)
    {
        cout << "WARNING - Reading of hSpectra failed." << endl;
        return 2;
    }

    // Draw TH2 histogramm for Pixel Spectra
    d->AddTab("Signal");
    hsignal->SetYTitle("Counts [a.u.]");
    hsignal->DrawCopy("colz");


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

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



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

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

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

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

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

    gLog << endl;
    gLog.Separator("recalculation of gain from pixel spectra");

    MTaskList tlist1;

    MParList plist1;
    plist1.AddToList(&tlist1);
    plist1.AddToList(&badpixels);

    MEvtLoop loop1(file);
    loop1.SetDisplay(d);
    loop1.SetParList(&plist1);

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

    MRawFitsRead read1;
    read1.LoadMap(map);

    MGeomApply apply1;

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

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

    tlist1.AddToList(&mytask);

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

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

    gStyle->SetOptFit(kTRUE);

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

    MGeomCamFACT fact;
    MHCamera hcam(fact);
    MHCamera hcam2(fact);
    MHCamera hNormCam(fact);

    // ======================================================
    // create histograms to plot spectra and parameter distributions

    TH1F hAmplitude ("Rate",      "Average number of dark counts per event",
                     200,  0,  20);
    TH1F hGain      ("Gain",      "Gain distribution",
                     50,  150,   350);
    TH1F hGain2     ("NormGain",      "Normalized gain distribution",
                     51,  0.5, 1.5);
    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 hNoise     ("Noise",    "Noise distribution",
                     80, -10,  30);
    TH1F hXtalkPot  ("XtalkPot",   "Crosstalk Potential Distribution",
                     200,   0,  2);

    hAmplitude.SetXTitle("Amplitude [a.u.]");
    hGain.SetXTitle("gain [a.u.]");
    hGain2.SetXTitle("gain [a.u.]");
    hGainRMS.SetXTitle("gainRMS [a.u.]");
    //hGainRMS.SetStats(false);
    hCrosstlk.SetXTitle("crosstalk [a.u.]");
    hOffset.SetXTitle("baseline offset [a.u.]");
    hFitProb.SetXTitle("FitProb p-value [a.u.]");
    //hFitProb.SetStats(false);

    TH1F hSumScale("Sum2", "Signal spectrum of all pixels",  100, 0.05, 10.05);
    hSumScale.SetXTitle("pulse integral [pe]");
    hSumScale.SetYTitle("Rate");
    hSumScale.SetStats(false);
    hSumScale.Sumw2();

    // define fit function for Amplitudespectrum
    TF1 func("spektrum", fcn, 0, hSum->GetXaxis()->GetXmax(), 7);
    func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "Noise", "XTalkPot");
    func.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;

    // ======================================================
    //-------------------_fit sum spectrum-------------------

    gLog << endl;
    gLog.Separator("Fitting Sum Spectrum");

    //Set Status Display Output and Prograssbar position
    d->SetStatusLine( "...fitting Sum Spectrum", 0);
    d->SetProgressBarPosition( 0.0 , 1);

    const Int_t    maxbin   = hSum->GetMaximumBin();
    const Double_t maxpos   = hSum->GetBinCenter(maxbin);
    const Double_t binwidth = hSum->GetBinWidth(maxbin);
    const Double_t ampl     = hSum->GetBinContent(maxbin);

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

    // full width half maximum of sumspectra
    double fwhmSum = 0;

    //Calculate full width half Maximum
    for (int i=1; i<maxbin; i++)
        if (hSum->GetBinContent(maxbin-i)+hSum->GetBinContent(maxbin+i)<ampl*2)
        {
            fwhmSum = (2*i+1)*binwidth;
            break;
        }

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

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

    //Set fit parameters
    Double_t sum_par[7] =
    {
        ampl, maxpos, fwhmSum*1.4, 0.1, -1, 0, 1
    };

    Double_t sum_par_err[7];

    Double_t fitmin = maxpos-fwhmSum*3;
    Double_t fitmax = hSum->GetXaxis()->GetXmax();

    //set status display prograssbar position
    d->SetProgressBarPosition(0.05 , 1);

    //Fit and draw spectrum
    func.SetParameters(sum_par);
    func.SetRange(fitmin, fitmax);
    hSum->Fit(&func, "INOQSR", "", fitmin, fitmax);
    func.GetParameters(sum_par);

    //Save Parameter Errors for output
    for (int i = 0; i < 7; i++){
        sum_par_err[i] = func.GetParError(i);
    }

    //Set Parameters for following fits
    const float  gain     = sum_par[1];
    const float  GainRMS  = sum_par[2];
    const float  Crosstlk = sum_par[3];
    const float  Offset   = sum_par[4];
    const float  Noise    = sum_par[5];
    const float  XtalkPot = sum_par[6];

    cout << "--------------------" << endl;
    cout << "MaxPos: " << maxpos   << "  \t +/- " << sum_par_err[0] << endl;
    cout << "FWHM:   " << fwhmSum  << endl;
    cout << "GAIN:   " << gain     << "  \t +/- " << sum_par_err[1] << endl;
    cout << "RMS:    " << GainRMS  << "  \t +/- " << sum_par_err[2] << endl;
    cout << "xTalk:  " << Crosstlk << "  \t +/- " << sum_par_err[3] << endl;
    cout << "XtPot:  " << XtalkPot << "  \t +/- " << sum_par_err[6] << endl;
    cout << "Noise:  " << Noise    << "  \t +/- " << sum_par_err[5] << endl;
    cout << "Offset: " << Offset   << "  \t +/- " << sum_par_err[4] << endl;
    cout << "--------------------" << endl;

    gROOT->SetSelectedPad(0);

    //set status display prograssbar position
    d->SetProgressBarPosition(0.1 , 1);

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

    TCanvas &c11 = d->AddTab("SumHist");
    c11.cd();
    gPad->SetLogy();
    gPad->SetGridx();
    gPad->SetGridy();
    hSum->DrawCopy("hist");
    func.DrawCopy("SAME");

    // ======================================================
    //----------Gain Calculation for each Pixel--------------

    gLog << endl;
    gLog.Separator("pixelwise gain calculation");

    // counter for number of processed pixel
    int     count_ok = 0;

    // marker for pixel with abnormal behavior
    bool    suspicous = 0;

    // number of pixels to process
    UInt_t  nPixels = hsignal->GetNbinsX();

    //set status display status
    d->SetStatusLine( "...fitting pixels", 0);

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

    // Loop over Pixels
    for (UInt_t pixel = 0; pixel < nPixels; pixel++)
    {
        //set status display prograssbar position
        d->SetProgressBarPosition(0.1 + 0.8*( (double)(pixel+1)/(double)nPixels ), 1);
        suspicous = 0;

        // jump to next pixel if the current is marked as broken
        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);

        //Rebin Projection
        hist->Rebin(2);

        // Callculate values of the bin with maximum entries
        const Int_t     maxbin   = hist->GetMaximumBin();
        const Double_t  maxPos   = hist->GetBinCenter(maxbin);

        // Calculate range for the fit
        const double    fit_min = maxPos-GainRMS*2;
        const double    fit_max = maxPos+gain*(maxOrder-0.5);

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

        // Set fit parameter start values
        sum_par[0] = hist->GetBinContent(maxbin);
//        sum_par[2] = first_gaus_RMS;
//        sum_par[3] = second_gaus_ampl/first_gaus_ampl;
//        sum_par[4] = -1*(second_gaus_ampl-2*first_gaus_ampl);
        func.SetParameters(sum_par);
//        func.FixParameter(0, first_gaus_ampl);
//        func.SetParLimits(0, 0.8*sum_par[0], 1.1*sum_par[0]);
//        func.SetParLimits(1, maxPos - 2*GainRMS, maxPos + 2*GainRMS);
//        func.SetParLimits(2, -first_gaus_RMS, +first_gaus_RMS);
//        func.FixParameter(5, sum_par[5]);
//        func.SetParLimits(5, 0.7, 1.1);
//        func.FixParameter(6, sum_par[6]);
//        func.SetParLimits(6, 0.1, 1.9);

        // ======================================================
        // -------------- Fit Pixels spectrum -------------------

        // Save fit result to a variable
        const TFitResultPtr rc = hist->Fit(&func, "LLN0QS", "", fit_min, fit_max);

        // marker to mark if fit was successfull or not
        const bool ok = int(rc)==0;

        // Save Parametervalues for statistics and bug fixing
        const double fit_prob   = rc->Prob();
//        const float  fMaxAmpl   = func.GetParameter(0);
        const float  fGain      = func.GetParameter(1);
        const float  fAmplitude = func.Integral(Offset, maxOrder*fGain+Offset); //func.GetParameter(0);
        const float  f2GainRMS  = func.GetParameter(2);
        const float  fCrosstlk  = func.GetParameter(3);
        const float  fOffset    = func.GetParameter(4)/30;
        const float  fNoise     = func.GetParameter(5);
        const float  fXtlkPot   = func.GetParameter(6);

        // Fill Parameter value hgistograms
        hAmplitude.Fill(fAmplitude);
        hGain.Fill(fGain);
        hGain2.Fill(fGain/gain);
        hCrosstlk.Fill(fCrosstlk*100);
        hOffset.Fill(fOffset);
        hFitProb.Fill(100*fit_prob);
        hGainRMS.Fill(100*f2GainRMS/fGain);
        hNoise.Fill(fNoise);
        hXtalkPot.Fill(fXtlkPot);

        // Plot calculated Gain values to Camera Display
        hcam.SetBinContent(pixel+1, fGain);
        hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);
        hNormCam.SetBinContent(pixel+1, fGain/gain);

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

        // cancel out pixel where the fit was not succsessfull
        usePixel[pixel] = ok;

        // mark pixels suspicious with negative noise
        if ( fNoise < 0)
            suspicous = 1;

        // mark pixels suspicious with negative GainRMS
        if ( f2GainRMS < 0)
            suspicous = 1;

        // mark pixels suspicious with with large deviation from mean gain
        if ( fGain < gain - 1.4*GainRMS || fGain > gain + 1.4*GainRMS)
            suspicous = 1;

        //plott especialy marked pixel
        if (count_ok<3 || suspicous == 1 || !ok)
        {
            hist->SetName("Pix%d");
            TCanvas &c = d->AddTab(Form("Pix%d", pixel));
            c.cd();
            gPad->SetLogy();
            gPad->SetGridx();
            gPad->SetGridy();

            rc->Print("V");

//            hist->SetStats(false);
            hist->SetYTitle("Counts [a.u.]");
            hist->SetName(Form("Pix%d", pixel));
            hist->DrawCopy()->SetDirectory(0);
            func.DrawCopy("SAME")->SetRange(fit_min, fit_max);
        }

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

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

        //Fill sum spectrum
        for (int b=1; b<=hist->GetNbinsX(); b++)
            hSumScale.Fill(
                        (hist->GetBinCenter(b)-fOffset)/fGain,
                        hist->GetBinContent(b)
                        );

        count_ok++;

        delete hist;
    }

    // define a general gaus function
    TF1 fgaus("fgaus", "gaus");
    hGain.Fit(&fgaus, "N0QS");

    // ======================================================
    //------------------Draw histograms----------------------

    gLog << endl;
    gLog.Separator("Drawing Histograms");

    d->SetStatusLine( "...drawing Histograms", 0);
    d->SetProgressBarPosition(0.9 , 1);

    // cancel out pixel in camera view
    hcam.SetUsed(usePixel);
    hcam2.SetUsed(usePixel);

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

    TCanvas &canv = d->AddTab("Gain");
    canv.cd();
    canv.Divide(2,2);

    canv.cd(1);
    hcam.SetNameTitle("gain", "Gain distribution over whole camera");
    hcam.DrawCopy();

    canv.cd(2);
    hcam2.SetNameTitle("gainRMS", "GainRMS distribution over whole camera");
    hcam2.DrawCopy();

    TCanvas &cam_canv = d->AddTab("Camera_Gain");

    // ======================================================
    //---------- Draw gain corrected sum specetrum ----------

    gROOT->SetSelectedPad(0);
    TCanvas &c12 = d->AddTab("GainHist");
    c12.cd();
    gPad->SetLogy();
    gPad->SetGridx();
    gPad->SetGridy();

    hSumScale.DrawCopy("hist");

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

    // Callculate values of the bin with maximum entries
    const Int_t    maxbin2 = hSumScale.GetMaximumBin();
    const Double_t ampl2   = hSumScale.GetBinContent(maxbin2);

    //Set fit start parameters
    Double_t par2[7] =
    {
        ampl2, 1, GainRMS/fgaus.GetParameter(1), Crosstlk, 0, 0, 1
    };
    func.SetParameters(par2);

    // fix gain to 1 as it was normalized allready
    func.FixParameter(1, 1);

    // fix offset parameter
    // func.FixParameter(4, 0);

    // set progressbar possition
    d->SetStatusLine( "...fitting corr. spectrum", 0);

    // Fit corrected Spectrum with Likelyhood Method
    const TFitResultPtr rc = hSumScale.Fit(&func, "LLEN0QS", "", 0.75, 10);

    func.DrawCopy("SAME")->SetRange(0.5, 10);

    // Print fit results
    rc->Print("V");

    // set progressbar possition
    d->SetProgressBarPosition(0.95 , 1);

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

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

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

    //plot normailzed camera view
    cam_canv.cd();
    hNormCam.SetStats(false);

    hNormCam.SetNameTitle("GainCam", "Normalized gain distribution");
    hNormCam.SetUsed(usePixel);
    hNormCam.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();
    hNoise.DrawCopy();

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

    TCanvas &c3 =d->AddTab("Prop");
    c3.Divide(2,1);
    c3.cd(1);
    gPad->SetGrid();
    gPad->SetLogy();
    hFitProb.DrawCopy();

    c3.cd(2);
    gPad->SetGrid();
    gPad->SetLogy();
    hAmplitude.DrawCopy();

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

    d->AddTab("Norm");
    gPad->SetGrid();
    gPad->SetLogy();
    TH1 *hh = hGain2.DrawCopy();
    hh->Fit("gaus");

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

    //save results to file
    if (save){
    d->SetStatusLine( "...saving results to rootfile", 0);
    d->SaveAs(outfile);
    cout << "..success!" << endl;
    }
    d->SetProgressBarPosition(1 , 1);
    d->SetStatusLine( "...done", 0);
    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 noise      = par[5];
    Double_t crossPot   = par[6];

    Double_t xTalk = 1;

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

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

        y += xTalk*gauss;
        xTalk = Power(cross, Power(order, crossPot) );
    }

    return ampl*y;
}

// end of PlotMedianEachSliceOfPulse
//----------------------------------------------------------------------------