Changeset 12392 for fact/tools/rootmacros

Timestamp:

11/05/11 13:32:59 (13 years ago)

Author:

neise

Message:

tidied up tpeak.C and put an entry into README.txt

Location:

fact/tools/rootmacros

Files:

• README.txt (modified) (1 diff)
• tpeak.C (modified) (20 diffs)

fact/tools/rootmacros/README.txt

@@ -105 +105 @@
 
 
+-----------------------------------------------------------------------------------------------
+tpeak.C ROOT macros to produce overlays of data around a found peak in a given
+window.
+
+int tpeak(
+  char *datafilename    = "data/20111016_013.fits.gz",
+  const char *drsfilename = "../../20111016_011.drs.fits.gz",
+  const char *OutRootFileName = "../analysis/fpeak_cdf.Coutput.root",
+  int firstevent      = 0,
+  int nevents       = -1,
+  int firstpixel      = 0,
+  int npixel        = -1,
+  bool spikeDebug = false,
+  int avg1    = 14,
+  int avg2    = 8,
+  int OverlayWindowLeft = 50,
+  int OverlayWindowRight = 150,
+  int verbosityLevel = 1, // different verbosity levels can be implemented
+here
+  bool ProduceGraphic = true
+  )
+
+call it like this, if you want to overlay a certain number of peaks for a
+single channel
+tpeak("data/20111029_017.fits.gz","data/20111029_013.drs.fits.gz","test.root",0,1000,333,1,true,0,0,50,150,1,true)
+
+set the first bool 'spikeDebug' to false in order to overlay quicker ... every
+50th event the canvas is updated.
+depending on what kind of peaks you like to detect the sliding average filters
+should be set to e.g. 14,8 for singles in a quiet G-APD pedestal run
+or 0,0 if you want to see just every thing that might be a signal...
+
+the window size 50,150 means .. 50 slices to the left of the maximum and 150
+to the right.
+
+there is a bit of cleaning included in the file, maybe you want to switch it
+off, then just search for the calls of these methods... 
+ removeMaximaBelow( *zXings, 3.0);
+ removeRegionWithMaxOnEdge( *zXings, 2);
+ removeRegionOnFallingEdge( *zXings, 100);
+
+and comment them out or play with the settings.
+
+they are defined in zerosearch.C, which is maybe a bad choice...
+
+
+

fact/tools/rootmacros/tpeak.C

@@ -42 +42 @@
 #include "factfir.C"
 
+#include "FOpenDataFile.h"
+#include "FOpenDataFile.c"
+
+#include "DrsCalibration.C"
+#include "DrsCalibration.h"
+
+#include "SpikeRemoval.h"
+#include "SpikeRemoval.C"
 
 //-----------------------------------------------------------------------------
 // Decleration of variables
 //-----------------------------------------------------------------------------
-
+bool breakout=false;
+
+int NEvents;
 vector<int16_t> AllPixelDataVector;
 vector<int16_t> StartCellVector;
-
 unsigned int CurrentEventID;
-
-bool breakout=false;
-
-size_t ROIxNOP;
+size_t PXLxROI;
+UInt_t RegionOfInterest;
 UInt_t NumberOfPixels;
-UInt_t RegionOfInterest;
-int NEvents;
 
 size_t drs_n;
@@ -64 +69 @@
 vector<float> drs_triggeroffsetmean;
 
-int FOpenDataFile( fits & );
-
-
 vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
-vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
 vector<float> Vcorr(FAD_MAX_SAMPLES);  // corrected Values
-vector<float> Vdiff(FAD_MAX_SAMPLES);  // numerical derivative
-
 vector<float> Vslide(FAD_MAX_SAMPLES);  // sliding average result
 vector<float> Vcfd(FAD_MAX_SAMPLES);    // CDF result
 vector<float> Vcfd2(FAD_MAX_SAMPLES);    // CDF result + 2nd sliding average
-
-
-float getValue( int, int );
-float correctDrsOffset( int slice, int pixel );
-void computeN1mean( int );
-void removeSpikes( int );
 
 // histograms
