Index: fact/tools/rootmacros/fpeak_discri.C
===================================================================
--- fact/tools/rootmacros/fpeak_discri.C	(revision 12262)
+++ fact/tools/rootmacros/fpeak_discri.C	(revision 12262)
@@ -0,0 +1,689 @@
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TProfile.h>
+#include <TTimer.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <Getline.h>
+#include <TLine.h>
+#include <TBox.h>
+#include <TMath.h>
+#include <TFile.h>
+#include <TStyle.h>
+
+#include <stdio.h>
+#include <stdint.h>
+#include <cstdio>
+
+#define NPIX  1440
+#define NCELL 1024
+#define FAD_MAX_SAMPLES 1024
+
+#define HAVE_ZLIB
+#include "fits.h"
+#include "FOpenCalibFile.c"
+
+#include "discriminator.h"
+#include "discriminator.C"
+#include "zerosearch.h"
+#include "zerosearch.C"
+#include "factfir.C"
+
+
+vector<int16_t> AllPixelDataVector;
+vector<int16_t> StartCellVector;
+
+unsigned int CurrentEventID;
+
+bool breakout=false;
+
+size_t ROIxNOP;
+UInt_t NumberOfPixels;
+UInt_t RegionOfInterest;
+int NEvents;
+
+size_t drs_n;
+vector<float> drs_basemean;
+vector<float> drs_gainmean;
+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
+
+// not needed?
+//vector<float> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors
+
+
+float getValue( int, int );
+float correctDrsOffset( int slice, int pixel );
+void computeN1mean( int );
+void removeSpikes( int );
+
+// histograms
+const int Ntypes = 7;
+const unsigned int  // arranged by Dominik
+	tAmeas 	= 0, 
+	tN1mean = 1, 
+	tVcorr  = 2,
+	tVtest  = 3,
+	tVslide = 4,
+	tVcfd   = 5,
+	tVcfd2  = 6;
+
+TH1F* h;
+TH1F *oldh;
+TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
+                    // x = pixel id, y = DRS cell id
+TH2F hPixelCellData(
+		"PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
+TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
+TH1F *hMeanBsl, *hpltMeanBsl;
+TH1F *hRmsBsl, *hpltRmsBsl; 
+TH2F * hAmplSpek_cfd;
+TH2F * hAmplSpek_discri;
+TObjArray hList;
+TObjArray hListBaseline;
+
+void BookHistos( );
+void SaveHistograms( char * );
+
+// Create a canvas
+TCanvas* CW;
+TCanvas* cFilter;
+
+
+int searchSinglesPeaks = 0;
+
+int fpeak_discri( 
+	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
+	)
+{
+
+gStyle->SetPalette(1,0);
+gROOT->SetStyle("Plain");
+// read FACT raw data 
+// 	* remove spikes
+//	* calculate baseline
+//	* find peaks (CFD and normal discriminator)
+//	* compute pulse height and pulse integral spektrum of the peaks
+
+// sliding window filter settings
+	int k_slide = 16;
+	vector<double> a_slide(k_slide, 1);
+	double b_slide = k_slide;
+
+// CFD filter settings
+	int k_cfd = 10;
+	vector<double> a_cfd(k_cfd, 0);
+	double b_cfd = 1.;
+	a_cfd[0]=-0.75;
+	a_cfd[k_cfd-1]=1.;
+
+// 2nd slinding window filter
+	int ks2 = 16;
+	vector<double> as2(ks2, 1);
+	double bs2 = ks2;
+	
+// 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);
+  if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
+
+//Get the DRS calibration
+	FOpenCalibFile(	drsfilename,
+					drs_basemean,
+					drs_gainmean,
+					drs_triggeroffsetmean,
+					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( );
+	TCanvas * cSpektrum = new TCanvas();
+	cSpektrum->Divide(1,2);
+	cSpektrum->cd(1);
+	hAmplSpek_discri->Draw("COLZ");
+	cSpektrum->cd(2);
+	hAmplSpek_cfd->Draw("COLZ");
+//	TCanvas * cStartCell = new TCanvas();
+//	cStartCell->cd();
+//	hStartCell->Draw();
+//	hPixelCellData.Draw();
+
+	float calibratedVoltage;
+	for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
+		datafile->GetRow( ev );
+
+
+		for ( int pix = 0; pix < 1440; pix++ ){
+			// This means: GetEvent
+
+			// this is a stupid hack ... these is more code at the
+			// end of this loop to complete this hack ...
+			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;
+				}
+			}
+
+			//hStartCell->Fill( pix, StartCellVector[pix] );
+
+			for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
+				//calibratedVoltage = getValue( sl, pix);
+				calibratedVoltage = correctDrsOffset( sl, pix);
+				if (verbosityLevel > 10){
+					printf("calibratedVoltage = %f\n", calibratedVoltage);
+				}
+				Ameas[ sl ] = calibratedVoltage;
+				if (spikeDebug){
+					oldh[tAmeas].SetBinContent(sl, calibratedVoltage);
+				}
+			}
+
+			computeN1mean( RegionOfInterest );
+			// operates on Ameas[] and writes to N1mean[]
+			removeSpikes( RegionOfInterest );
+			// operates on Ameas[] and N1mean[], then writes to Vcorr[]
+
+			if (spikeDebug){
+				for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
+					hPixelCellData.Fill(sl, Vcorr[ sl ]);
+				}
+			}
+
+			// 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);
+
+			float myTHR = 3.5;
+			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);
+			}
+
+			// peaks in Ameas[] are found by searching for zero crossings
+			// in Vcfd2
+			// first Argument 1 means ... *rising* edge
+			// second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
+			vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 1);
+
+			// zXings means "zero cross ings"
+			ShiftRegionBy(*zXings, -ks2/2);
+			EnlargeRegion(*zXings, 10, 10);
+			findAbsMaxInRegions(*zXings, Vslide);
+
+			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++){
+					oldh[tVslide].SetBinContent( sl, Vslide[sl] );
+					oldh[tVcfd].SetBinContent( sl, Vcfd[sl] );
+					oldh[tVcfd2].SetBinContent( sl, Vcfd2[sl] );
+				}
+
+				CW->cd( 1);
+				gPad->SetGrid();
+				oldh[tAmeas].Draw();
+
+//				CW->cd( tN1mean + 1);
+//				gPad->SetGrid();
+//				oldh[tN1mean].Draw();
+
+				CW->cd( 2);
+				gPad->SetGrid();
+				oldh[tVcorr].Draw();
+
+				//cFilter->cd( Ntypes - tVslide );
+				cFilter->cd(1);
+				gPad->SetGrid();
+				oldh[tVslide].Draw();
+				TLine thrLine(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;
+				for (unsigned int i=0; i<zXings->size(); i++){
+					OneBox = new TBox(
+						zXings->at(i).maxPos -10 ,
+						zXings->at(i).maxVal -0.5,
+						zXings->at(i).maxPos +10 ,
+						zXings->at(i).maxVal +0.5);
+					OneBox->SetLineColor(kBlue);
+					OneBox->SetLineWidth(1);
+					OneBox->SetFillStyle(0);
+					OneBox->SetFillColor(kRed);
+					MyBoxes.push_back(OneBox);
+					OneBox->Draw();
+				}
+
+//				cFilter->cd( Ntypes - tVcfd );
+//				cFilter->cd(2);
+//				gPad->SetGrid();
+//				oldh[tVcfd].Draw();
+//				TLine zeroline(0, 0, 1024, 0);
+//				zeroline.SetLineColor(kBlue);
+//				zeroline.Draw();
+
+				//cFilter->cd( Ntypes - tVcfd2 );
+				cFilter->cd(2);
+				gPad->SetGrid();
+				oldh[tVcfd2].Draw();
+				TLine zeroline(0, 0, 1024, 0);
+				zeroline.SetLineColor(kBlue);
+				zeroline.Draw();
+
+
+
+
+
+				CW->Update();
+				cFilter->Update();
+
+				//Process gui events asynchronously during input
+				TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
+				timer.TurnOn();
+				TString input = Getline("Type 'q' to exit, <return> to go on: ");
+				timer.TurnOff();
+				if (input=="q\n") {
+					breakout=true;
+					delete OneLine;
+//					for (unsigned int i=0; i<MyLines.size(); i++){
+//						delete MyLines[i];
+//					}
+//					for (unsigned int i=0; i<MyBoxes.size(); i++){
+//						delete MyBoxes[i];
+//					}
+				}
+			}// end of if(spikeDebug)
+			delete Regions;
+			delete zXings;
+
+			if (PixelID != -1){
+				pix = 2000;
+			}
+
+			if (breakout)
+				break;
+
+		} // end of loop over pixels
+
+	if (ev % 10 ==0){
+		cSpektrum->Modified();
+		cSpektrum->Update();
+//		cStartCell->Modified();
+//		cStartCell->Update();
+	}
+		if (breakout)
+			break;
+	}	// end of loop over pixels
+
+
+//	delete cStartCell;
+	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++ ) oldh[ tVcorr ].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
+
+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 );
+}
+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 );
+}
+
+void BookHistos( ){
+// booking and parameter settings for all histos
+
+	// 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 );
+
+    oldh = new TH1F[ Ntypes ];
+
+    for ( int type = 0; type < Ntypes; type++){
+
+        oldh[ type ].SetBins(1024, 0, 1024);
+        oldh[ type ].SetLineColor(1);
+        oldh[ type ].SetLineWidth(2);
+
+        // set X axis paras
+        oldh[ type ].GetXaxis()->SetLabelSize(0.1);
+        oldh[ type ].GetXaxis()->SetTitleSize(0.1);
+        oldh[ type ].GetXaxis()->SetTitleOffset(1.2);
+        oldh[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
+
+        // set Y axis paras
+        oldh[ type ].GetYaxis()->SetLabelSize(0.1);
+        oldh[ type ].GetYaxis()->SetTitleSize(0.1);
+        oldh[ type ].GetYaxis()->SetTitleOffset(0.3);
+        oldh[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
+    }
+
+    CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
+	CW->Divide(1, 2);
+    cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
+	cFilter->Divide(1, 2);
+}
+
+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
+
+	// create the histogram file (replace if already existing)
+	TFile tf( 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.)");
+    }
+    CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
+    CW->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);
+
+}
+*/
+