#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 <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 "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
{
    TH2F       fSignal;     //changed from MH2F to TH2F
    TH2F       fTime;       //changed from MH2F to TH2F
    TProfile2D fPulse;      //changed from MProfile2D to TProfile2D


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

    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);
Double_t fcn_cross(Double_t *x, Double_t *par);

void FitGaussiansToSpectrum( UInt_t nfunc,  TH1 &hSource,   TH1 &hDest,
                             Int_t maxbin,  double fwhm, Double_t gain );

Double_t MedianOfH1 (
  TH1*            inputHisto
);

int singlepe(
        const char *runfile, const char *drsfile, const char *outfile
)
{

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

//    const char *drsfile = "/fact/raw/2012/06/25/20120625_112.drs.fits.gz";

  MDirIter iter;
    iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(runfile));

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

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

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

  MEvtLoop loop1(Title);
  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->SetYTitle("Counts [a.u.]");
  ht->Scale(1./1440);
  ht->Draw();

  d->AddTab("Pulse");
  TH1 *hp = hpulse->ProjectionY();
  hp->SetYTitle("Counts [a.u.]");
  hp->Scale(1./1440);
  hp->Draw();


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


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

  //built an array of Amplitude spectra
  TH1F hAmplitude ("Rate",      "Average number of dark counts per event",
                   200,  0,  20);
  hAmplitude.SetXTitle("Amplitude [a.u.]");
  hAmplitude.SetYTitle("Rate");

  TH1F hGain      ("Gain",      "Gain distribution",
                   100,  0,   400);
  hGain.SetXTitle("gain [a.u.]");
  hGain.SetYTitle("Rate");

  TH1F hGainRMS   ("RelSigma",  "Distribution of rel. Sigma",
                   100,   0,   30);
  hGainRMS.SetXTitle("gainRMS [a.u.]");
  hGainRMS.SetYTitle("Rate");
  hGainRMS.SetStats(false);

  TH1F hCrosstlk  ("Crosstalk", "Distribution of crosstalk probability",
                   100,   0,    25);
  hCrosstlk.SetXTitle("crosstalk [a.u.]");
  hCrosstlk.SetYTitle("Rate");

  TH1F hOffset    ("Offset",    "Baseline offset distribution",
                   50, -7.5,  2.5);
  hOffset.SetXTitle("baseline offset [a.u.]");
  hOffset.SetYTitle("Rate");

  TH1F hFitProb   ("FitProb",   "FitProb distribution",
                   50,   0,  100);
  hFitProb.SetXTitle("FitProb p-value [a.u.]");
  hFitProb.SetYTitle("Rate");
  hFitProb.SetStats(false);


  TH1F hSum       ("Sum1",      "Sum of all pixel's' integral spectra",
                   100,    0,  2000);
  hSum.SetXTitle("pulse integral [a.u.]");
  hSum.SetYTitle("Rate");
  hSum.SetStats(false);
  hSum.Sumw2();


  TH1F hSumScale  ("Sum2",
                   "Sum of all pixel's' integral spectra (scaled with gain)",
                   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, 1000, 5);
  func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
  func.SetLineColor(kBlue);

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

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

  /*
  TF1 funcScaledBL("gain_sum_spektrum", fcn, 0, 10, 5);
  funcScaledBL.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
  funcScaledBL.SetLineColor(kRed);
  */


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

