Index: fact/tools/rootmacros/fpeak_cdf.C
===================================================================
--- fact/tools/rootmacros/fpeak_cdf.C	(revision 12357)
+++ fact/tools/rootmacros/fpeak_cdf.C	(revision 12369)
@@ -31,16 +31,19 @@
 #include "factfir.C"
 
-
+#include "FOpenDataFile.h"
+#include "FOpenDataFile.c"
+
+#include "DrsCalibration.C"
+#include "DrsCalibration.h"
+
+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;
@@ -49,7 +52,4 @@
 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
@@ -62,5 +62,5 @@
 
 
-float getValue( int, int );
+//float getValue( int, int );
 void computeN1mean( int );
 void removeSpikes( int );
@@ -79,6 +79,6 @@
 TH1F* h;
 TH1F *debugHistos;
-TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
-                    // x = pixel id, y = DRS cell id
+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
@@ -92,18 +92,28 @@
 
 void BookHistos( );
-void SaveHistograms( char * );
-
-int searchSinglesPeaks = 0;
-
+void SaveHistograms(const char * );
+
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// read FACT raw data
+// 	* remove spikes
+//	* calculate baseline
+//	* find peaks (CFD and normal discriminator)
+//	* compute pulse height and pulse integral spektrum of the peaks
 int fpeak(
 	char *datafilename 		= "../../20111011_055.fits.gz",
 	const char *drsfilename = "../../20111011_054.drs.fits.gz",
+	const char *OutRootFileName = "./fpeak_cdf.Coutput.root",
+	int firstevent 			= 0,
 	int nevents 			= -1,
-	int firstevent 			= 0,
-	int PixelID				= -1,
+	int firstpixel			= 0,
+	int npixel				= -1,
 	bool spikeDebug = false,
 	int verbosityLevel = 1 // different verbosity levels can be implemented here
 	)
 {
+gStyle->SetPalette(1,0);
+gROOT->SetStyle("Plain");
 
 // Create (pointer to) Canvases, which are used in every run,
@@ -130,16 +140,5 @@
 	}
 
-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
@@ -150,22 +149,29 @@
 	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;
+	fits * datafile;
+	NEvents = OpenDataFile( datafilename, &datafile,
+		AllPixelDataVector, StartCellVector, CurrentEventID,
+		RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
+	if (NEvents == 0){
+		cout << "return code of FOpenDataFile:" << datafilename<< endl;
+		cout << "is zero -> aborting." << endl;
+		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
+
+	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 <<"Totel # of Pixel: "<< NumberOfPixels << "\t";
+	if ( npixel == -1 || npixel > (int)NumberOfPixels  ) npixel = (int)NumberOfPixels; // -1 means all!
+	if (verbosityLevel > 0)
+		cout <<"of, which "<< nevents << "will be processed"<< endl;
+
+
+
+	//Get the DRS calibration
 	FOpenCalibFile(	drsfilename,
 					drs_basemean,
@@ -174,24 +180,5 @@
 					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 calibratedVoltage;
 
 	for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
@@ -199,19 +186,5 @@
 		datafile->GetRow( ev );
 
-		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;
-				}
-			}
-
+		for ( int pix = firstpixel; pix < firstpixel+npixel; pix++ ){
 			if (verbosityLevel > 0){
 				if (pix % 20 ==0){
@@ -221,18 +194,8 @@
 			}
 
-			// 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);
-				}
-			}
+			applyDrsCalibration( Ameas,pix,
+				drs_basemean, drs_gainmean, drs_triggeroffsetmean,
+				RegionOfInterest, AllPixelDataVector, StartCellVector);
+
 			// operates on Ameas[] and writes to N1mean[]
 			computeN1mean( RegionOfInterest );
@@ -260,7 +223,8 @@
 			EnlargeRegion(*zXings, 10, 10);
 			findAbsMaxInRegions(*zXings, Vslide);
-			removeMaximaBelow( *zXings, 11.0, 0);
-			removeMaximaAbove( *zXings, 14.0, 0);
-
+			//removeMaximaBelow( *zXings, 11.0, 0);
+			//removeMaximaAbove( *zXings, 14.0, 0);
+
+			// fill maxima in Histogram
 			if (zXings->size() != 0 ){
 				for (unsigned int i=0; i<zXings->size(); i++){
@@ -280,4 +244,5 @@
 			if ( spikeDebug ){
 				for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
+					debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
 					debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
 					debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
@@ -344,12 +309,4 @@
 
 			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 = 1440;
-			}
-
-
 			if (breakout)
 				break;
@@ -373,6 +330,6 @@
 
 		cTemplate->cd();
-		hTemp_Array[PixelID]->GetYaxis()->SetRangeUser(-10,40);
-		hTemp_Array[PixelID]->Draw("COLZ");
+		hTemp_Array[firstpixel]->GetYaxis()->SetRangeUser(-10,40);
+		hTemp_Array[firstpixel]->Draw("COLZ");
 		cTemplate->Modified();
 		cTemplate->Update();
@@ -385,5 +342,5 @@
 
 
-	SaveHistograms( datafilename );
+	SaveHistograms( OutRootFileName );
 
 	return( 0 );
@@ -401,5 +358,4 @@
     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];
@@ -411,12 +367,8 @@
             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; 
             }
@@ -425,6 +377,4 @@
                 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
@@ -458,26 +408,4 @@
 } // 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 );
-}
 
 // booking and parameter settings for all histos
@@ -487,11 +415,16 @@
 	char hTitle[500];
 
+	TString name,title;
+	title = "all events all slices of pixel ";
+	name = "base";
 	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 );
+		title += i;
+		name += i;
+//		sprintf(&hTitle[0],"all events all slices of pixel %d", i);
+//		sprintf(&hName[0],"base%d", i);
+
+		h = new TH1F( name, title, 400, -99.5 ,100.5 );
 
 		h->GetXaxis()->SetTitle( "Sample value (mV)" );
@@ -530,9 +463,4 @@
 	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 ];
@@ -571,16 +499,7 @@
 }
 
-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
@@ -590,33 +509,7 @@
 	tf.cd("..");
 	tf.mkdir("PulseTempplates");
-	tf.cd("PulseTempplates");
+	tf.cd("PulseTemplates");
 	hListTemplates.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() ;
-}
+} // end of SaveHistograms(const char * loc_fname )
