Index: /fact/tools/marsmacros/singlepe_final.C
===================================================================
--- /fact/tools/marsmacros/singlepe_final.C	(revision 14164)
+++ /fact/tools/marsmacros/singlepe_final.C	(revision 14164)
@@ -0,0 +1,666 @@
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TSystem.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <TFitResultPtr.h>
+#include <TFitResult.h>
+
+#include "MH.h"
+#include "MLog.h"
+#include "MLogManip.h"
+#include "MDirIter.h"
+#include "MFillH.h"
+#include "MEvtLoop.h"
+#include "MCamEvent.h"
+#include "MGeomApply.h"
+#include "MTaskList.h"
+#include "MParList.h"
+#include "MContinue.h"
+#include "MBinning.h"
+#include "MHDrsCalib.h"
+#include "MDrsCalibApply.h"
+#include "MRawFitsRead.h"
+#include "MBadPixelsCam.h"
+#include "MStatusDisplay.h"
+#include "MTaskInteractive.h"
+#include "MPedestalSubtractedEvt.h"
+#include "MHCamera.h"
+#include "MGeomCamFACT.h"
+
+using namespace std;
+
+MPedestalSubtractedEvt *fEvent = 0;
+MLog *fLog = &gLog;
+
+struct Single
+{
+    float fSignal;
+    float fTime;
+};
+
+class MSingles : public MParContainer, public MCamEvent
+{
+    vector<vector<Single>> fData;
+
+public:
+    MSingles(const char *name=NULL, const char *title=NULL)
+    {
+        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(); }
+
+    Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
+    {
+    }
+    void   DrawPixelContent(Int_t num) const { }
+
+    ClassDef(MSingles, 1)
+};
+
+
+class MHSingles : public MH
+{
+    MSingles *fSingles;
+
+    TH2F fHist;
+
+public:
+    MHSingles() : fSingles(0)
+    {
+        fName = "MHSingles";
+
+        fHist.SetName("Name");
+        fHist.SetTitle("Title");
+        fHist.SetXTitle("X title");
+        fHist.SetYTitle("Y title");
+        fHist.SetDirectory(NULL);
+    }
+
+    Bool_t SetupFill(const MParList *plist)
+    {
+        fSingles = (MSingles*)plist->FindObject("MSingles");
+        if (!fSingles)
+        {
+            *fLog << /* err << */ "MSingles not found... abort." << endl;
+            return kFALSE;
+        }
+
+        return kTRUE;
+    }
+    Bool_t ReInit(MParList *)
+    {
+        MBinning binsx, binsy;
+        binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
+        binsy.SetEdges(300, -25, 600-25);
+
+        MH::SetBinning(fHist, binsx, binsy);
+
+        return kTRUE;
+    }
+
+    Int_t Fill(const MParContainer *par, const Stat_t weight=1)
+    {
+        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++)
+                fHist.Fill(i, result[j].fSignal);
+        }
+
+        return kTRUE;
+    }
+
+    TH2 *GetHist() { return &fHist; }
+
+    void Draw(Option_t *opt)
+    {
+        TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
+
+        AppendPad("");
+
+        pad->Divide(1,1);
+
+        pad->cd(1);
+
+        gPad->SetBorderMode(0);
+        gPad->SetFrameBorderMode(0);
+
+        fHist.Draw("colz");
+    }
+
+    ClassDef(MHSingles, 1)
+};
+
+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;
+    }
+
+    return kTRUE;
+}
+
+Int_t Process()
+{
+
+    for (int i=0; i<fEvent->GetNumPixels(); i++)
+    {
+        vector<Single> &result = fSingles->GetVector(i);
+        result.clear();
+
+        const UInt_t roi = fEvent->GetNumSamples();
+
+        const Float_t *ptr = fEvent->GetSamples(i);
+
+    //---------------------------
+        // ...Extract pulses here....
+    //---------------------------
+
+        //sliding average parameter (center weight = 3)
+    int avg1 = 8;
+        vector<float> Vslide(roi-2*avg1);//result of first sliding average
+
+        //Sliding window
+        for (UInt_t k=avg1; k<roi-avg1; k++)
+        {
+            Vslide[k-avg1]=ptr[k];
+            for(unsigned int offset=0; offset<avg1; offset++)
+            {
+                Vslide[k-avg1] += ptr[k+offset];
+                Vslide[k-avg1] += ptr[k-offset];
+            }
+            Vslide[k-avg1] /= float(2*avg1+1);
+        }
+
+    //peak finding via threshold
+        float threshold = 5;
+        Single single;
+        UInt_t start_sample = 5;
+        UInt_t end_sample = roi-2*avg1-21;
+        UInt_t min_dist = 20;
+        UInt_t integration_size = 5;
+        UInt_t integration_delay = 8;
+
+        for (UInt_t k=start_sample; k<end_sample; k++)
+        {
+      //search for threshold crossings
+            if(     (Vslide[k-start_sample]<threshold)
+                    &&(Vslide[k-1]<threshold)
+                    &&(Vslide[k]>threshold)
+                    &&(Vslide[k+5]>threshold)
+                    )
+            {
+                //average baseline before crossing
+                //double baseline = (Vslide[k-start_sample]+Vslide[k-start_sample+1]+Vslide[k-start_sample-1])/3.;
+                single.fSignal = 0;
+                single.fTime = k+integration_delay;
+                for(UInt_t l=integration_delay; l<integration_delay+integration_size; l++)
+                {
+                    single.fSignal+=Vslide[k+l];
+                }
+                //single.fSignal-=baseline;//baseline subtraction
+                result.push_back(single);
+                k+=min_dist;
+            }
+        }
+    }
+    return kTRUE;
+}
+
+Int_t PostProcess()
+{
+    return kTRUE;
+}
+
+Double_t
+spek_fit(
+        Double_t*   x,
+        Double_t*   par
+        );
+
+int singlepe_final()
+{
+    // ======================================================
+
+    const char *drsfile = "/fact/raw/2012/02/25/20120225_004.drs.fits.gz";
+
+    MDirIter iter;
+    iter.AddDirectory("/fact/raw/2012/06/01", "20120601_005.fits.gz");
+
+    // ======================================================
+
+    // true:  Display correctly mapped pixels in the camera displays
+    //        but the value-vs-index plot is in software/spiral indices
+    // false: Display pixels in hardware/linear indices,
+    //        but the order is the camera display is distorted
+    bool usemap = true;
+
+    // map file to use (get that from La Palma!)
+    const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
+    gStyle->SetPalette(1,0);
+    // ------------------------------------------------------
+
+    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;
+
+    gStyle->SetOptFit(kTRUE);
+
+    // ======================================================
+
+    gLog << endl;
+    gLog.Separator("Processing external light pulser run");
+
+    MTaskList tlist1;
+
+    MSingles singles;
+
+    MParList plist1;
+    plist1.AddToList(&tlist1);
+    plist1.AddToList(&drscalib300);
+    plist1.AddToList(&badpixels);
+    plist1.AddToList(&singles);
+
+    MEvtLoop loop1("DetermineCalConst");
+    loop1.SetDisplay(d);
+    loop1.SetParList(&plist1);
+
+    // ------------------ Setup the tasks ---------------
+
+    MRawFitsRead read1;
+    read1.LoadMap(map);
+    read1.AddFiles(iter);
+
+    MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
+
+    MGeomApply apply1;
+
+    MDrsCalibApply drsapply1;
+
+    MTaskInteractive mytask;
+    mytask.SetPreProcess(PreProcess);
+    mytask.SetProcess(Process);
+    mytask.SetPostProcess(PostProcess);
+
+    MFillH fill("MHSingles");
+
+    // ------------------ Setup eventloop and run analysis ---------------
+
+    tlist1.AddToList(&read1);
+//    tlist1.AddToList(&cont1);
+    tlist1.AddToList(&apply1);
+    tlist1.AddToList(&drsapply1);
+    tlist1.AddToList(&mytask);
+    tlist1.AddToList(&fill);
+
+    if (!loop1.Eventloop(max1))
+        return 10;
+
+    if (!loop1.GetDisplay())
+        return 11;
+
+    MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
+    if (!h)
+        return 0;
+
+    TH2 *hist = h->GetHist();
+
+    // ------------------ Spectrum Fit ---------------
+    // AFTER this the Code for the spektrum fit follows
+    // access the spektrumhistogram by the pointer *hist
+
+    hist->SetTitle("Amplitude Spectrum vs. Pixel ID");
+    hist->SetXTitle("Pixel SoftID");
+    hist->SetYTitle("Amplitude [mV]");
+    hist->SetZTitle("counts");
+
+    MGeomCamFACT fact;
+    MHCamera hcam(fact);
+
+    struct fit
+    {
+        float fAmplitude;
+        float fGain;
+        float f2GainRMS;
+        float fCrosstlk;
+        float fOffset;
+    };
+    fit fitParameters[singles.GetNumPixels()];
+
+    TCanvas &canv1 = d->AddTab("Fits");
+    TCanvas &canv2 = d->AddTab("Distributions");
+    canv2.Divide(3,2);
+    TCanvas &gain_canv = d->AddTab("Gain Distribution");
+
+    //built an array of Amplitude spectra
+    TH1D*       hAmplSpek[singles.GetNumPixels()];
+    TString     histotitle;
+
+    TF1 gaus4("g4","gaus");
+    gaus4.SetLineColor(2);
+
+    TH1F        hAmplitude(
+                    "hAmplitude",
+                    "Amplitude of Single Preak",
+                    200,
+                    50,
+                    250
+                    );
+    hAmplitude.GetXaxis()->SetTitle("Amplitude (from Integral 5 Slices) [mV]");
+    hAmplitude.GetYaxis()->SetTitle("counts");
+
+    TH1F        hGain(
+                    "hGain",
+                    "Gain",
+                    40*5,
+                    30,
+                    70
+                    );
+    hGain.GetXaxis()->SetTitle("Gain (from Integral 5 Slices) [mV]");
+    hGain.GetYaxis()->SetTitle("counts");
+
+    TH1F        hGainRMS(
+                    "hSigma",
+                    "Sigma of gauss peak in gainfit",
+                    20*5,
+                    0,
+                    40
+                    );
+    hGainRMS.GetXaxis()->SetTitle("Sigma (from Integral 5 Slices) [mV]");
+    hGainRMS.GetYaxis()->SetTitle("counts");
+
+    TH1F        hCrosstlk(
+                    "hCrosstlk",
+                    "Parameter proportional to Crosstalk",
+                    30*5,
+                    0,
+                    4
+                    );
+    hCrosstlk.GetXaxis()->SetTitle("~Crosstalk (from Integral 5 Slices) [mV]");
+    hCrosstlk.GetYaxis()->SetTitle("counts");
+
+    TH1F        hOffset(
+                    "hOffset",
+                    "X-Offset of Spectrum",
+                    100*2,
+                    0,
+                    4
+                    );
+    hOffset.GetXaxis()->SetTitle(" X-Offset of Spectrum (from Integral 5 Slices) [mV]");
+    hOffset.GetYaxis()->SetTitle("counts");
+
+    TH1F        hFitProb(
+                    "hFitProb",
+                    "Fit Probability",
+                    100*2,
+                    0,
+                    100
+                    );
+    hFitProb.GetXaxis()->SetTitle("fit-probability [%]");
+    hFitProb.GetYaxis()->SetTitle("counts");
+
+    // define fit function for Amplitudespectrum
+    TF1         fspek_fit("fspek_fit", spek_fit, 10, 600, 5);
+    fspek_fit.SetParNames("Ampl 1.Max", "Gain", "Sigma", "Ampl 2.Max", "Offset");
+
+    // fit parameters
+    Double_t    par[5];
+
+    // Map for which pixel shall be plotted and which not
+    TArrayC usePixel(1440);
+
+    // List of Pixel that should be ignored in camera view
+    Int_t doNotUse[12]= {
+        424,
+        923,
+        1208,
+        583,
+        830,
+        1399,
+        113,
+        115,
+        354,
+        423,
+        1195,
+        1393
+    };
+
+    // Begin of Loop over Pixels
+    for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
+    {
+        usePixel[pixel] = 0;
+
+        //define projection histogram title
+        histotitle = "hPixelPro_";
+        histotitle += pixel;
+
+        //Projectipon of a certain Pixel out of Ampl.Spectrum
+        hAmplSpek[pixel] = hist->ProjectionY( histotitle , pixel+1, pixel+1);
+
+        //set histogram and Axis title
+        histotitle = "Amplitude Spectrum of Pixel #";
+        histotitle += pixel;
+        hAmplSpek[pixel]->SetTitle(histotitle);
+
+        // define appearance of fit
+        fspek_fit.SetLineColor(kBlue);
+        fspek_fit.SetLineStyle(2);
+
+        par[0] = hAmplSpek[pixel]->GetBinContent(hAmplSpek[pixel]->GetMaximumBin());
+        // par[1] + par[4] is the position of the first maximum
+        par[1] = 0.8*hAmplSpek[pixel]->GetBinCenter(hAmplSpek[pixel]->GetMaximumBin());
+        par[2] = 0.2*par[1];
+        par[3] = 0.1*par[0];
+        par[4] = 0.2*hAmplSpek[pixel]->GetBinCenter(hAmplSpek[pixel]->GetMaximumBin());
+        fspek_fit.SetParameters(par);
+
+        //Fit Pixels spectrum
+        canv1.cd();
+        TFitResultPtr fitStatus(hAmplSpek[ pixel ]->Fit(&fspek_fit,"+S","",0.8*par[1],5*par[1]));
+        double fit_prob = fitStatus->Prob();
+
+        //cout << "pixel: " << pixel << "\t status: " << fitStatus << endl;
+
+        fitParameters[pixel].fAmplitude = fspek_fit.GetParameter(0);
+        fitParameters[pixel].fGain      = fspek_fit.GetParameter(1);
+        fitParameters[pixel].f2GainRMS  = fspek_fit.GetParameter(2);
+        fitParameters[pixel].fCrosstlk  = fspek_fit.GetParameter(3)/fspek_fit.GetParameter(0);
+        fitParameters[pixel].fOffset    = fspek_fit.GetParameter(4);
+
+        gain_canv.cd();
+        gPad->SetLogy(1);
+        hGain.Fill(fitParameters[pixel].fGain);
+        hGain.Fit(&gaus4);
+        hGain.DrawCopy();
+
+        canv2.cd(2);
+        gPad->SetLogy(1);
+        hGain.DrawCopy();
+
+        canv2.cd(3);
+        gPad->SetLogy(1);
+        hGainRMS.Fill(fitParameters[pixel].f2GainRMS);
+        hGainRMS.DrawCopy();
+
+        canv2.cd(1);
+        gPad->SetLogy(1);
+        hAmplitude.Fill(fitParameters[pixel].fAmplitude);
+        hAmplitude.DrawCopy();
+
+        canv2.cd(4);
+        gPad->SetLogy(1);
+        hCrosstlk.Fill(fitParameters[pixel].fCrosstlk);
+        hCrosstlk.DrawCopy();
+
+        canv2.cd(5);
+        gPad->SetLogy(1);
+        hOffset.Fill(fitParameters[pixel].fOffset);
+        hOffset.DrawCopy();
+
+        canv2.cd(6);
+        gPad->SetLogy(1);
+        hFitProb.Fill(100*fit_prob);
+        hFitProb.DrawCopy();
+
+
+        hcam.SetBinContent(pixel+1, fitParameters[pixel].fGain);
+        if (!fitStatus)
+        {
+            usePixel[pixel] = 1;
+            cout << "Pixel: " << pixel << " success!" << endl << endl;
+        }
+        else
+        {
+            cout << "Pixel: " << pixel << " failed!" << endl << endl;
+        }
+
+    }
+//    TCanvas &canvBroken[12];
+//    TString name;
+    for  (int i = 0; i < 12; i++)
+    {
+//        name    = "Pixel";
+//        name   += doNotUse[i];
+//        canvBroken[i] = d->AddTab(name);
+//        hAmplSpek[doNotUse[ i ]]->Draw();
+
+        usePixel[doNotUse[ i ]] = 0;
+    }
+
+    hcam.SetUsed(usePixel);
+
+    TCanvas &canv = d->AddTab("Gain");
+    canv.cd();
+    hcam.SetNameTitle( "gain","Gain distribution over whole camera");
+    hcam.DrawCopy();
+
+
+    canv1.cd();
+
+    /*
+    TString title = "--  Calibrated signal #";
+    title += seq.GetSequence();
+    title += " (";
+    title += drsfile;
+    title += ")  --";
+    d->SetTitle(title, kFALSE);
+
+    TString path;
+    path += Form("%s/20%6d_%03d-calibration.root", outpath,
+                 seq.GetSequence()/1000, seq.GetSequence()%1000);
+
+    d->SaveAs(path);
+    */
+
+    return 0;
+}
+
+Double_t
+spek_fit(
+        Double_t*   x,
+        Double_t*   par
+        )
+{
+    Double_t arg        = 0;
+    Double_t cross      = 1;
+    Double_t peak       = 0;
+//    Double_t arg2     = 0;
+//    Double_t arg3     = 0;
+//    Double_t arg4     = 0;
+
+//  if (par[0] != 0) cross   = par[3]/par[0];
+//  if (par[2] != 0) arg     = (x[0] - par[1] + par[4])/par[2];
+//  if (par[2] != 0) arg2    = (x[0] - 2*par[1] + par[4])/par[2];
+//  if (par[2] != 0) arg3    = (x[0] - 3*par[1] + par[4])/par[2];
+//  if (par[2] != 0) arg4    = (x[0] - 4*par[1] + par[4])/par[2];
+//    if (par[0] = 0) return 0;
+//    if (par[2] = 0) return 0;
+
+    if (par[0] != 0) cross   = par[3]/par[0];
+    for (int order = 1; order < 5; order++)
+    {
+        arg = (x[0] - order*par[1] - par[4] )/par[2];
+        Double_t xTalk    = TMath::Power(cross, order-1);
+//        Double_t gauss    = TMath::Exp(-0.5 * arg * arg/order);
+        Double_t gauss    = TMath::Exp(-0.5 * arg * arg);
+        peak += xTalk*par[0]*gauss;
+    }
+
+//  peak += par[0]*TMath::Exp(-0.5 * arg * arg);
+//  peak += cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2);
+//  peak += cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3);
+//  peak += cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4);
+
+    return peak;
+}