Index: ct/tools/rootmacros/fpeak_cdf.C
===================================================================
--- /fact/tools/rootmacros/fpeak_cdf.C	(revision 12385)
+++ 	(revision )
@@ -1,341 +1,0 @@
-#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
-}