@@ -95 +88 @@
 TH1F* h;
 TH1F *debugHistos;
-TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
-                    // x = pixel id, y = DRS cell id
-
-TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
-TH1F *hMeanBsl, *hpltMeanBsl;
-TH1F *hRmsBsl, *hpltRmsBsl;
-TH2F * hAmplSpek_cfd;
-TH2F * hAmplSpek_discri;   
-
 TH2F *hPeakOverlay;  //histogrammm for overlay of detected Peaks
+int hPeakOverlayXaxisLeft,hPeakOverlayXaxisRight; 
 
 TObjArray hList;
@@ -110 +95 @@
 
 void BookHistos( );
-void SaveHistograms( char * );
-
-int searchSinglesPeaks = 0;
+void SaveHistograms( const char * );
 
 //-----------------------------------------------------------------------------
 // Main
 //-----------------------------------------------------------------------------
-
 int tpeak(
-        char *datafilename      = "../../20111011_055.fits.gz",
-        const char *drsfilename = "../../20111011_054.drs.fits.gz",
-        int PixelID             = -1,
-        int firstevent          = 0,
-        int nevents             = -1,
-        bool spikeDebug         = false,
-        int verbosityLevel      = 1 // different verbosity levels can be implemented here
-        )
+  char *datafilename    = "data/20111016_013.fits.gz",
+  const char *drsfilename = "../../20111016_011.drs.fits.gz",
+  const char *OutRootFileName = "../analysis/fpeak_cdf.Coutput.root",
+  int firstevent      = 0,
+  int nevents       = -1,
+  int firstpixel      = 0,
+  int npixel        = -1,
+  bool spikeDebug = false,
+  int avg1    = 14,
+  int avg2    = 8,
+        int OverlayWindowLeft = 50,
+        int OverlayWindowRight = 150,
+  int verbosityLevel = 1, // different verbosity levels can be implemented here
+  bool ProduceGraphic = true
+  )
+
 {
+hPeakOverlayXaxisLeft = OverlayWindowLeft;
+hPeakOverlayXaxisRight = OverlayWindowRight;
+
+        gStyle->SetPalette(1,0);
+        gROOT->SetStyle("Plain");
 
 //-----------------------------------------------------------------------------
@@ -133 +128 @@
 // also in 'non-debug' runs
 //-----------------------------------------------------------------------------
-
-        TCanvas * cSpektrum;
-        TCanvas * cStartCell;
-        cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400);
-        cSpektrum->Divide(1,2);
-        cStartCell = new TCanvas("cStartCell ", "The Startcells of this run", 10,410,400,400);
-
         // Canvases only need if spike Debug, but I want to deklare
         // the pointers anyway ...
-        TCanvas *cRawAndSpikeRemoval = NULL;
         TCanvas *cFiltered = NULL;
         if (spikeDebug){
-                cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",410,10,400,400);
-                cRawAndSpikeRemoval->Divide(1, 2);
-                cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
-                cFiltered->Divide(1, 2);
+                cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",1,10,400,300);
+                cFiltered->Divide(1, 3);
         }
         // Canvases to show the peak template
         TCanvas *cPeakOverlay = NULL;
-        cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 410, 10, 600, 600);
-
-gStyle->SetPalette(1,0);
-gROOT->SetStyle("Plain");
+        cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 1, 310, 400, 400);
+
 
 //-----------------------------------------------------------------------------
@@ -169 +152 @@
 // Filter-Settings
 //-----------------------------------------------------------------------------
-
-// sliding window filter settings
-        int k_slide = 16;
-        vector<double> a_slide(k_slide, 1); //erzeugt vector der länge k_slide, und füllt die elemente mit 1
-        double b_slide = k_slide;
-
 // CFD filter settings
         int k_cfd = 10;
@@ -182 +159 @@
         a_cfd[k_cfd-1]=1.;
 
-// 2nd slinding window filter
-        int ks2 = 16;
-        vector<double> as2(ks2, 1);
-        double bs2 = ks2;
-
 //-----------------------------------------------------------------------------
 // prepare datafile
