Index: /fact/tools/marsmacros/singlepe.C
===================================================================
--- /fact/tools/marsmacros/singlepe.C	(revision 14229)
+++ /fact/tools/marsmacros/singlepe.C	(revision 14230)
@@ -8,4 +8,9 @@
 #include <TFitResultPtr.h>
 #include <TFitResult.h>
+#include <Getline.h>
+
+#include <cstdio>
+#include <stdio.h>
+#include <stdint.h>
 
 #include "MH.h"
@@ -407,4 +412,8 @@
 
 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 );
 
 int singlepe()
@@ -412,8 +421,14 @@
     // ======================================================
 
-    const char *drsfile = "/fact/raw/2012/06/01/20120601_003.drs.fits.gz";
+    const char *drsfile = "/fact/raw/2012/03/09/20120309_012.drs.fits.gz";
+//    const char *drsfile = "/fact/raw/2012/01/13/20120113_003.drs.fits.gz";
 
     MDirIter iter;
-    iter.AddDirectory("/fact/raw/2012/06/01", "20120601_005.fits.gz");
+//    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_014.fits.gz");
+//    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_032.fits.gz");
+//    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_038.fits.gz");
+    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_017.fits.gz");
+//    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_018.fits.gz");
+//    iter.AddDirectory("/fact/raw/2012/01/13", "20120113_007.fits.gz");
 
     // ======================================================
@@ -566,4 +581,6 @@
     // AFTER this the Code for the spektrum fit follows
     // access the spektrumhistogram by the pointer *hist
+    UInt_t maxOrder     = 6;
+
 
     MGeomCamFACT fact;
@@ -571,20 +588,21 @@
 
     //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",
+    TH1F hAmplitude ("Rate",      "Average number of dark counts per event",
+                     200,  0,  20);
+    TH1F hGain      ("Gain",      "Gain distribution",
+                     50,  100,   300);
+    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);
+    TH1F hSum       ("Sum1",      "Sum of all pixel amplitude Spectra",
+                     100,    0,  2000);
+    TH1F hSumScale  ("Sum2",
+                     "Sum of all pixel amplitude Spectra (scaled with gain)",
+                     100,    0,   10);
 
     // define fit function for Amplitudespectrum
@@ -594,6 +612,6 @@
 
     // define fit function for Amplitudespectrum
-    TF1 func2("sum_spektrum", fcn, 0, 2000, 5);
-    func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
+    TF1 func2("sum_spektrum", fcn_cross, 0, 2000, 6);
+    func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
     func2.SetLineColor(kBlue);
 
@@ -623,9 +641,141 @@
     usePixel[1393] = 0;
 
-    int count_ok = 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;
+    suspicous[543]     = true;
+    suspicous[465]     = true;
+    suspicous[1154]     = true;
+    suspicous[342]     = true;
+    suspicous[1319]      = true;
+    suspicous[1401]      = true;
+    suspicous[1318]      = true;
+    suspicous[1400]     = true;
+    suspicous[243]      = true;
+    suspicous[1299]     = true;
+    suspicous[764]     = true;
+    suspicous[762]     = true;
+    suspicous[1427]      = true;
+    suspicous[1319]     = true;
+    suspicous[118]     = 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));
+        }
+    }
+
+    gROOT->SetSelectedPad(0);
+    TCanvas &c11 = d->AddTab("SumHist");
+    gPad->SetLogy();
+    hSum.DrawCopy();
+    //------------------------fit sum spectrum-------------------------------
+    const Int_t    maxbin   = hSum.GetMaximumBin();
+    const Double_t binwidth = hSum.GetBinWidth(maxbin);
+    const Double_t ampl     = hSum.GetBinContent(maxbin);
+
+    double fwhm2 = 0;
+
+    //Calculate full width half Maximum
+    for (int i=1; i<maxbin; i++)
+        if (hSum.GetBinContent(maxbin-i)+hSum.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 sum_par[6] =
+    {
+        ampl, hSum.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10, 1
+    };
+    func2.SetParameters(sum_par);
+
+    //calculate fit range
+    const double min = sum_par[1]-fwhm2*1.4;
+    const double max = sum_par[1]*(maxOrder+1);
+    func2.SetRange(min-2*fwhm2, max);
+
+    //Fit and draw spectrum
+    hSum.Fit(&func2, "NOS", "", min, max);
+    func2.DrawCopy("SAME");
+    func2.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*(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,
+                hSum,
+                hMaxHeightsSum,
+                maxbin,
+                fwhm2,
+                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, "N0S" );
+    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;
@@ -638,4 +788,5 @@
         const Double_t binwidth = hist->GetBinWidth(maxbin);
         const Double_t ampl     = hist->GetBinContent(maxbin);
+        const Double_t maxPos   = hist->GetBinCenter(maxbin);
 
         double fwhm = 0;
@@ -657,14 +808,33 @@
         Double_t par[5] =
         {
-            ampl, hist->GetBinCenter(maxbin), fwhm*1.4, 0.1, -10
+            ampl,
+            hist->GetBinCenter(maxbin),
+//            gain,
+            fwhm*2,
+//            GainRMS,
+            Crosstlk,
+            Offset
         };
 
+
         func.SetParameters(par);
-
-        const double min = par[1]-fwhm*1.4;
-        const double max = par[1]*5;
-
+        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;
+        }
         //Fit Pixels spectrum