//    usePixel[990] = 0;

  // Map for which pixel which were suspicous
  bool suspicous[1440];

  for (int i=0; i<1440; i++) suspicous[i]=false;

  // List of Pixel that should be ignored in camera view
  //    suspicous[802]       = true;

  //------------------------------------------------------------------------
  //  Fill and calculate sum spectrum

  // Begin of Loop over Pixels
  for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
    {
      //jump to next pixel if the current is marked as faulty
      if (usePixel[pixel]==0)
        continue;

      TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
      hist->SetDirectory(0);

      //loop over slices
      for (int b=1; b<=hist->GetNbinsX(); b++)
        {
          //Fill sum spectrum
          hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
        }
    }

  for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
    {
      hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
    }

  gROOT->SetSelectedPad(0);
  TCanvas &c11 = d->AddTab("SumHist");
  c11.cd();
  gPad->SetLogy();
  gPad->SetGridx();
  gPad->SetGridy();
  hSum.DrawCopy("hist");
  //------------------------fit sum spectrum-------------------------------
  const Int_t    maxbin   = hSum.GetMaximumBin();
  const Double_t binwidth = hSum.GetBinWidth(maxbin);
  const Double_t ampl     = hSum.GetBinContent(maxbin);

  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[6] =
  {
    ampl, hSum.GetBinCenter(maxbin), fwhmSum*1.4, 0.1, -10, 1
  };
  funcSum.SetParameters(sum_par);

  //calculate fit range
  const double min = sum_par[1]-fwhmSum*2;
  const double max = sum_par[1]*(maxOrder);
  funcSum.SetRange(min, max);
  funcSum.FixParameter(5,0.1);

  //Fit and draw spectrum
  hSum.Fit(&funcSum, "NOQS", "", min, max);
  funcSum.DrawCopy("SAME");
  funcSum.GetParameters(sum_par);

  //Set Parameters for following fits
//    const float  Amplitude      = sum_par[0];
  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  PowCrosstlk    = sum_par[5];


  //--------fit gausses to peaks of sum specetrum -----------

  // create histo for height of Maximum amplitudes vs. pe
  double min_bin  =   hSum.GetBinCenter(maxbin);
  double max_bin  =   hSum.GetBinCenter((maxOrder-2)*(maxbin));
  double bin_width=   (max_bin - min_bin)/10;

  TH1D hMaxHeightsSum("MaxHeightSum",      "peak heights",
                      10,  min_bin-bin_width, max_bin*2-bin_width/2 );

  FitGaussiansToSpectrum
  (
    maxOrder-2,
    hSum,
    hMaxHeightsSum,
    maxbin,
    fwhmSum,
    gain
  );

  //EOF: fit sum spectrum-----------------------------

  hMaxHeightsSum.SetLineColor(kRed);
//    hMaxHeightsSum.DrawCopy("SAME");

  //Fit Peak heights
  TF1 fExpo1( "fExpo1","expo", min, max);
  fExpo1.SetLineColor(kRed);
  hMaxHeightsSum.Fit(&fExpo1, "N0QS" );
//    hMaxHeightsSum.DrawCopy("SAME");
//    fExpo1.DrawCopy("SAME");

  //  EOF:    Fill and calculate sum spectrum
  //------------------------------------------------------------------------



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

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

      for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
        {
          hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1);
        }

      const Int_t    maxbin   = hist->GetMaximumBin();

      const Double_t binwidth = hist->GetBinWidth(maxbin);

      const Double_t ampl     = hist->GetBinContent(maxbin);

      const Double_t maxPos   = hist->GetBinCenter(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),
//            gain,
        fwhm*2,
//            GainRMS,
        Crosstlk,
        Offset
      };


      func.SetParameters(par);
      const double fit_min = par[1]-fwhm*1.4;
      const double fit_max = par[1]*maxOrder;
      func.SetRange(fit_min-fwhm*2, fit_max);
      func.SetParLimits(1, maxPos - fwhm, maxPos + fwhm);

      if (suspicous[pixel] == true)
        {
          cout << "Maxbin\t" << maxbin << endl;
          cout << "min\t" << fit_min << endl;
          cout << "max\t" << fit_max << endl;
          cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
               << par[0] << "\t"
               << par[1] << "\t"
               << par[2] << "\t"
               << fwhm << "\t" << endl;
        }

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

      //Fit Pixels spectrum
      const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_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();