@@ -192 +164 @@
 
 // Open the data file
-        fits *datafile = new fits( datafilename );
-        if (!datafile) {
-            printf( "Could not open the file: %s\n", datafilename );
-                  return 1;
-        }
-
-// access data
-        NEvents = FOpenDataFile( *datafile );
-        printf("number of events in file: %d\n", NEvents);
+
+  fits * datafile;
+  // Opens the raw data file and 'binds' the variables given as
+  // Parameters to the data file. So they are filled with
+  // raw data as soon as datafile->GetRow(int) is called.
+  NEvents = OpenDataFile( datafilename, &datafile,
+    AllPixelDataVector, StartCellVector, CurrentEventID,
+    RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
+  if (NEvents == 0){
+    cout << "return code of OpenDataFile:" << datafilename<< endl;
+    cout << "is zero -> aborting." << endl;
+    return 1;
+  }
+
+  if (verbosityLevel > 0)
+    cout <<"number of events in file: "<< NEvents << "\t";
   if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
+  if (verbosityLevel > 0)
+    cout <<"of, which "<< nevents << "will be processed"<< endl;
+
+  if (verbosityLevel > 0)
+    cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
+  if ( npixel == -1 || npixel > (int)NumberOfPixels  ) npixel = (int)NumberOfPixels; // -1 means all!
+  if (verbosityLevel > 0)
+    cout <<"of, which "<< npixel << "will be processed"<< endl;
 
 //Get the DRS calibration
@@ -210 +197 @@
                                         drs_n);
 
-//Check the sizes of the data columns
-        if (drs_n != 1474560) {
-                cerr << "error: DRS calib file has wrong ...erm...size ... drs_n is: "
-                        << drs_n << endl;
-                cerr << " Aborting." << endl;
-                return 1;
-        }
-
-        if(ROIxNOP != 1474560)
-        {
-                cout << "warning: data_n should better be 1440x1024=1474560, but is "
-                        << ROIxNOP << endl;
-                cout << "this script is not guaranteed to run under these "
-                        <<" circumstances....any way ... it is never guaranteed." << endl;
-        }
-
 // Book the histograms
-
         BookHistos( );
 
-        float myTHR = 3.5;
-        float calibratedVoltage;
 //-----------------------------------------------------------------------------
 // Loops over Every Event and Pixel
@@ -243 +211 @@
         //-------------------------------------
 
