Index: /fact/tools/marsmacros/singlepe.C
===================================================================
--- /fact/tools/marsmacros/singlepe.C	(revision 14241)
+++ /fact/tools/marsmacros/singlepe.C	(revision 14242)
@@ -153,19 +153,19 @@
 
         fSignal.SetName("Signal");
-        fSignal.SetTitle("Title");
-        fSignal.SetXTitle("X title");
-        fSignal.SetYTitle("Y title");
+        fSignal.SetTitle("pulse integral spectrum");
+        fSignal.SetXTitle("pixel [SoftID]");
+        fSignal.SetYTitle("time [a.u.]");
         fSignal.SetDirectory(NULL);
 
         fTime.SetName("Time");
-        fTime.SetTitle("Title");
-        fTime.SetXTitle("X title");
-        fTime.SetYTitle("Y title");
+        fTime.SetTitle("arival time spectrum");
+        fTime.SetXTitle("pixel [SoftID]");
+        fTime.SetYTitle("time [a.u.]");
         fTime.SetDirectory(NULL);
 
         fPulse.SetName("Pulse");
-        fPulse.SetTitle("Title");
-        fPulse.SetXTitle("X title");
-        fPulse.SetYTitle("Y title");
+        fPulse.SetTitle("average pulse shape spectrum");
+        fPulse.SetXTitle("pixel [SoftID]");
+        fPulse.SetYTitle("time [a.u.]");
         fPulse.SetDirectory(NULL);
     }
@@ -417,13 +417,35 @@
                              Int_t maxbin,  double fwhm, Double_t gain );
 
-int singlepe(const char *runfile, const char *drsfile, const char *outfile)
+Double_t MedianOfH1 (
+        TH1*            inputHisto
+        );
+
+int singlepe(
+//        const char *runfile, const char *drsfile, const char *outfile
+        )
 {
+
     // ======================================================
 
-    //const char *drsfile = "/fact/raw/2012/05/18/20120518_012.drs.fits.gz";
+    const char *drsfile = "/fact/raw/2012/05/18/20120518_003.drs.fits.gz";
+//    const char *drsfile = "/fact/raw/2012/06/25/20120625_112.drs.fits.gz";
 
     MDirIter iter;
-    //iter.AddDirectory("/fact/raw/2012/05/18", "20120518_005.fits.gz");
-    iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(runfile));
+    iter.AddDirectory("/fact/raw/2012/05/18", "20120518_005.fits.gz");
+//    iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(runfile));
+
+//    TString filename;
+//    for (UInt_t fileNr = 116; fileNr < 147; fileNr++)
+//    {
+//        filename    = "20120625_";
+//        filename   += fileNr;
+//        filename   += ".fits.gz";
+//        iter.AddDirectory("/fact/raw/2012/06/25", filename);
+//    }
+
+//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_116.fits.gz");
+//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_117.fits.gz");
+//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_118.fits.gz");
+//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_119.fits.gz");
 
     // ======================================================
@@ -570,4 +592,5 @@
     d->AddTab("Time");
     TH1 *ht = htime->ProjectionY();
+    ht->SetYTitle("Counts [a.u.]");
     ht->Scale(1./1440);
     ht->Draw();
@@ -575,4 +598,5 @@
     d->AddTab("Pulse");
     TH1 *hp = hpulse->ProjectionY();
+    hp->SetYTitle("Counts [a.u.]");
     hp->Scale(1./1440);
     hp->Draw();
@@ -582,30 +606,59 @@
     // AFTER this the Code for the spektrum fit follows
     // access the spektrumhistogram by the pointer *hist
-    UInt_t maxOrder     = 6;
+    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",
                      50,  100,   300);
+    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);
-    TH1F hSum       ("Sum1",      "Sum of all pixel amplitude Spectra",
+    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 amplitude Spectra (scaled with gain)",
+                     "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();
 
@@ -616,12 +669,12 @@
 
     // define fit function for Amplitudespectrum
-    TF1 func2("sum_spektrum", fcn_cross, 0, 2000, 6);
-    func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
-    func2.SetLineColor(kBlue);
+    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 func3("gain_sum_spektrum", fcn, 0, 10, 5);
-    func3.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
-    func3.SetLineColor(kBlue);
+    TF1 funcScaled("gain_sum_spektrum", fcn, 0, 10, 5);
+    funcScaled.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
+    funcScaled.SetLineColor(kBlue);
 
 
@@ -673,5 +726,5 @@
         }
     }