//        const float  fCrossOrd  = func.GetParameter(5);

      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);
      hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);

      usePixel[pixel] = ok;

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

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

      //plott special pixel
      if (
        count_ok==0 ||
        suspicous[pixel]
      )
        {
          TCanvas &c = d->AddTab(Form("Pix%d", pixel));
          c.cd();
          gPad->SetLogy();
          gPad->SetGridx();
          gPad->SetGridy();

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

      count_ok++;

      delete hist;
//        if (suspicous[pixel] == true) usePixel[pixel]=0;
    }

  // End of Loop over Pixel

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

  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();

  for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
    {
      hSumScale.SetBinError(bin, hSumScale.GetBinContent(bin)*0.1);
    }

  hSumScale.DrawCopy("hist");


  //-------------------- fit gain corrected sum spectrum -------------------
  const Int_t    maxbin2   = hSumScale.GetMaximumBin();
  const Double_t binwidth2 = hSumScale.GetBinWidth(maxbin2);
  const Double_t ampl2     = hSumScale.GetBinContent(maxbin2);

  //Calculate full width half Maximum
  double fwhmScaled = 0;

  for (int i=1; i<maxbin; i++)
    if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
      {
        fwhmScaled = (2*i+1)*binwidth2;
        break;
      }

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

  //Set fit parameters
  Double_t par2[6] =
  {
    ampl2, 0.9, fwhmScaled*1.4, Crosstlk, 0, 1.1
  };

  funcScaled.SetParameters(par2);
//  funcScaledBL.SetParameters(par2);
//    funcScaled.FixParameter(1,0.9);
  funcScaled.FixParameter(4,0);
//  funcScaledBL.FixParameter(4,0);
  funcScaled.FixParameter(5,Crosstlk);
//  funcScaledBL.FixParameter(5,Crosstlk);

  const double minScaled = par2[1]-fwhmScaled;
  const double maxScaled = par2[1]*maxOrder;
  cout << "par2[1]" << par2[1] << endl;
  cout << "maxScaled\t"           << maxScaled
       << " \t\t minScaled\t"     << minScaled << endl;

  funcScaled.SetRange(minScaled, maxScaled);
//  funcScaledBL.SetRange(minScaled, maxScaled);
  hSumScale.Fit(&funcScaled, "N0QS", "", minScaled, maxScaled);
//  hSumScale.Fit(&funcScaledBL, "WLN0QS", "", minScaled, maxScaled);
  funcScaled.DrawCopy("SAME");
//  funcScaledBL.DrawCopy("SAME");
  //--------fit gausses to peaks of gain corrected sum specetrum -----------

  // create histo for height of Maximum amplitudes vs. pe
  TH1D hMaxHeights("MaxHeights",      "peak heights",
                   10,  hSumScale.GetBinCenter(maxbin-1)-0.5,   10+hSumScale.GetBinCenter(maxbin-1)-0.5);

  FitGaussiansToSpectrum
  (
    maxOrder-2,
    hSumScale,
    hMaxHeights,
    maxbin2,
    fwhmScaled,
    1
  );



//    TCanvas &c16 = d->AddTab("PeakHeights");
//    gPad->SetLogy();
  hMaxHeights.SetLineColor(kRed);
//    hMaxHeights.DrawCopy("SAME");
  TF1 fExpo( "fExpo","expo", 0,  maxOrder-2);
  fExpo.SetLineColor(kRed);
  hMaxHeights.Fit(&fExpo, "N0QS" );
//    fExpo.DrawCopy("SAME");
//    c12.cd();
//    fExpo.DrawCopy("SAME");

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

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


  TF1 GainGaus( "GainGaus", "gaus", 0, hGain.GetXaxis()->GetXmax());
  hGain.Rebin(2);
  hGain.Fit(&GainGaus, "N0QS");

//    float gain_mean     = GainGaus.GetParameter(1);
  float gain_median   = MedianOfH1(&hGain);
  TH1F hNormGain          ("NormGain",      "median normalized Gain distribution",
                           hGain.GetNbinsX(),
                           (hGain.GetXaxis()->GetXmin())/gain_median,
                           (hGain.GetXaxis()->GetXmax())/gain_median
                          );
  hNormGain.SetXTitle("gain [fraction of median gain value]");
  hNormGain.SetYTitle("Counts");


  hNormGain.Add(&hGain);