-                for ( int pix = 0; pix < 1440; pix++ ){
-
-                        hStartCell->Fill( pix, StartCellVector[pix] );
-
-                        // this is a stupid hack ... there is more code at the
-                        // end of this loop to complete this hack ...
-                        // beginning with if (PixelID != -1) as well.
-                        if (PixelID != -1) {
-                                pix = PixelID;
-                                if (verbosityLevel > 0){
-                                        cout << "Processing Event number: " << CurrentEventID << "\t"
-                                                << "Pixel number: "<< pix << endl;
-                                }
-                        }
-
-                        if (verbosityLevel > 0){
-                                if (pix % 20 ==0){
-                                        cout << "Processing Event number: " << CurrentEventID << "\t"
-                                                << "Pixel number: "<< pix << endl;
-                                }
-                        }
-
-                        // compute the DRs calibrated values and put them into Ameas[]
-                        for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
-                                calibratedVoltage = getValue( sl, pix);
-                                if (verbosityLevel > 10){
-                                        printf("calibratedVoltage = %f\n", calibratedVoltage);
-                                }
-                                Ameas[ sl ] = calibratedVoltage;
-
-                                // in case we want to plot this ... we need to put it into a
-                                // Histgramm
-                                if (spikeDebug){
-                                        debugHistos[Ameas_].SetBinContent(sl, calibratedVoltage);
-                                }
-                        }
-                        // operates on Ameas[] and writes to N1mean[]
-                        computeN1mean( RegionOfInterest );
-
-                        // operates on Ameas[] and N1mean[], then writes to Vcorr[]
-                        removeSpikes( RegionOfInterest );
-
-                        // filter Vcorr with sliding average using FIR filter function
-                        factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
-                        // filter Vslide with CFD using FIR filter function
-                        factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
-                        // filter Vcfd with sliding average using FIR filter function
-                        factfir(bs2 , as2, ks2, Vcfd, Vcfd2);
-
-
-                        vector<Region> *Regions =  discriminator( Vslide, myTHR,true , 120);
-                        for (unsigned int p=0; p<Regions->size(); p++ ){
-                                        hAmplSpek_discri->Fill(pix , Regions->at(p).maxVal);
-                        }
+                for ( int pix = firstpixel; pix < firstpixel + npixel; pix++ ){
+
+                                                                                                if (verbosityLevel > 0){
+                                                                                if (pix == firstpixel){
+                                                                                  cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl;
+                                                                                }
+                                                                              }
+
+                                                                                        applyDrsCalibration( Ameas,pix,12,12,
+                                                                        drs_basemean, drs_gainmean, drs_triggeroffsetmean,
+                                                                        RegionOfInterest, AllPixelDataVector, StartCellVector);
+
+                                                                              // finds spikes in the raw data, and interpolates the value
+                                                                              // spikes are: 1 or 2 slice wide, positive non physical artifacts
+                                                                              removeSpikes (Ameas, Vcorr);
+
+                                                                                     // filter Vcorr with sliding average using FIR filter function
+                                                                              sliding_avg(Vcorr, Vslide, avg1);
+                                                                              // filter Vslide with CFD using FIR filter function
+                                                                              factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
+                                                                              // filter Vcfd with sliding average using FIR filter function
+                                                                              sliding_avg(Vcfd, Vcfd2, avg2);
 
                         // peaks in Ameas[] are found by searching for zero crossings
@@ -304 +240 @@
                         vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 1);
                         // zXings means "zero cross ings"
-                        ShiftRegionBy(*zXings, -ks2/2);
                         EnlargeRegion(*zXings, 10, 10);
                         findAbsMaxInRegions(*zXings, Vslide);
-
-                        for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
-                           // hPixelCellData.Fill( sl, Vcorr[sl] );
-                           hBaseline[pix]->Fill( Vslide[sl] ) ;
-                        }
+                                                                              removeMaximaBelow( *zXings, 3.0);
+                                                                              removeRegionWithMaxOnEdge( *zXings, 2);
+                                                                              removeRegionOnFallingEdge( *zXings, 100);
 
                         //following Code should produce the Overlay of peaks