-        const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", min, max);
+        const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max);
 
         const bool ok = int(rc)==0;
@@ -676,4 +846,14 @@
         const float  fCrosstlk  = func.GetParameter(3);
         const float  fOffset    = func.GetParameter(4)/singles.GetIntegrationWindow();
+//        const float  fCrossOrd  = func.GetParameter(5);
+
+        if (suspicous[pixel] == true)
+        {
+            cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
+                 << fAmplitude << "\t"
+                 << fGain << "\t"
+                 << f2GainRMS << "\t"
+                 << fwhm << "\t" << endl;
+        }
 
         hAmplitude.Fill(fAmplitude);
@@ -695,11 +875,16 @@
         }
 
+        //Fill sum spectrum
         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)
+            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));
@@ -707,13 +892,14 @@
             hist->SetName(Form("Pix%d", pixel));
             hist->DrawCopy()->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-----------------------
@@ -726,92 +912,19 @@
     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
+    //------------------ 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);
+    hSumScale.DrawCopy();
+
+    //-------------------- 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
     fwhm2 = 0;
     for (int i=1; i<maxbin; i++)
-        if (hSum2.GetBinContent(maxbin2-i)+hSum2.GetBinContent(maxbin2+i)<ampl2*2)
+        if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
         {
             fwhm2 = (2*i+1)*binwidth2;
@@ -826,5 +939,5 @@
     Double_t par2[5] =
     {
-        ampl2, hSum2.GetBinCenter(maxbin2), fwhm2*1.4, 0.1, 0
+        ampl2, hSumScale.GetBinCenter(maxbin2), fwhm2*1.4, 0.1, -0.2
     };
 
@@ -832,13 +945,37 @@
 
     const double min2 = par2[1]-fwhm2*1.4;
-    const double max2 = par2[1]*5;
-
-    hSum2.Fit(&func3, "N0QS", "", min2, max2);
+    const double max2 = par2[1]*maxOrder;
+//    func.SetRange(min2-fwhm2, max);
+    hSumScale.Fit(&func3, "N0QS", "", min2, max2);
     func3.DrawCopy("SAME");
-    //-------------------END OF: fit sum spectrum-----------------------------
-
-    TCanvas &c15 = d->AddTab("PeakHeights");
-    gPad->SetLogy();
-    hMaxHeights.DrawCopy();
+
+    //--------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,
+                hSumScale,
+                hMaxHeights,
+                maxbin2,
+                fwhm2,
+                1
+                );
+
+
+
+//    TCanvas &c16 = d->AddTab("PeakHeights");
+//    gPad->SetLogy();
+    hMaxHeights.SetLineColor(kRed);
+//    hMaxHeights.DrawCopy("SAME");
+    TF1 fExpo( "fExpo","expo", 0,  maxOrder-1);
+    fExpo.SetLineColor(kRed);
+    hMaxHeights.Fit(&fExpo, "N0QS" );
+    fExpo.DrawCopy("SAME");
+//    c12.cd();
+//    fExpo.DrawCopy("SAME");
 
     TCanvas &c2 = d->AddTab("Result");
@@ -880,9 +1017,10 @@
     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 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 k      = par[5];
 
     Double_t xTalk = 1;
@@ -897,7 +1035,82 @@
         y += xTalk*gauss;
 
-        xTalk *= cross;
+//        for (int j = 1; j < k; j++)
+//        {
+            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;
+
+//        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;
+}
+