-    for(UInt_t bin = 1; bin < hSum.GetNbinsX()+1; bin++)
+    for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
     {
         hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
@@ -680,4 +733,5 @@
     gROOT->SetSelectedPad(0);
     TCanvas &c11 = d->AddTab("SumHist");
+    c11.cd();
     gPad->SetLogy();
     gPad->SetGridx();
@@ -689,5 +743,5 @@
     const Double_t ampl     = hSum.GetBinContent(maxbin);
 
-    double fwhm2 = 0;
+    double fwhmSum = 0;
 
     //Calculate full width half Maximum
@@ -695,8 +749,8 @@
         if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl*2)
         {
-            fwhm2 = (2*i+1)*binwidth;
+            fwhmSum = (2*i+1)*binwidth;
             break;
         }
-    if (fwhm2==0)
+    if (fwhmSum==0)
     {
         cout << "Could not determine start value for sigma." << endl;
@@ -706,25 +760,26 @@
     Double_t sum_par[6] =
     {
-        ampl, hSum.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10, 1
+        ampl, hSum.GetBinCenter(maxbin), fwhmSum*1.4, 0.1, -10, 1
     };
-    func2.SetParameters(sum_par);
+    funcSum.SetParameters(sum_par);
 
     //calculate fit range
-    const double min = sum_par[1]-fwhm2*2;
-    const double max = sum_par[1]*(maxOrder+1);
-    func2.SetRange(min, max);
+    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(&func2, "NOS", "", min, max);
-    func2.DrawCopy("SAME");
-    func2.GetParameters(sum_par);
+    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  Amplitude      = sum_par[0];
     const float  gain           = sum_par[1];
-    const float  GainRMS        = sum_par[2];
+//    const float  GainRMS        = sum_par[2];
     const float  Crosstlk       = sum_par[3];
     const float  Offset         = sum_par[4];
-    const float  PowCrosstlk    = sum_par[5];
+//    const float  PowCrosstlk    = sum_par[5];
 
 
@@ -733,5 +788,5 @@
     // create histo for height of Maximum amplitudes vs. pe
     double min_bin  =   hSum.GetBinCenter(maxbin);
-    double max_bin  =   hSum.GetBinCenter(maxOrder*(maxbin));
+    double max_bin  =   hSum.GetBinCenter((maxOrder-2)*(maxbin));
     double bin_width=   (max_bin - min_bin)/10;
 
@@ -741,9 +796,9 @@
     FitGaussiansToSpectrum
             (
-                maxOrder,
+                maxOrder-2,
                 hSum,
                 hMaxHeightsSum,
                 maxbin,
-                fwhm2,
+                fwhmSum,
                 gain
                 );
@@ -757,7 +812,7 @@
     TF1 fExpo1( "fExpo1","expo", min, max);
     fExpo1.SetLineColor(kRed);
-    hMaxHeightsSum.Fit(&fExpo1, "WLN0S" );
-    hMaxHeightsSum.DrawCopy("SAME");
-    fExpo1.DrawCopy("SAME");
+    hMaxHeightsSum.Fit(&fExpo1, "N0QS" );
+//    hMaxHeightsSum.DrawCopy("SAME");
+//    fExpo1.DrawCopy("SAME");
 
     //  EOF:    Fill and calculate sum spectrum
@@ -781,5 +836,5 @@
         hist->SetDirectory(0);
 
-        for(UInt_t bin = 1; bin < hSum.GetNbinsX()+1; bin++)
+        for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
         {
             hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1);
@@ -841,5 +896,5 @@
 
         //Fit Pixels spectrum
-        const TFitResultPtr rc = hist->Fit(&func, "WLN0QS", "", fit_min, fit_max);
+        const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max);
 
         const bool ok = int(rc)==0;
@@ -853,13 +908,4 @@
 //        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);
         hGain.Fill(fGain);
@@ -871,4 +917,5 @@
 
         hcam.SetBinContent(pixel+1, fGain);