@@ -317 +250 @@
 
                                 for (it = zXings->begin() ; it < zXings->end() ; it++){
-                                  int a = it->maxPos-100;
-                                  int b = it->maxPos+100;
-                                  if (a < 0)
-                                    a =0;
-                                  if (b>1023)
-                                    b=1023;
-                                        for ( int pos = a; pos <= b; pos++){
+                                  int Left = it->maxPos - OverlayWindowLeft;
+                                  int Right = it->maxPos + OverlayWindowRight;
+                                  if (Left < 0)
+                                    Left =0;
+                                  if (Right > (int)Vcorr.size() )
+                                    Right=Vcorr.size() ;
+                                        for ( int pos = Left; pos < Right; pos++){
 
                                             hPeakOverlay->Fill( pos - (it->maxPos), Vcorr[pos]) ;
@@ -330 +263 @@
 
 
-
-                        if (zXings->size() != 0 ){
-                                for (unsigned int i=0; i<zXings->size(); i++){
-                                        if (verbosityLevel > 1){
-                                                cout << zXings->at(i).maxPos << ":\t"
-                                                        << zXings->at(i).maxVal <<endl;
-                                        }
-                                        hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);
-                                }
-                        }
-
                         if ( spikeDebug ){
-                                for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
-                                        debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
-                                        debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
-                                        debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
-                                }
-
-                                cRawAndSpikeRemoval->cd( 1);
-                                gPad->SetGrid();
-                                debugHistos[Ameas_].Draw();
-
-                                cRawAndSpikeRemoval->cd( 2);
+
+
+                                                                                                          // TODO do this correct. The vectors should be the rigt ones... this is just luck 
+                                                                                  debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
+                                                                                  debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
+                                                                                  debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
+                                                                                  debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
+                                                                                  debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
+
+                                                                                        for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
+                                                                                          debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
+                                                                                          debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
+                                                                                          debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
+                                                                                          debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
+                                                                                          debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
+                                                                                        }
+
+
+                                cFiltered->cd(1);
                                 gPad->SetGrid();
                                 debugHistos[Vcorr_].Draw();
 
-                                cFiltered->cd(1);
+                                cFiltered->cd(2);
                                 gPad->SetGrid();
                                 debugHistos[Vslide_].Draw();
-                                TLine * thrLine = new TLine(0, myTHR, 1024, myTHR);
-                                thrLine->SetLineColor(kRed);
-                                thrLine->SetLineWidth(2);
-                                thrLine->Draw();
-                                TLine * OneLine;
-                                vector<TLine*> MyLines;
-                                for (unsigned int p=0; p<Regions->size(); p++ ){
-                                        OneLine = new TLine(
-                                                Regions->at(p).begin ,
-                                                Regions->at(p).maxVal,
-                                                Regions->at(p).end,
-                                                Regions->at(p).maxVal);
-                                        OneLine->SetLineColor(kRed);
-                                        OneLine->SetLineWidth(2);
-                                        MyLines.push_back(OneLine);
-                                        OneLine->Draw();
-                                }
+
                                 TBox *OneBox;
                                 vector<TBox*> MyBoxes;
@@ -392 +306 @@
                                 }
 
-                                cFiltered->cd(2);
+                                cFiltered->cd(3);
                                 gPad->SetGrid();
                                 debugHistos[Vcfd2_].Draw();
@@ -399 +313 @@
                                 zeroline->Draw();
 
-                                cRawAndSpikeRemoval->Update();
                                 cFiltered->Update();
 
@@ -423 +336 @@
                         }// end of if(spikeDebug)
 
-                        delete Regions;
                         delete zXings;
-
-                        // this is the 2nd part of the ugly hack in the beginning.
-                        // this makes sure, that the loop ends
-                        if (PixelID != -1){
-                                pix = 2000;
-                        }
-
                         if (breakout)
                                 break;
@@ -439 +344 @@
         // end of loop over pixels
         //-------------------------------------
-if (ev % 50 == 0){
-                cSpektrum->cd(1);
-                hAmplSpek_cfd->Draw("COLZ");
-                cSpektrum->cd(2);
-                hAmplSpek_discri->Draw("COLZ");
-                cSpektrum->Modified();
-                cSpektrum->Update();
-
-                // updating seems not to work ..
-                // debug cout
-                cStartCell->cd();
-                hStartCell->Draw();
-                cStartCell->Modified();
-                cStartCell->Update();
-
+                                if (ProduceGraphic){
+                                                if (ev % 50 == 0){
                 //OVERLAY PEAKS
                 cPeakOverlay->cd();
@@ -459 +351 @@
                 cPeakOverlay->Modified();
                 cPeakOverlay->Update();
-
-}
-
-
+                                                }
+                                }
                 if (breakout)
                         break;
@@ -469 +359 @@
 // end of loop over Events
 //-------------------------------------
-
-
-        for (int pix = 0; pix < NumberOfPixels; pix++) {
-            //printf("Maximum bin pix[%d] %f\n", pix ,hBaseline[pix]->GetMaximumBin() );
-
-            float vmaxprob = hBaseline[pix]->GetXaxis()->GetBinCenter(
-                hBaseline[pix]->GetMaximumBin() );
-
-            printf("%8.3f", vmaxprob );
-
-            hMeanBsl->Fill( vmaxprob );
-            hpltMeanBsl->SetBinContent(pix+1, vmaxprob );
-
-            hRmsBsl->Fill(hBaseline[pix]->GetRMS() );
-            hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() );
-        }
-
-
-        SaveHistograms( datafilename );
+                                if (ProduceGraphic){
+                //OVERLAY PEAKS
+                cPeakOverlay->cd();
+                hPeakOverlay->Draw("COLZ");
+                cPeakOverlay->Modified();
+                cPeakOverlay->Update();
+                                }
+
+        SaveHistograms( OutRootFileName );
 
         return( 0 );
 }
 
