Index: /fact/tools/rootmacros/fpeak_cfd.C
===================================================================
--- /fact/tools/rootmacros/fpeak_cfd.C	(revision 12385)
+++ /fact/tools/rootmacros/fpeak_cfd.C	(revision 12385)
@@ -0,0 +1,341 @@
+#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>
+#include <deque>
+
+#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"
+
+#include "FOpenDataFile.h"
+#include "FOpenDataFile.c"
+
+#include "DrsCalibration.C"
+#include "DrsCalibration.h"
+
+#include "SpikeRemoval.h"
+#include "SpikeRemoval.C"
+
+bool breakout=false;
+
+int NEvents;
+vector<int16_t> AllPixelDataVector;
+vector<int16_t> StartCellVector;
+unsigned int CurrentEventID;
+size_t PXLxROI;
+UInt_t RegionOfInterest;
+UInt_t NumberOfPixels;
+
+size_t drs_n;
+vector<float> drs_basemean;
+vector<float> drs_gainmean;
+vector<float> drs_triggeroffsetmean;
+
+vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
+vector<float> Vcorr(FAD_MAX_SAMPLES);  // corrected Values
+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
+
+// histograms
+const int NumberOfDebugHistoTypes = 7;
+const unsigned int
+	Ameas_ 	= 0,
+	N1mean_ = 1,
+	Vcorr_  = 2,
+	Vtest_  = 3,
+	Vslide_ = 4,
+	Vcfd_   = 5,
+	Vcfd2_  = 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 * hAmplSpek_cfd;
+TObjArray hList;
+
+void BookHistos( );
+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 		= "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 verbosityLevel = 1, // different verbosity levels can be implemented here
+	bool ProduceGraphic = true
+	)
+{
+gStyle->SetPalette(1,0);
+gROOT->SetStyle("Plain");
+
+	TCanvas *cFiltered = NULL;
+	if (spikeDebug){
+		cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
+		cFiltered->Divide(1, 4);
+	}
+
+
+
+// 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.;
+
+	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 FOpenDataFile:" << 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
+	FOpenCalibFile(	drsfilename,
+					drs_basemean,
+					drs_gainmean,
+					drs_triggeroffsetmean,
+					drs_n);
+
+	BookHistos( );
+
+	for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
+		// Get an Event --> consists of 1440 Pixel ...erm....data
+		datafile->GetRow( ev );
+
+		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
+			// in Vcfd2
+			// first Argument 1 means ... *rising* edge
+			// second Argument 1 means ... search with stepsize 1 ... 8 is okay as well
+			vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 8);
+			// zXings means "zero cross ings"
+			EnlargeRegion( *zXings, 10, 10);
+			findAbsMaxInRegions( *zXings, Vslide);
+			removeMaximaBelow( *zXings, 3.0, 0);
+			removeRegionWithMaxOnEdge( *zXings, 2);
+			removeRegionOnFallingEdge( *zXings, 100);
+
+			// fill maxima in Histogram
+			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 ){
+		
+					// 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[Ameas_].Draw();
+				
+				cFiltered->cd(2);
+				gPad->SetGrid();
+				debugHistos[Vcorr_].Draw();
+
+				cFiltered->cd(3);
+				gPad->SetGrid();
+				debugHistos[Vslide_].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();
+				}
+
+				cFiltered->cd(4);
+				gPad->SetGrid();
+				debugHistos[Vcfd2_].Draw();
+				TLine *zeroline = new TLine(0, 0, 1024, 0);
+				zeroline->SetLineColor(kBlue);
+				zeroline->Draw();
+
+				cFiltered->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;
+				}
+
+				//TODO!!!!!!!!!
+				// do some Garbage collection ...
+			}// end of if(spikeDebug)
+
+			delete zXings;
+			if (breakout)
+				break;
+
+		} // end of loop over pixels
+
+		if (breakout)
+			break;
+	}	// end of loop over pixels
+if (ProduceGraphic){
+	TCanvas * cSpektrum;
+	cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400);
+	cSpektrum->Divide(1,1);
+  cSpektrum->cd(1);
+  hAmplSpek_cfd->Draw("COLZ");
+  cSpektrum->Modified();
+  cSpektrum->Update();
+}
+	SaveHistograms( OutRootFileName );
+	return( 0 );
+}
+
+// booking and parameter settings for all histos
+void BookHistos( ){
+	// histograms for baseline extraction
+
+	TString name,title;
+	title = "all events all slices of pixel ";
+	name = "base";
+
+	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 );
+
+
+	debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
+	for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
+		debugHistos[ type ].SetBins(1024, 0, 1024);
+		debugHistos[ type ].SetLineColor(1);
+		debugHistos[ type ].SetLineWidth(2);
+
+		// set X axis paras
+		debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
+		debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
+		debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
+		debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
+
+		// set Y axis paras
+		debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
+		debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
+		debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
+		debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
+	}
+}
+
+void SaveHistograms( const char * loc_fname ){
+	// create the histogram file (replace if already existing)
+	TFile tf( loc_fname, "RECREATE");
+
+	hList.Write(); // write the major histograms into the top level directory
+
+	tf.Close(); // close the file
+}