+        hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);
 
         usePixel[pixel] = ok;
@@ -884,5 +931,5 @@
         {
 //            hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
-            hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
+            hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, (hist->GetBinContent(b)));
         }
 
@@ -895,4 +942,10 @@
             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);
@@ -911,17 +964,29 @@
 
     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 < hSum.GetNbinsX()+1; bin++)
+
+    for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
     {
         hSumScale.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
@@ -936,12 +1001,12 @@
 
     //Calculate full width half Maximum
-    fwhm2 = 0;
+    double fwhmScaled = 0;
     for (int i=1; i<maxbin; i++)
         if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
         {
-            fwhm2 = (2*i+1)*binwidth2;
+            fwhmScaled = (2*i+1)*binwidth2;
             break;
         }
-    if (fwhm2==0)
+    if (fwhmScaled==0)
     {
         cout << "Could not determine start value for sigma." << endl;
@@ -949,19 +1014,23 @@
 
     //Set fit parameters
-    Double_t par2[5] =
-    {
-        ampl2, 1, fwhm2*1.4, Crosstlk, 0
+    Double_t par2[6] =
+    {
+        ampl2, 0.9, fwhmScaled*1.4, Crosstlk, 0, 1.1
     };
 
-    func3.SetParameters(par2);
-    func3.FixParameter(1,1);
-    func3.FixParameter(4,0);
-
-    const double min2 = par2[1]-fwhm2*1.4;
-    const double max2 = par2[1]*maxOrder;
-//    func.SetRange(min2-fwhm2, max);
-    hSumScale.Fit(&func3, "WLN0QS", "", min2, max2);
-    func3.DrawCopy("SAME");
-
+    funcScaled.SetParameters(par2);
+//    funcScaled.FixParameter(1,0.9);
+    funcScaled.FixParameter(4,0);
+    funcScaled.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);
+    hSumScale.Fit(&funcScaled, "N0QS", "", minScaled, maxScaled);
+    funcScaled.DrawCopy("SAME");
     //--------fit gausses to peaks of gain corrected sum specetrum -----------
 
@@ -972,9 +1041,9 @@
     FitGaussiansToSpectrum
             (
-                maxOrder,
+                maxOrder-2,
                 hSumScale,
                 hMaxHeights,
                 maxbin2,
-                fwhm2,
+                fwhmScaled,
                 1
                 );
@@ -986,8 +1055,8 @@
     hMaxHeights.SetLineColor(kRed);
 //    hMaxHeights.DrawCopy("SAME");
-    TF1 fExpo( "fExpo","expo", 0,  maxOrder-1);
+    TF1 fExpo( "fExpo","expo", 0,  maxOrder-2);
     fExpo.SetLineColor(kRed);
     hMaxHeights.Fit(&fExpo, "N0QS" );
-    fExpo.DrawCopy("SAME");
+//    fExpo.DrawCopy("SAME");
 //    c12.cd();
 //    fExpo.DrawCopy("SAME");
@@ -1001,9 +1070,47 @@
     c2.cd(2);
     gPad->SetLogy();
-    TF1 GainGaus( "GainGaus", "gaus");
+//    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));
-    hGain.DrawCopy();
-    GainGaus.DrawCopy("SAME");
+//    GainGaus.DrawCopy("SAME");
 
     c2.cd(3);
@@ -1022,4 +1129,20 @@
     gPad->SetLogy();
     hFitProb.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");
+
 
     /*
@@ -1029,5 +1152,5 @@
     */
 
-    d->SaveAs(outfile);
+//    d->SaveAs(outfile);
     return 0;
 }
@@ -1073,5 +1196,5 @@
     Double_t cross  = par[3];
     Double_t shift  = par[4];
-    Double_t powOfX = par[5];
+//    Double_t powOfX = par[5];
 
     Double_t xTalk = 1;
@@ -1086,4 +1209,5 @@
         y += xTalk*gauss;
 
+//        xTalk = pow(cross, order*powOfX);
 //        for (int j = 1; j < powOfX; j++)
 //        {
@@ -1135,2 +1259,21 @@
 }
 
+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
+//----------------------------------------------------------------------------