-void removeSpikes(int Samples){
-
-    const float fract = 0.8;
-    float x, xp, xpp, x3p;
-
-    // assume that there are no spikes
-    for ( int i = 0; i <  Samples; i++) Vcorr[i] = Ameas[i];
-
-// find the spike and replace it by mean value of neighbors
-    for ( int i = 0; i < Samples; i++) {
-
-    // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ i ] );
-
-    x = Ameas[i] - N1mean[i];
-
-        if ( x < -5. ){ // a spike candidate
-            // check consistency with a single channel spike
-            xp = Ameas[i+1] - N1mean[i+1];
-            xpp = Ameas[i+2] - N1mean[i+2];
-            x3p = Ameas[i+3] - N1mean[i+3];
-
-            // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p);
-
-            if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
-                // printf("double spike candidate\n");
-                Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
-                Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
-                // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]);
-                // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]);
-                i = i + 3;
-            }
-            else{
-
-                if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
-                    Vcorr[i+1] = N1mean[i+1];
-                    // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
-                    // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
-                    N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
-                    i = i + 2;//do not care about the next sample it was the spike
-                }
-                // treatment for the end of the pipeline must be added !!!
-            }
-        }
-        else{
-             // do nothing
-        }
-    } // end of spike search and correction
-        for ( int i = 0; i < Samples; i++ ) debugHistos[ Vcorr_ ].SetBinContent( i, Vcorr[i] );
-}
-/*
-void computeN1mean( int Samples ){
-cout << "In compute N1mean" << endl;
-// compute the mean of the left and right neighbors of a channel
-
-    for( int i = 0; i < Samples; i++){
-        if (i == 0){ // use right sample es mean
-            N1mean[i] = Ameas[i+1];
-        }
-        else if ( i == Samples-1 ){ //use left sample as mean
-            N1mean[i] = Ameas[i-1];
-        }
-        else{
-            N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
-        }
-        h[tN1mean].SetBinContent(i, Ameas[i] - N1mean[i]);
-    }
-} // end of computeN1mean computation
-*/
-
-void computeN1mean( int Samples ){
-// compute the mean of the left and right neighbors of a channel
-
-    for( int i = 2; i < Samples - 2; i++){
-/*        if (i == 0){ // use right sample es mean
-            N1mean[i] = Ameas[i+1];
-        }
-        else if ( i == Samples-1 ){ //use left sample as mean
-            N1mean[i] = Ameas[i-1];
-        }
-        else{
-            N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
-        }
-*/
-        N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
-    }
-} // end of computeN1mean computation
-
-
-
-/*
-// this function takes a vector<float> and tries to apply the DRS calibration to it.
-// therfor it uses a DRS calibration file, which must be given as a parameter.
-// in case the calibration file can not be opened the function generates a warning on cerr
-// and returns an empty vector
-//
-// the function assumes the first parameter conatins the data of all Pixel
-// the RegionOfInterest is deduced from its size().
-
-        //
-
-vector<float> * calibrateDrsData(
-        vector<float> & AllPixelData,
-        vector<float> & AllStartCells,
-        const char * DrsCalibFilePath,
-        bool ApplyCalibToSourceVector = false)
-{
-        vector<float> * CalibratedDRSData;
-
-        size_t DrsLengthOfVectors;
-        vector<float> DrsOffset;
-        vector<float> DrsGain;
-        vector<float> DrsOffsetSpecialTrigger;
-
-        // the Total Number of Pixel, we have data of
-        // is deduced from the size() of AllPixelStartCell
-        unsigned int NumberOfPixel = AllStartCells.size();
-        //sanity Check:
-        if (NumberOfPixel < 1) {
-                cerr << "calibrateDrsData - Error: NumberOfPixel=AllStartCells.size() ia:"
-                        << AllStartCells.size() << " must be >= 1"  << endl;
-                return CalibratedDRSData
-        }
-
-        // The RegionOfInterest is deduced from the size() of AllPixel Data
-        unsigned int RegionOfInterest = AllPixelData.size() / NumberOfPixel;
-        // sanity Check
-        if ( RegionOfInterest < 1){
-                cerr << "calibrateDrsData - Error: RegionOfInterest = AllPixelData.size() / NumberOfPixel is:"
-                        << RegionOfInterest << " must be >= 1"  << endl;
-                return CalibratedDRSData
-        }
-
-        // now lets see if we can find the drs calib file
-        FOpenCalibFile( DrsCalibFilePath,
-                DrsOffset,
-                DrsGain,
-                DrsOffsetSpecialTrigger,
-                DrsLengthOfVectors);
-        //sanity check is done in FOpenCalibFile, I guess ... check and maybe TODO
-        // maybe at least the return code of FOpenCalibFile should be checked here ... TODO
-
-
-
-
-}
-
-*/
-
-float getValue( int slice, int pixel ){
-        const float dconv = 2000/4096.0;
-
-        float vraw, vcal;
-
-        unsigned int pixel_pt;
-        unsigned int slice_pt;
-        unsigned int cal_pt;
-        unsigned int drs_cal_offset;
-
-        // printf("pixel = %d, slice = %d\n", slice, pixel);
-
-        pixel_pt = pixel * RegionOfInterest;
-        slice_pt = pixel_pt + slice;
-        drs_cal_offset = ( slice + StartCellVector[ pixel ] )%RegionOfInterest;
-        cal_pt    = pixel_pt + drs_cal_offset;
-
-        vraw = AllPixelDataVector[ slice_pt ] * dconv;
-        vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
-
-        return( vcal );
-}
-
-// this is a faster version of getValue, that does not do a full calibration
-// was used for a performace test...
-float correctDrsOffset( int slice, int pixel ){
-
-        const float dconv = 2000/4096.0;
-
-        // here 1024 is not the RegionOfInterest, but really the lenth of the pipeline
-        unsigned int physical_slice = ( slice + StartCellVector[ pixel ] ) % 1024;
-
-        unsigned int slice_pt;
-        unsigned int physical_slice_pt;
-        slice_pt                        = pixel * RegionOfInterest + slice;
-        physical_slice_pt       = pixel * RegionOfInterest + physical_slice;
-
-        float vcal = AllPixelDataVector[ slice_pt ] *
-                dconv - drs_basemean[ physical_slice_pt ];
-        return( vcal );
-}
-
 // booking and parameter settings for all histos
 void BookHistos( ){
-        // histograms for baseline extraction
-        char hName[500];
-        char hTitle[500];
-
-        TH1F *h;
-
-        for( int i = 0; i < NPIX; i++ ) {
-                sprintf(&hTitle[0],"all events all slices of pixel %d", i);
-                sprintf(&hName[0],"base%d", i);
-
-                h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );
-
-                h->GetXaxis()->SetTitle( "Sample value (mV)" );
-                h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
-                hListBaseline.Add( h );
-                hBaseline[i] = h;
-        }
-
-    hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
-    hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
-    hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
-    hList.Add( hMeanBsl );
-
-    hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
-    hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
-    hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
-    hList.Add( hpltMeanBsl );
-
-    hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
-    hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
-    hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
-    hList.Add( hRmsBsl );
-
-    hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
-    hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
-    hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
-    hList.Add( hpltRmsBsl );
-
-        hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",1440,-0.5,1439.5, 256, -27.5, 100.5);
-        hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" );
-        hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" );
-        hList.Add( hAmplSpek_cfd );
-
-        hAmplSpek_discri = new TH2F("hAmplSpek_discri","amplitude spektrum - std discriminator",1440,-0.5,1439.5, 256, -27.5, 100.5);
-        hAmplSpek_discri->GetXaxis()->SetTitle( "pixel" );
-        hAmplSpek_discri->GetXaxis()->SetTitle( "amplitude in mV" );
-        hList.Add( hAmplSpek_discri );
-
-        hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
-        hStartCell->GetXaxis()->SetTitle( "pixel" );
-        hStartCell->GetXaxis()->SetTitle( "slice" );
-        hList.Add( hStartCell );
 
         debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