//    hNormGain.SetStats(false);
  hNormGain.GetXaxis()->SetRangeUser(
    1-(hNormGain.GetXaxis()->GetXmax()-1),
    hNormGain.GetXaxis()->GetXmax()
  );
  hNormGain.DrawCopy();
  hNormGain.Fit(&GainGaus, "N0QS");
  GainGaus.DrawCopy("SAME");

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

  for (UInt_t pixel =1; pixel <= 1440; pixel++)
    {
      hNormCam.SetBinContent(pixel, hcam.GetBinContent(pixel)/gain_median);
    }

  hNormCam.SetNameTitle( "gain","normalized [Median] Gain distribution over whole camera");
  hNormCam.SetZTitle("Horst");
  hNormCam.SetUsed(usePixel);
  hNormCam.DrawCopy();
//    hGain.Scale(1/GainGaus.GetParameter(1));
//    GainGaus.DrawCopy("SAME");

  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();
  hGain.DrawCopy();

/*
  TCanvas &c25 = d->AddTab("Spectra");
  c25.Divide(2,1);
  c25.cd(1);
  gPad->SetLogy();
  gPad->SetGrid();
  hSum.DrawCopy("hist");
  funcSum.DrawCopy("SAME");

  c25.cd(2);
  gPad->SetLogy();
  gPad->SetGrid();
  hSumScale.DrawCopy("hist");
  funcScaled.DrawCopy("SAME");
  funcScaledBL.DrawCopy("SAME");
*/

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

  d->SaveAs(outfile);
  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;
}

Double_t fcn_cross(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 powOfX = par[5];

  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 = pow(cross, order*powOfX);
//        for (int j = 1; j < powOfX; j++)
//        {
//            xTalk *= cross;
//        }
//        cross *= TMath::Exp(-powOfX*cross);
      xTalk *= cross;
    }

  return ampl*y;
}


void
FitGaussiansToSpectrum(
  UInt_t      maxOrder,
  TH1        &hSource,
  TH1        &hDest,
  Int_t       maxbin,
  double      fwhm,
  Double_t    gain
)
{

  Double_t    peakposition = hSource.GetBinCenter(maxbin);
  char        fname[20];

  // fit gauss functions to spectrum
  for (UInt_t nr = 1; nr<maxOrder; nr++)
    {
      sprintf(fname,"gaus%d",nr);

      //create gauss functions
      TF1 gaus( fname,"gaus", peakposition-fwhm,     peakposition+fwhm);
      gaus.SetParLimits(1, peakposition-fwhm, peakposition+fwhm);
      gaus.SetParLimits(2, -fwhm, fwhm );
      //fit and draw gauss functions
      hSource.Fit(  &gaus,   "N0QS", "", peakposition-fwhm, peakposition+fwhm);
//        gaus.DrawCopy("SAME");

      peakposition = gaus.GetParameter(1);

      //fill histo for height of Maximum amplitudes vs. pe
      hDest.SetBinContent(nr, gaus.GetParameter(0));

      peakposition += gain;
    }

  return;
}

Double_t MedianOfH1 (
  TH1*            inputHisto
)
{
  //compute the median for 1-d histogram h1
  Int_t nbins = inputHisto->GetXaxis()->GetNbins();
  Double_t *x = new Double_t[nbins];
  Double_t *y = new Double_t[nbins];

  for (Int_t i=0; i<nbins; i++)
    {
      x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
      y[i] = inputHisto->GetBinContent(i+1);
    }

  Double_t median = TMath::Median(nbins,x,y);
  delete [] x;
  delete [] y;
  return median;
}
// end of PlotMedianEachSliceOfPulse
//----------------------------------------------------------------------------