@@ -755 +393 @@
                 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
         }
-        hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks", 201, -100.5, 100.5, 1024, -99.5, 100.5);
+        hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks", 
+                                        hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
+                                        4096, -48.5, 1999.5 );
         hPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
         hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
-        //hList.Add( hPeakOverlay h);
+        hList.Add( hPeakOverlay );
 
 }
 
-void SaveHistograms( char * loc_fname ){
-
-        TString fName; // name of the histogram file
-
-        // create the filename for the histogram file
-        fName = loc_fname; // use the name of the tree file
-        // TODO ... next statement doesn't work for ".fits.gz"
-        fName.Remove(fName.Length() - 5); // remove the extension .fits
-        fName += "_discri.root"; // add the new extension
-
+void SaveHistograms(const char * loc_fname ){
         // create the histogram file (replace if already existing)
-        TFile tf( fName, "RECREATE");
+        TFile tf( loc_fname, "RECREATE");
 
         hList.Write(); // write the major histograms into the top level directory
-        tf.mkdir("BaselineHisto");
-        tf.cd("BaselineHisto"); // go to new subdirectory
-        hListBaseline.Write(); // write histos into subdirectory
 
         tf.Close(); // close the file
 } // end of SaveHistograms( char * loc_fname )
 
-int FOpenDataFile(fits &datafile){
-
-//      cout << "-------------------- Data Header -------------------" << endl;
-//      datafile.PrintKeys();
-//      cout << "------------------- Data Columns -------------------" << endl;
-//      datafile.PrintColumns();
-
-        //Get the size of the data column
-        RegionOfInterest        = datafile.GetUInt("NROI");
-        NumberOfPixels          = datafile.GetUInt("NPIX");
-        //Size of column "Data" = #Pixel x ROI
-        ROIxNOP                         = datafile.GetN("Data");
-
-        //Set the sizes of the data vectors
-        AllPixelDataVector.resize(ROIxNOP,0);
-        StartCellVector.resize(NumberOfPixels,0);
-
-        //Link the data to variables
-        datafile.SetRefAddress("EventNum",  CurrentEventID);
-        datafile.SetVecAddress("Data", AllPixelDataVector);
-        datafile.SetVecAddress("StartCellData", StartCellVector);
-        datafile.GetRow(0);
-
-        return datafile.GetNumRows() ;
-}
-
-
-/////////////////////////////////////////////////////////////////////////////
-/// old BookHistos
-/*
-void BookHistos( int Samples ){
-// booking and parameter settings for all histos
-
-    h = new TH1F[ Ntypes ];
-
-    for ( int type = 0; type < Ntypes; type++){
-
-        h[ type ].SetBins(Samples, 0, Samples);
-        h[ type ].SetLineColor(1);
-        h[ type ].SetLineWidth(2);
-
-        // set X axis paras
-        h[ type ].GetXaxis()->SetLabelSize(0.1);
-        h[ type ].GetXaxis()->SetTitleSize(0.1);
-        h[ type ].GetXaxis()->SetTitleOffset(1.2);
-        h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
-
-        // set Y axis paras
-        h[ type ].GetYaxis()->SetLabelSize(0.1);
-        h[ type ].GetYaxis()->SetTitleSize(0.1);
-        h[ type ].GetYaxis()->SetTitleOffset(0.3);
-        h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
-    }
-        cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",10,10,800,600);
-        cRawAndSpikeRemoval->Divide(1, 3);
-    cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
-    cFilter->Divide(1, 3);
-
-    hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
-
-}
-*/
-