Index: fact/tools/rootmacros/FIntFixedPosAllPx.c
===================================================================
--- fact/tools/rootmacros/FIntFixedPosAllPx.c	(revision 12166)
+++ fact/tools/rootmacros/FIntFixedPosAllPx.c	(revision 12166)
@@ -0,0 +1,54 @@
+#include <cstdio>
+int FIntFixedPosAllPx(fits &datafile, vector<int16_t> &data, vector<int16_t> &data_offset, unsigned int &data_num, UInt_t data_px, vector<float> &drs_basemean, vector<float> &drs_gainmean, vector<float> &drs_triggeroffsetmean, UInt_t data_roi, UInt_t integration_size, UInt_t integration_delay, TH1F* spectrum[])
+//Data, calibration, data_roi, data_num, data_px, threshold, two histograms
+{
+//	UInt_t integration_size = 5;
+//	UInt_t integration_delay = 238;
+	
+//	float sample, last_sample, integral;
+	float integral;
+	UInt_t drs_calib_offset;
+
+//**********************************************************************************	
+	for (size_t i=0; (i<datafile.GetNumRows()) && (i<100000); i++)
+	{
+		datafile.GetRow(i);
+		cout << "Event number: " << data_num << endl;
+		
+		//Iterate over the pixels
+		for (int j=0; j<data_px; j++)
+		{
+//		size_t j=pixelnr; //Fix the Pixel to a SoftID
+		
+			//Integrate
+			integral = 0;
+			for(UInt_t l=integration_delay; l<integration_delay+integration_size; l++)
+			{
+				drs_calib_offset = (l+data_offset[j])%data_roi;
+				integral+=(data[j*data_roi+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+			}
+			spectrum[j]->Fill(integral);
+//			for(Int_t l=-10; l<60; l++)
+//			{
+//				drs_calib_offset = (l+data_offset[j])%data_roi;
+//				sample = (data[j*data_roi+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+//			}
+//			k+=min_dist;
+//			
+//			drs_calib_offset = (start_sample-1+data_offset[j])%data_roi;
+//			last_sample = (data[j*data_roi+start_sample-1]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+start_sample-1])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+//			for (UInt_t k=start_sample; k<end_sample; k++)
+//			{
+//				drs_calib_offset = (k+data_offset[j])%data_roi;
+//				sample = (data[j*data_roi+k]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+//				if((last_sample<threshold)&&(sample>threshold))
+//				{
+//				}
+//				//TBD: mistake that after the deadtime the last_sample must be set new...
+//				//Solution: first process full pipeline
+//				last_sample = sample;
+//			}
+		}
+	}
+	return 0;
+}
Index: fact/tools/rootmacros/FOpenCalibFile.c
===================================================================
--- fact/tools/rootmacros/FOpenCalibFile.c	(revision 12166)
+++ fact/tools/rootmacros/FOpenCalibFile.c	(revision 12166)
@@ -0,0 +1,45 @@
+#include <cstdio>
+int FOpenCalibFile(const char *drsname, vector<float> &drs_basemean, vector<float> &drs_gainmean, vector<float> &drs_triggeroffsetmean, size_t &drs_n)
+{
+//-------------------------------------------
+//Open the file
+//-------------------------------------------
+	fits drsfile(drsname);
+	
+	if (!drsfile)
+	{
+	cout << "Couldn't properly open the drsfile." << endl;
+	return 1;
+	}
+	
+//-------------------------------------------
+//Print the header
+//-------------------------------------------
+	cout << "----------------- DRS Calib Header -----------------" << endl;
+	drsfile.PrintKeys();
+	cout << "---------------- DRS Calib Columns -----------------" << endl;
+	drsfile.PrintColumns();
+	
+//-------------------------------------------
+//Get the sizes of the data column
+//-------------------------------------------
+	drs_n = drsfile.GetN("BaselineMean");
+	
+//-------------------------------------------
+//Set the sizes of the DRS data vectors
+//-------------------------------------------
+	drs_basemean.resize(drs_n,0);
+	drs_gainmean.resize(drs_n,0);
+	drs_triggeroffsetmean.resize(drs_n,0);
+	
+//-------------------------------------------
+//Link the data to variables
+//-------------------------------------------
+	drsfile.SetVecAddress("BaselineMean", drs_basemean);
+	drsfile.SetVecAddress("GainMean", drs_gainmean);
+	drsfile.SetVecAddress("TriggerOffsetMean", drs_triggeroffsetmean);
+	drsfile.GetRow(0); //Read the calibration data
+	
+	cout << "Reading calibration file successful..." << endl;
+	return 0;
+}
Index: fact/tools/rootmacros/FOpenDataFile.c
===================================================================
--- fact/tools/rootmacros/FOpenDataFile.c	(revision 12166)
+++ fact/tools/rootmacros/FOpenDataFile.c	(revision 12166)
@@ -0,0 +1,35 @@
+#include <cstdio>
+int FOpenDataFile(fits &datafile, vector<int16_t> &data, vector<int16_t> &data_offset, unsigned int &data_num, size_t &data_n, UInt_t &data_roi, UInt_t &data_px)
+{
+//-------------------------------------------
+//Print the header
+//-------------------------------------------
+/*	cout << "-------------------- Data Header -------------------" << endl;
+	datafile.PrintKeys();
+	cout << "------------------- Data Columns -------------------" << endl;
+	datafile.PrintColumns();
+	*/
+//-------------------------------------------
+//Get the size of the data column
+//-------------------------------------------
+	data_roi = datafile.GetUInt("NROI");  // Value from header
+	data_px = datafile.GetUInt("NPIX");
+	data_n = datafile.GetN("Data");         //Size of column "Data" = #Pixel x ROI
+	
+//-------------------------------------------
+//Set the sizes of the data vectors
+//-------------------------------------------
+	data.resize(data_n,0);
+	data_offset.resize(data_px,0);
+	
+//-------------------------------------------
+//Link the data to variables
+//-------------------------------------------
+	datafile.SetRefAddress("EventNum", data_num);
+	datafile.SetVecAddress("Data", data);
+	datafile.SetVecAddress("StartCellData", data_offset);
+	datafile.GetRow(0);
+	
+	cout << "Opening data file successful..." << endl;
+	return 0;
+}
Index: fact/tools/rootmacros/FOscilloscope.c
===================================================================
--- fact/tools/rootmacros/FOscilloscope.c	(revision 12166)
+++ fact/tools/rootmacros/FOscilloscope.c	(revision 12166)
@@ -0,0 +1,59 @@
+#include <cstdio>
+int FOscilloscope(fits &datafile, vector<int16_t> &data, vector<int16_t> &data_offset, unsigned int &data_num, UInt_t data_px, vector<float> &drs_basemean, vector<float> &drs_gainmean, vector<float> &drs_triggeroffsetmean, UInt_t data_roi, float threshold, TH2* pulseshape, TH1* spectrum)
+//Data, calibration, data_roi, data_num, data_px, threshold, two histograms
+{
+	UInt_t start_sample = 10;
+	UInt_t end_sample = data_roi-60;
+	UInt_t min_dist = 150;
+	UInt_t integration_size = 10;
+	UInt_t integration_delay = 5;
+	
+	float sample, last_sample, integral;
+	UInt_t drs_calib_offset;
+
+//**********************************************************************************	
+	for (size_t i=0; i<datafile.GetNumRows(); i++)
+	{
+		datafile.GetRow(i);
+		cout << "Event number: " << data_num << endl;
+		
+		//Iterate over the pixels
+		for (int j=0; j<data_px; j++)
+		{
+//		size_t j=pixelnr; //Fix the Pixel to a SoftID
+		
+			//Iterate over the slices
+			drs_calib_offset = (start_sample-1+data_offset[j])%data_roi;
+			last_sample = (data[j*data_roi+start_sample-1]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+start_sample-1])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+			for (UInt_t k=start_sample; k<end_sample; k++)
+			{
+				drs_calib_offset = (k+data_offset[j])%data_roi;
+				sample = (data[j*data_roi+k]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+				if((last_sample<threshold)&&(sample>threshold))
+				{
+					integral = 0;
+					for(UInt_t l=integration_delay; l<integration_delay+integration_size; l++)
+					{
+						drs_calib_offset = (k+l+data_offset[j])%data_roi;
+						integral+=(data[j*data_roi+k+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+					}
+//					if((integral>60)&&(integral<=300)) {
+//					if((integral>60)&&(integral<=300)) {
+						for(Int_t l=-10; l<60; l++)
+						{
+							drs_calib_offset = (k+l+data_offset[j])%data_roi;
+							sample = (data[j*data_roi+k+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+							pulseshape->Fill(l,sample);
+						}
+//					}
+					spectrum->Fill(integral);
+					k+=min_dist;
+				}
+				//TBD: mistake that after the deadtime the last_sample must be set new...
+				//Solution: first process full pipeline
+				last_sample = sample;
+			}
+		}
+	}
+	return 0;
+}
Index: fact/tools/rootmacros/FOscilloscopeAllPx.c
===================================================================
--- fact/tools/rootmacros/FOscilloscopeAllPx.c	(revision 12166)
+++ fact/tools/rootmacros/FOscilloscopeAllPx.c	(revision 12166)
@@ -0,0 +1,59 @@
+#include <cstdio>
+int FOscilloscopeAllPx(fits &datafile, vector<int16_t> &data, vector<int16_t> &data_offset, unsigned int &data_num, UInt_t data_px, vector<float> &drs_basemean, vector<float> &drs_gainmean, vector<float> &drs_triggeroffsetmean, UInt_t data_roi, float threshold, TH1F* spectrum[])
+//Data, calibration, data_roi, data_num, data_px, threshold, two histograms
+{
+	UInt_t start_sample = 10;
+	UInt_t end_sample = data_roi-60;
+	UInt_t min_dist = 150;
+	UInt_t integration_size = 10;
+	UInt_t integration_delay = 5;
+	
+	float sample, last_sample, integral;
+	UInt_t drs_calib_offset;
+
+//**********************************************************************************	
+	for (size_t i=0; i<datafile.GetNumRows(); i++)
+	{
+		datafile.GetRow(i);
+		cout << "Event number: " << data_num << endl;
+		
+		//Iterate over the pixels
+		for (int j=0; j<data_px; j++)
+		{
+//		size_t j=pixelnr; //Fix the Pixel to a SoftID
+		
+			//Iterate over the slices
+			drs_calib_offset = (start_sample-1+data_offset[j])%data_roi;
+			last_sample = (data[j*data_roi+start_sample-1]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+start_sample-1])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+			for (UInt_t k=start_sample; k<end_sample; k++)
+			{
+				drs_calib_offset = (k+data_offset[j])%data_roi;
+				sample = (data[j*data_roi+k]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+				if((last_sample<threshold)&&(sample>threshold))
+				{
+					integral = 0;
+					for(UInt_t l=integration_delay; l<integration_delay+integration_size; l++)
+					{
+						drs_calib_offset = (k+l+data_offset[j])%data_roi;
+						integral+=(data[j*data_roi+k+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+					}
+//					if((integral>60)&&(integral<=300)) {
+//					if((integral>60)&&(integral<=300)) {
+					//Iterate for the pulseshape
+//					for(Int_t l=-10; l<60; l++)
+//					{
+//						drs_calib_offset = (k+l+data_offset[j])%data_roi;
+//						sample = (data[j*data_roi+k+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+//					}
+//					}
+					spectrum[j]->Fill(integral);
+					k+=min_dist;
+				}
+				//TBD: mistake that after the deadtime the last_sample must be set new...
+				//Solution: first process full pipeline
+				last_sample = sample;
+			}
+		}
+	}
+	return 0;
+}
Index: fact/tools/rootmacros/FPedestal.c
===================================================================
--- fact/tools/rootmacros/FPedestal.c	(revision 12166)
+++ fact/tools/rootmacros/FPedestal.c	(revision 12166)
@@ -0,0 +1,35 @@
+#include <cstdio>
+int FPedestal(char *title_base, fits &datafile, vector<int16_t> &data, vector<int16_t> &data_offset, vector<float> &drs_basemean, vector<float> &drs_gainmean, vector<float> &drs_triggeroffsetmean, UInt_t data_roi, size_t pixelnr, float &pedestal_mean, float &pedestal_rms)
+{
+//	TCanvas *canv_base = new TCanvas( "canv_base", "Determine the baseline", 100, 520, 700, 500 );
+	TH1F *baseline = new TH1F("baseline", title_base, 1000,-10,40);
+//	canv_base->SetLogy();
+	baseline->GetXaxis()->SetTitle("Sample value (mV)");
+	baseline->GetYaxis()->SetTitle("Entries");
+	
+//-------------------------------------------
+//Find the baseline
+//-------------------------------------------
+	for (size_t i=0; i<datafile.GetNumRows(); i++)
+	{
+		datafile.GetRow(i);
+		
+		//Iterate over the slices
+		for (UInt_t k=0; k<data_roi; k++)
+		{
+			UInt_t drs_calib_offset = (k+data_offset[pixelnr])%data_roi;
+			float sample = (data[pixelnr*data_roi+k]*2000/4096.-drs_basemean[pixelnr*data_roi+drs_calib_offset]-drs_triggeroffsetmean[pixelnr*data_roi+k])/drs_gainmean[pixelnr*data_roi+drs_calib_offset]*1907.35;
+			baseline->Fill(sample);
+		}
+	}
+	pedestal_mean = baseline->GetXaxis()->GetBinCenter(baseline->GetMaximumBin());
+	pedestal_rms = baseline->GetRMS();
+	std::cout << "Value of maximal probability: " << pedestal_mean << " +- " << pedestal_rms << std::endl;
+	
+//	canv_base->cd();
+//	baseline->Draw();
+//	canv_base->Modified();
+//	canv_base->Update();
+	
+	return 0;
+}
Index: fact/tools/rootmacros/FPedestalAllPx.c
===================================================================
--- fact/tools/rootmacros/FPedestalAllPx.c	(revision 12166)
+++ fact/tools/rootmacros/FPedestalAllPx.c	(revision 12166)
@@ -0,0 +1,62 @@
+#include <cstdio>
+int FPedestalAllPx(fits &datafile, vector<int16_t> &data, vector<int16_t> &data_offset, vector<float> &drs_basemean, vector<float> &drs_gainmean, vector<float> &drs_triggeroffsetmean, UInt_t data_px, UInt_t data_roi, float pedestal_mean[], float pedestal_rms[])
+{
+//	TCanvas *canv_base = new TCanvas( "canv_base", "Determine the baseline", 100, 520, 700, 500 );
+//	canv_base->SetLogy();
+
+	char title_base[500];
+	char name_base[50];
+	TH1F* baseline[data_px];
+	for(int i=0; i<data_px; i++) {
+		std::sprintf(title_base,"Baseline of Px %i (data: %s, DRS: %s)",i);
+		std::sprintf(name_base,"base%i",i);
+		baseline[i] = new TH1F(name_base,title_base,2000,-99.5,100.5);
+		baseline[i]->GetXaxis()->SetTitle("Sample value (mV)");
+		baseline[i]->GetYaxis()->SetTitle("Entries");
+	}
+	
+//	TH1F *baseline = new TH1F("baseline", title_base, 1000,-10,40);
+//	baseline[i]->GetXaxis()->SetTitle("Sample value (mV)");
+//	baseline[i]->GetYaxis()->SetTitle("Entries");
+	
+//-------------------------------------------
+//Find the baseline
+//-------------------------------------------
+	for (size_t i=0; i<datafile.GetNumRows(); i++)
+	{
+		datafile.GetRow(i);
+		cout << "Loop variable: " << i << endl;
+		
+		//Iterate over the pixels
+		for (int j=0; j<data_px; j++)
+		{
+//		size_t j=pixelnr; //Fix the Pixel to a SoftID
+
+			//Iterate over the slices
+			for (UInt_t k=0; k<data_roi; k++)
+			{
+				UInt_t drs_calib_offset = (k+data_offset[j])%data_roi;
+				float sample = (data[j*data_roi+k]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+				baseline[j]->Fill(sample);
+				
+				//Original code
+				//UInt_t drs_calib_offset = (k+data_offset[j])%data_roi;
+				//float sample = (data[j*data_roi+k]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+				//baseline->Fill(sample);
+			}
+		}
+	}
+	
+	for (int j=0; j<data_px; j++) {
+		pedestal_mean[j] = baseline[j]->GetXaxis()->GetBinCenter(baseline[j]->GetMaximumBin());
+		pedestal_rms[j] = baseline[j]->GetRMS();
+	}
+	std::cout << "Value of maximal probability: " << pedestal_mean[1348] << " +- " << pedestal_rms[1348] << std::endl;
+	
+//	canv_base->cd();
+//	baseline->Draw();
+//	canv_base->Modified();
+//	canv_base->Update();
+	
+	return 0;
+}
Index: fact/tools/rootmacros/TPKplotevent.c
===================================================================
--- fact/tools/rootmacros/TPKplotevent.c	(revision 12166)
+++ fact/tools/rootmacros/TPKplotevent.c	(revision 12166)
@@ -0,0 +1,20 @@
+#include <cstdio>
+
+int TPKplotevent(char *title, vector<int16_t> &data, vector<int16_t> &data_offset, vector<float> &drs_basemean, vector<float> &drs_gainmean, vector<float> &drs_triggeroffsetmean, UInt_t data_roi, size_t pixelnr)
+{
+	TCanvas *canv = new TCanvas( "canv", "Mean values of the first event", 100, 10, 700, 500 );
+	TProfile *pix = new TProfile("pix", title, 1024, -0.5, 1023.5);
+
+	for (UInt_t k=0; k<data_roi; k++)
+	{
+		UInt_t drs_calib_offset = (k+data_offset[pixelnr])%data_roi;
+		float sample = (data[pixelnr*data_roi+k]*2000/4096.-drs_basemean[pixelnr*data_roi+drs_calib_offset]-drs_triggeroffsetmean[pixelnr*data_roi+k])/drs_gainmean[pixelnr*data_roi+drs_calib_offset]*1907.35;
+		pix->Fill(k,sample);
+	}
+	pix->Draw();
+	canv->Modified();
+	canv->Update();
+	
+	cout << "Plotting successful..." << endl;
+	return 0;
+}
Index: fact/tools/rootmacros/calscope.C
===================================================================
--- fact/tools/rootmacros/calscope.C	(revision 12166)
+++ fact/tools/rootmacros/calscope.C	(revision 12166)
@@ -0,0 +1,84 @@
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TProfile.h>
+
+#include <stdint.h>
+#include <cstdio>
+
+#include "fits.h"
+#include "TPKplotevent.c"
+#include "FOpenDataFile.c"
+#include "FOpenCalibFile.c"
+
+int calscope(const char *name = "../raw/20110916_025.fits", const char *drsname = "../raw/20110916_024.drs.fits", size_t eventnr = 0, size_t pixelnr = 0)
+{
+//******************************************************************************
+//Read a datafile and plot the DRS-calibrated data
+//ATTENTION: only works for ROI=1024
+// (array indices of the calibration wrong otherwise)
+//Example call in ROOT:
+//root [74] .x calscope.C++("20110804_024.fits","20110804_023.drs.fits",10,1348)
+//T. Kr�henb�hl, August 2011, tpk@phys.ethz.ch
+//******************************************************************************
+
+	gROOT->SetStyle("Plain");
+	
+//-------------------------------------------
+//Open the file
+//-------------------------------------------
+	fits datafile(name);
+	if (!datafile)
+	{
+	cout << "Couldn't properly open the datafile." << endl;
+	return 1;
+	}
+	
+//-------------------------------------------
+//Get the data
+//-------------------------------------------
+	vector<int16_t> data;
+	vector<int16_t> data_offset;
+	unsigned int data_num;
+	size_t data_n;
+	UInt_t data_px;
+	UInt_t data_roi;
+	FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
+	
+//-------------------------------------------
+//Get the DRS calibration
+//-------------------------------------------
+	size_t drs_n;
+	vector<float> drs_basemean;
+	vector<float> drs_gainmean;
+	vector<float> drs_triggeroffsetmean;
+	FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
+	
+//-------------------------------------------
+//Check the sizes of the data columns
+//-------------------------------------------
+	if(drs_n!=data_n)
+	{
+		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
+		return 1;
+	}
+	
+	
+//-------------------------------------------
+//Create the title
+//-------------------------------------------
+	char title[500];
+	std::sprintf(title,"Data: %s, DRS: %s, Px %i Ev %i",name,drsname,pixelnr,eventnr);
+	
+//-------------------------------------------
+//Get the event
+//-------------------------------------------
+	cout << "--------------------- Data --------------------" << endl;
+	datafile.GetRow(eventnr);
+	cout << "Event number: " << data_num << endl;
+		
+//-------------------------------------------
+//Draw the data
+//-------------------------------------------
+	TPKplotevent(title, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr);
+	return 0;
+}
Index: fact/tools/rootmacros/calscope_persistent.C
===================================================================
--- fact/tools/rootmacros/calscope_persistent.C	(revision 12166)
+++ fact/tools/rootmacros/calscope_persistent.C	(revision 12166)
@@ -0,0 +1,153 @@
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TProfile.h>
+#include <TH2.h>
+#include <stdint.h>
+#include <cstdio>
+#include "fits.h"
+
+int calscope_persistent(const char *name, const char *drsname, size_t pixelnr, UInt_t maxevents)
+{
+//******************************************************************************
+//Read a datafile and plot the DRS-calibrated data
+//ATTENTION: only works for ROI=1024
+// (array indices of the calibration wrong otherwise)
+//Example call in ROOT:
+//root [74] .x calscope.C++("20110804_024.fits","20110804_023.drs.fits",10,1348)
+//T. Krähenbühl, August 2011, tpk@phys.ethz.ch
+//******************************************************************************
+
+	gROOT->SetStyle("Plain");
+	
+//-------------------------------------------
+//Open the files
+//-------------------------------------------
+	fits datafile(name);
+	fits drsfile(drsname);
+	
+	if (!datafile)
+	{
+	cout << "Couldn't properly open the datafile." << endl;
+	return 1;
+	}
+	
+	if (!drsfile)
+	{
+	cout << "Couldn't properly open the drsfile." << endl;
+	return 1;
+	}
+
+//-------------------------------------------
+//Print the headers
+//-------------------------------------------
+	cout << "-------------------- Data Header -------------------" << endl;
+	datafile.PrintKeys();
+	cout << "----------------- DRS Calib Header -----------------" << endl;
+	drsfile.PrintKeys();
+	cout << "------------------- Data Columns -------------------" << endl;
+	datafile.PrintColumns();
+	cout << "---------------- DRS Calib Columns -----------------" << endl;
+	drsfile.PrintColumns();
+	cout << "--------------------- Data --------------------" << endl;
+	
+//-------------------------------------------
+//Get the sizes of the data columns
+//-------------------------------------------
+	const size_t data_n = datafile.GetN("Data"); //Size of column "Data" = #Pixel x ROI
+	const size_t drs_n = drsfile.GetN("BaselineMean");
+	if(drs_n!=data_n)
+	{
+		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
+		return 1;
+	}
+	
+	const UInt_t data_roi = datafile.GetUInt("NROI"); // Value from header
+	const UInt_t data_px = datafile.GetUInt("NPIX");
+	
+//-------------------------------------------
+//Link the data to variables
+//-------------------------------------------
+	unsigned int data_num;
+	datafile.SetRefAddress("EventNum", data_num);
+	vector<int16_t> data(data_n);
+	datafile.SetVecAddress("Data", data);
+	vector<int16_t> dataOffset(data_px);
+	datafile.SetVecAddress("StartCellData", dataOffset);
+	
+	vector<float> drs_basemean(drs_n);
+	drsfile.SetVecAddress("BaselineMean", drs_basemean);
+	vector<float> drs_gainmean(drs_n);
+	drsfile.SetVecAddress("GainMean", drs_gainmean);
+	vector<float> drs_triggeroffsetmean(drs_n);
+	drsfile.SetVecAddress("TriggerOffsetMean", drs_triggeroffsetmean);
+	drsfile.GetRow(0); //Read the calibration data
+	
+//-------------------------------------------
+//Create the canvas, plot and title
+//-------------------------------------------
+	char title[500];
+	std::sprintf(title,"Data: %s, DRS: %s, Px %i",name,drsname,pixelnr);
+	TCanvas *canv = new TCanvas( "canv", "Mean values of the first event", 100, 10, 700, 500 );
+	TProfile *pix = new TProfile("pix", title, 1024, -0.5, 1023.5);
+	pix->GetXaxis()->SetTitle("DRS bin (@2 GHz)");
+	pix->GetYaxis()->SetTitle("Amplitude (mV)");
+	
+	char ptitle[500];
+	std::sprintf(ptitle,"Data: %s, DRS: %s, Px %i",name,drsname);
+	TCanvas *pcanv = new TCanvas( "pcanv", title, 800, 10, 700, 500 );
+	TH2F *pers = new TH2F("pers",ptitle,1024,-0.5,1023.5,200,-49.5,150.5);
+	pers->GetXaxis()->SetTitle("DRS bin (@2 GHz)");
+	pers->GetYaxis()->SetTitle("Amplitude (mV)");
+	
+//-------------------------------------------
+//Iterate over the events
+//-------------------------------------------
+if(maxevents==0) maxevents = datafile.GetNumRows();
+//	for (size_t i=0; i<datafile.GetNumRows(); i++)
+	for (size_t i=0; i<maxevents; i++)
+	{
+//	size_t i=eventnr; //Fix the Event
+		datafile.GetRow(i);
+		cout << "Event number: " << data_num << endl;
+		
+//-------------------------------------------
+//Iterate over the pixels
+//-------------------------------------------
+//		for (int j=0; j<data_px; j++)
+//		{
+		size_t j=pixelnr; //Fix the Pixel to a SoftID
+		
+//-------------------------------------------
+//Iterate over the slices
+//-------------------------------------------
+			for (UInt_t k=0; k<data_roi; k++)
+			{
+				UInt_t drs_calib_offset = (k+dataOffset[j])%data_roi;
+				//Example values for the various datasets (20110803_015 / 018?9)
+				//data = -1856 +- 30.76
+				//drs_basemean = -906.6 +- 14.73
+				//drs_gainmean = 1696 +- 4
+				//drs_triggeroffsetmean = -0.01 +- 1.532
+				
+				//Gain calibration to mV: 50000(DAC) / drs_gainmean(ADC) * 2500(mV) / 65536(DAC)
+				//Target value: 1.907 V für DRS-Calib files beim DAC-Wert 50000
+				float sample = (data[j*data_roi+k]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35;
+				pix->Fill(k,sample);
+				pers->Fill(k,sample);
+			}
+//		}
+	}
+
+//-------------------------------------------
+//Draw the data
+//-------------------------------------------
+	canv->cd();
+	pix->Draw();
+	canv->Modified();
+	canv->Update();
+	pcanv->cd();
+	pers->Draw("col");
+	pcanv->Modified();
+	pcanv->Update();
+	return 0;
+}
Index: fact/tools/rootmacros/factfir.C
===================================================================
--- fact/tools/rootmacros/factfir.C	(revision 12166)
+++ fact/tools/rootmacros/factfir.C	(revision 12166)
@@ -0,0 +1,21 @@
+#include <vector>
+
+#ifndef FAD_MAX_SAMPLES
+	#warning FAD_MAX_SAMPLES not defined. defining: FAD_MAX_SAMPLES to 1024
+	#define FAD_MAX_SAMPLES 1024
+#endif
+
+// source vector is
+// float Ameas[FAD_MAX_SAMPLES];
+void factfir(double b, vector<double> &a, int k, vector<float> &source, vector<float> &dest){
+	//dest.clear();
+	
+	for (int slice =0; slice < FAD_MAX_SAMPLES; slice++) {
+		float currentval = 0;
+		
+		for (int i=0; i < k ; i++){
+			currentval += a[i] * source[ (slice - i + FAD_MAX_SAMPLES) % FAD_MAX_SAMPLES ];
+		}	
+		dest[slice] = currentval / b;	
+	}
+}
Index: fact/tools/rootmacros/fana.C
===================================================================
--- fact/tools/rootmacros/fana.C	(revision 12166)
+++ fact/tools/rootmacros/fana.C	(revision 12166)
@@ -0,0 +1,444 @@
+#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 <TMath.h>
+
+#include <stdint.h>
+#include <cstdio>
+
+#include "fits.h"
+//#include "TPKplotevent.c"
+//#include "FOpenDataFile.c"
+#include "FOpenCalibFile.c"
+
+#include "zerosearch.C"
+
+#define NPIX  1440
+#define NCELL 1024
+
+// data access and treatment
+#define FAD_MAX_SAMPLES 1024
+
+vector<int16_t> data;
+vector<int16_t> data_offset;
+
+unsigned int data_num;
+
+size_t data_n;
+UInt_t data_px;
+UInt_t data_roi;
+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> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 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> Vcfd_slide(FAD_MAX_SAMPLES);  // CDF result sliding averaged again
+
+#include "factfir.C"
+
+float getValue( int, int );
+void computeN1mean( int );
+void removeSpikes( int );
+
+// histograms
+const int Ntypes = 6;
+const unsigned int  tAmeas = 0, 
+	tN1mean = 1, 
+	tVcorr = 2,  
+	tVslide = 3, 
+	tVcfd = 4, 
+	tVcfd_slide = 5;
+
+
+TH1F* h;
+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.);
+
+void BookHistos( int );
+
+// Create a canvas
+TCanvas* CW;
+
+int spikeDebug = 0;
+
+
+int fana( const char *datafilename = "/media/daten_platte/FACT/data/20111009_013.fits", 
+	const char *drsfilename = "/media/daten_platte/FACT/data/20111009_009.drs.fits", 
+	int pixelnr = 0, 
+	int firstevent = 0, 
+	int nevents = -1 )
+{
+// 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;
+
+	gROOT->SetStyle("Plain");
+	
+//-------------------------------------------
+// Open the file
+//-------------------------------------------
+	fits datafile( 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);
+
+    // compare the number of events in the data file with the nevents the
+    // the user would like to read. nevents = -1 means: read all
+    if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;
+
+	
+//-------------------------------------------
+//Get the DRS calibration
+//-------------------------------------------
+
+	FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
+	
+//-------------------------------------------
+//Check the sizes of the data columns
+//-------------------------------------------
+	if(drs_n!=data_n)
+	{
+		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
+		return 1;
+	}
+// Book the histograms
+    BookHistos( data_roi );
+
+// Loop over events
+	cout << "--------------------- Data --------------------" << endl;
+
+    float value;
+		float integral =0;
+
+		TH1F * sp = new TH1F("spektrum", "test of Stepktrum", 256, -0.5, 63.5);
+
+    for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
+
+	    datafile.GetRow( ev );
+	    if (ev % 50 ==0)
+				cout << "Event number: " << data_num << endl;
+
+        for ( int pix = 0; pix < 1440; pix++ ){
+            hStartCell->Fill( pix, data_offset[pix] );
+        }
+
+        for ( unsigned int sl = 0; sl < data_roi; sl++){
+            value = getValue(sl, pixelnr);
+            //printf("value = %f\n", value);
+            Ameas[ sl ] = value;
+            h[tAmeas].SetBinContent(sl, value);
+            
+        }
+
+        computeN1mean( data_roi );
+        removeSpikes( data_roi );
+
+        for ( unsigned int sl = 0; sl < data_roi; sl++){
+            hPixelCellData.Fill(sl, Vcorr[ sl ]);
+        }   
+
+				// calculate sliding mean using FIR filter function
+				factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
+				// calculate CDF
+				factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
+				// calculate i2nd sliding mean
+				factfir(b_slide , a_slide, k_slide, Vcfd, Vcfd_slide);
+        
+				for ( unsigned int sl = 0; sl < data_roi; sl++){
+            h[tVcfd].SetBinContent(sl, Vcfd[sl]);
+            h[tVslide].SetBinContent(sl, Vslide[sl]);
+            h[tVcfd_slide].SetBinContent(sl, Vcfd_slide[sl]);
+        } 
+
+
+				vector<int> * zeros = zerosearch(Vcfd_slide,-1,10,20);
+				if (zeros->size() ==0)
+					continue;
+
+				// check value of Vside at zero position	
+				for ( int i=0; i<zeros->size(); i++){
+					cout << zeros->at(i) << ":\t" << Vslide[zeros->at(i)]<<endl;
+					sp->Fill(Vslide[zeros->at(i)]);
+				}
+
+/*
+
+				vector<int> goodzeros;
+				if (zeros->size() > 1){
+				for ( int i=0; i<zeros->size(); i++){
+					if (i==0) {
+						if ( abs(zeros->at(i) - zeros->at(i+1)) > 90 )
+							goodzeros.push_back(zeros->at(i));
+					} else if ( i==zeros->size()-1){
+						if (abs(zeros->at(i) - zeros->at(i-1)) > 90)
+							goodzeros.push_back(zeros->at(i));
+					} else {
+						if (abs(zeros->at(i) - zeros->at(i-1)) > 90 && 
+							abs(zeros->at(i) - zeros->at(i+1)) > 90 ) 
+								goodzeros.push_back(zeros->at(i));
+					}
+				}
+				} else {
+					goodzeros.push_back(zeros->at(0));
+				}
+				
+				// compute integrals around the good zeros.
+				// from -30 to +80 around it.
+				for (int i=0; i<goodzeros.size(); i++){
+//					cout << goodzeros[i];
+					integral =0;
+					for (int sl = goodzeros[i]-30; sl < goodzeros[i]+80; sl++){
+						if ( Ameas[sl] > integral)
+							integral = Ameas[sl];
+					} 
+//					cout << "\t-->" << integral << endl;
+					sp->Fill(integral);	
+				}      
+
+*/			 
+ 
+        
+				if ( spikeDebug ){
+
+            CW->cd( tAmeas + 1);
+            gPad->SetGrid();
+            h[tAmeas].Draw();
+
+            CW->cd( tN1mean + 1);
+            gPad->SetGrid();
+            h[tN1mean].Draw();
+
+            CW->cd( tVcorr + 1);
+            gPad->SetGrid();
+            h[tVcorr].Draw();
+
+        		CW->cd( tVslide + 1);
+		        gPad->SetGrid();
+    		    h[tVslide].Draw();
+        
+        		CW->cd( tVcfd + 1);
+		        gPad->SetGrid();
+    		    h[tVcfd].Draw();
+						TLine zeroline(0,0,1024,0);
+						zeroline.SetLineColor(kBlue);
+						zeroline.Draw();
+        
+        		CW->cd( tVcfd_slide + 1);
+		        gPad->SetGrid();
+    		    h[tVcfd_slide].Draw();
+						zeroline.Draw();
+				    
+						CW->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") break;
+        }
+
+
+	}
+    
+    TCanvas * cSpektrum = new TCanvas();
+    cSpektrum->cd();
+    sp->Draw();
+    
+
+		TCanvas * cStartCell = new TCanvas();
+    cStartCell->cd();
+    hStartCell->Draw();
+    hPixelCellData.Draw();
+		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++ ) h[ tVcorr ].SetBinContent( i, Vcorr[i] );
+}
+
+void computeN1mean( int Samples ){
+
+// 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
+
+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 * data_roi;
+    slice_pt = pixel_pt + slice;
+    drs_cal_offset = ( slice + data_offset[ pixel ] )%data_roi;
+    cal_pt    = pixel_pt + drs_cal_offset;
+
+    vraw = data[ slice_pt ] * dconv;
+    vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
+
+    return( vcal );
+}
+
+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, Ntypes);
+
+    hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
+
+}
+int FOpenDataFile(fits &datafile){
+
+/*  cout << "-------------------- Data Header -------------------" << endl;
+    datafile.PrintKeys();
+    cout << "------------------- Data Columns -------------------" << endl;
+    datafile.PrintColumns();
+    */
+
+    //Get the size of the data column
+    data_roi = datafile.GetUInt("NROI");  // Value from header
+    data_px = datafile.GetUInt("NPIX");
+    data_n = datafile.GetN("Data");         //Size of column "Data" = #Pixel x ROI
+    
+    //Set the sizes of the data vectors
+    data.resize(data_n,0);
+    data_offset.resize(data_px,0);
+    
+    //Link the data to variables
+    datafile.SetRefAddress("EventNum", data_num);
+    datafile.SetVecAddress("Data", data);
+    datafile.SetVecAddress("StartCellData", data_offset);
+    datafile.GetRow(0);
+    
+    cout << "Opening data file successful..." << endl;
+
+    return datafile.GetNumRows() ;
+}
+
Index: fact/tools/rootmacros/fir.h
===================================================================
--- fact/tools/rootmacros/fir.h	(revision 12166)
+++ fact/tools/rootmacros/fir.h	(revision 12166)
@@ -0,0 +1,8 @@
+#ifndef __FIR_H
+#define __FIR_H
+
+	#define FAD_MAX_SAMPLES 1024
+//	extern float Vmeas[FAD_MAX_SAMPLES];
+//	extern float Vfir[FAD_MAX_SAMPLES];
+
+#endif
Index: fact/tools/rootmacros/fits.h
===================================================================
--- fact/tools/rootmacros/fits.h	(revision 12166)
+++ fact/tools/rootmacros/fits.h	(revision 12166)
@@ -0,0 +1,657 @@
+#ifndef MARS_fits
+#define MARS_fits
+
+#include <stdint.h>
+
+#include <map>
+#include <string>
+#include <sstream>
+#include <algorithm>
+
+#ifdef __EXCEPTIONS
+#include <stdexcept>
+#endif
+
+#ifdef __CINT__
+#define off_t size_t
+#endif
+
+#ifndef __MARS__
+#include <vector>
+#include <iomanip>
+#include <iostream>
+#define gLog cerr
+#define err ""
+#define all ""
+#else
+#include "MLog.h"
+#include "MLogManip.h"
+#endif
+
+#ifdef HAVE_ZLIB
+#include "izstream.h"
+#else
+#include <fstream>
+#define izstream ifstream
+#warning Support for zipped FITS files disabled.
+#endif
+
+namespace std
+{
+
+class fits : public izstream
+{
+public:
+    struct Entry
+    {
+        char   type;
+        string value;
+        string comment;
+
+        template<typename T>
+            T Get() const
+        {
+            T t;
+
+            istringstream str(value);
+            str >> t;
+
+            return t;
+        }
+    };
+
+    struct Table
+    {
+        off_t offset;
+
+        string name;
+        size_t bytes_per_row;
+        size_t num_rows;
+        size_t num_cols;
+
+        struct Column
+        {
+            size_t offset;
+            size_t num;
+            size_t size;
+            char   type;
+        };
+
+        typedef map<string, Entry>  Keys;
+        typedef map<string, Column> Columns;
+
+        Columns cols;
+        Keys    keys;
+
+        string Trim(const string &str, char c=' ') const
+        {
+            // Trim Both leading and trailing spaces
+            const size_t pstart = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces
+            const size_t pend   = str.find_last_not_of(c);  // Find the first character position from reverse af
+
+            // if all spaces or empty return an empty string
+            if (string::npos==pstart || string::npos==pend)
+                return string();
+
+            return str.substr(pstart, pend-pstart+1);
+        }
+
+        bool Check(const string &key, char type, const string &value="") const
+        {
+            const Keys::const_iterator it = keys.find(key);
+            if (it==keys.end())
+            {
+                ostringstream str;
+                str << "Key '" << key << "' not found.";
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << err << "ERROR - " << str.str() << endl;
+                return false;
+#endif
+            }
+
+            if (it->second.type!=type)
+            {
+                ostringstream str;
+                str << "Wrong type for key '" << key << "': expected " << type << ", found " << it->second.type << ".";
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << err << "ERROR - " << str.str() << endl;
+                return false;
+#endif
+            }
+
+            if (!value.empty() && it->second.value!=value)
+            {
+                ostringstream str;
+                str << "Wrong value for key '" << key << "': expected " << value << ", found " << it->second.value << ".";
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << err << "ERROR - " << str.str() << endl;
+                return false;
+#endif
+            }
+
+            return true;
+        }
+
+        Keys ParseBlock(const vector<string> &vec) const
+        {
+            map<string,Entry> rc;
+
+            for (unsigned int i=0; i<vec.size(); i++)
+            {
+                const string key = Trim(vec[i].substr(0,8));
+                // Keywords without a value, like COMMENT / HISTORY
+                if (vec[i].substr(8,2)!="= ")
+                    continue;
+
+                char type = 0;
+
+                string com;
+                string val = Trim(vec[i].substr(10));
+                if (val[0]=='\'')
+                {
+                    // First skip all '' in the string
+                    size_t p = 1;
+                    while (1)
+                    {
+                        const size_t pp = val.find_first_of('\'', p);
+                        if (pp==string::npos)
+                            break;
+
+                        p = val[pp+1]=='\'' ? pp+2 : pp+1;
+                    }
+
+                    // Now find the comment
+                    const size_t ppp = val.find_first_of('/', p);
+
+                    // Set value, comment and type
+                    com  = ppp==string::npos ? "" : Trim(val.substr(ppp+1));
+                    val  = Trim(val.substr(1, p-2));
+                    type = 'T';
+                }
+                else
+                {
+                    const size_t p = val.find_first_of('/');
+
+                    com = Trim(val.substr(p+2));
+                    val = Trim(val.substr(0, p));
+
+                    if (val.empty() || val.find_first_of('T')!=string::npos || val.find_first_of('F')!=string::npos)
+                        type = 'B';
+                    else
+                        type = val.find_last_of('.')==string::npos ? 'I' : 'F';
+                }
+
+                const Entry e = { type, val, com };
+                rc[key] = e;
+            }
+
+            return rc;
+        }
+
+        Table() : offset(0) { }
+        Table(const vector<string> &vec, off_t off) :
+            offset(off), keys(ParseBlock(vec))
+        {
+            if (!Check("XTENSION", 'T', "BINTABLE") ||
+                !Check("NAXIS",    'I', "2")        ||
+                !Check("BITPIX",   'I', "8")        ||
+                !Check("PCOUNT",   'I', "0")        ||
+                !Check("GCOUNT",   'I', "1")        ||
+                !Check("EXTNAME",  'T')             ||
+                !Check("NAXIS1",   'I')             ||
+                !Check("NAXIS2",   'I')             ||
+                !Check("TFIELDS",  'I'))
+                return;
+
+            bytes_per_row = Get<size_t>("NAXIS1");
+            num_rows      = Get<size_t>("NAXIS2");
+            num_cols      = Get<size_t>("TFIELDS");
+
+            size_t bytes = 0;
+            for (size_t i=1; i<=num_cols; i++)
+            {
+                ostringstream num;
+                num << i;
+
+                if (!Check("TTYPE"+num.str(), 'T') ||
+                    !Check("TFORM"+num.str(), 'T'))
+                    return;
+
+                const string id  = Get<string>("TTYPE"+num.str());
+                const string fmt = Get<string>("TFORM"+num.str());
+
+                istringstream sin(fmt);
+                int n = 0;
+                sin >> n;
+                if (!sin)
+                    n = 1;
+
+                const char type = fmt[fmt.length()-1];
+
+                size_t size = 0;
+                switch (type)
+                {
+                    // We could use negative values to mark floats
+                    // otheriwse we could just cast them to int64_t?
+                case 'L': size = 1; break; // logical
+                // case 'X': size = n; break; // bits (n=number of bytes needed to contain all bits)
+                case 'B': size = 1; break; // byte
+                case 'I': size = 2; break; // short
+                case 'J': size = 4; break; // int
+                case 'K': size = 8; break; // long long
+                case 'E': size = 4; break; // float
+                case 'D': size = 8; break; // double
+                // case 'C': size =  8; break; // complex float
+                // case 'M': size = 16; break; // complex double
+                // case 'P': size =  8; break; // array descriptor (32bit)
+                // case 'Q': size = 16; break; // array descriptor (64bit)
+                default:
+                    {
+                        ostringstream str;
+                        str << "FITS format TFORM='" << fmt << "' not yet supported.";
+#ifdef __EXCEPTIONS
+                        throw runtime_error(str.str());
+#else
+                        gLog << err << "ERROR - " << str.str() << endl;
+                        return;
+#endif
+                    }
+                }
+
+                const Table::Column col = { bytes, n, size, type };
+
+                cols[id] = col;
+                bytes  += n*size;
+            }
+
+            if (bytes!=bytes_per_row)
+            {
+#ifdef __EXCEPTIONS
+                throw runtime_error("Column size mismatch");
+#else
+                gLog << err << "ERROR - Column size mismatch" << endl;
+                return;
+#endif
+            }
+
+            name = Get<string>("EXTNAME");
+        }
+
+        void PrintKeys() const
+        {
+            for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++)
+                gLog << all << setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << endl;
+        }
+
+        void PrintColumns() const
+        {
+            typedef map<pair<size_t, string>, Column> Sorted;
+
+            Sorted sorted;
+
+            for (Columns::const_iterator it=cols.begin(); it!=cols.end(); it++)
+                sorted[make_pair(it->second.offset, it->first)] = it->second;
+
+            for (Sorted::const_iterator it=sorted.begin(); it!=sorted.end(); it++)
+                gLog << all << setw(6) << it->second.offset << "| " << it->first.second << " [" << it->second.type << '=' << it->second.num << 'x' << it->second.size << "byte]" << endl;
+        }
+
+        operator bool() const { return !name.empty(); }
+
+        bool HasKey(const string &key) const
+        {
+            return keys.find(key)!=keys.end();
+        }
+
+            // Values of keys are always signed
+        template<typename T>
+            T Get(const string &key) const
+        {
+            const map<string,Entry>::const_iterator it = keys.find(key);
+            if (it==keys.end())
+            {
+                ostringstream str;
+                str << "Key '" << key << "' not found." << endl;
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << err << "ERROR - " << str.str() << endl;
+                return 0;
+#endif
+            }
+            return it->second.Get<T>();
+        }
+
+        size_t GetN(const string &key) const
+        {
+            const Columns::const_iterator it = cols.find(key);
+            return it==cols.end() ? 0 : it->second.num;
+        }
+    };
+
+private:
+    Table fTable;
+
+    typedef pair<void*, Table::Column> Address;
+    typedef vector<Address> Addresses;
+    //map<void*, Table::Column> fAddresses;
+    Addresses fAddresses;
+
+    vector<char> fBufferRow;
+    vector<char> fBufferDat;
+
+    size_t fRow;
+
+    vector<string> ReadBlock(vector<string> &vec)
+    {
+        bool endtag = false;
+        for (int i=0; i<36; i++)
+        {
+            char c[81];
+            c[80] = 0;
+            read(c, 80);
+            if (!good())
+                break;
+
+            if (c[0]==0)
+                return vector<string>();
+
+            string str(c);
+
+//            if (!str.empty())
+//                cout << setw(2) << i << "|" << str << "|" << (endtag?'-':'+') << endl;
+
+            if (str=="END                                                                             ")
+            {
+                endtag = true;
+
+                // Make sure that no empty vector is returned
+                if (vec.size()%36==0)
+                    vec.push_back(string("END     = '' / "));
+            }
+
+            if (endtag)
+                continue;
+
+            vec.push_back(str);
+        }
+
+        return vec;
+    }
+
+    string Compile(const string &key, int16_t i=-1)
+    {
+        if (i<0)
+            return key;
+
+        ostringstream str;
+        str << key << i;
+        return str.str();
+    }
+
+public:
+    fits(const string &fname) : izstream(fname.c_str())
+    {
+        char simple[10];
+        read(simple, 10);
+        if (!good())
+            return;
+
+        if (memcmp(simple, "SIMPLE  = ", 10))
+        {
+            clear(rdstate()|ios::badbit);
+#ifdef __EXCEPTIONS
+            throw runtime_error("File is not a FITS file.");
+#else
+            gLog << err << "ERROR - File is not a FITS file." << endl;
+            return;
+#endif
+        }
+
+        seekg(0);
+
+        while (good())
+        {
+            vector<string> block;
+            while (1)
+            {
+                // FIXME: Set limit on memory consumption
+                ReadBlock(block);
+                if (!good())
+                {
+                    clear(rdstate()|ios::badbit);
+#ifdef __EXCEPTIONS
+                    throw runtime_error("FITS file corrupted.");
+#else
+                    gLog << err << "ERROR - FITS file corrupted." << endl;
+                    return;
+#endif
+                }
+
+                if (block.size()%36)
+                    break;
+            }
+
+            if (block.size()==0)
+                break;
+
+            if (block[0].substr(0, 9)=="SIMPLE  =")
+                continue;
+
+            if (block[0].substr(0, 9)=="XTENSION=")
+            {
+                // FIXME: Check for table name
+
+                fTable = Table(block, tellg());
+                fRow   = (size_t)-1;
+
+                //fTable.PrintKeys();
+
+                if (!fTable)
+                {
+                    clear(rdstate()|ios::badbit);
+                    return;
+                }
+
+                fBufferRow.resize(fTable.bytes_per_row);
+                fBufferDat.resize(fTable.bytes_per_row);
+
+                /*
+                 // Next table should start at:
+                 const size_t size = fTable.bytes_per_row*fTable.num_rows;
+                 const size_t blks = size/(36*80);
+                 const size_t rest = size%(36*80);
+
+                seekg((blks+(rest>0?1:0))*(36*80), ios::cur);
+                if (!good())
+                     gLog << err << "File seems to be incomplete (less data than expected from header)." << endl;
+
+                fRow   = fTable.num_rows;
+                 */
+
+                break;
+            }
+        }
+    }
+
+    void ReadRow(size_t row)
+    {
+        // if (row!=fRow+1) // Fast seeking is ensured by izstream
+        seekg(fTable.offset+row*fTable.bytes_per_row);
+
+        fRow = row;
+
+        read(fBufferRow.data(), fBufferRow.size());
+        //fin.clear(fin.rdstate()&~ios::eofbit);
+    }
+
+    template<size_t N>
+        void revcpy(char *dest, const char *src, int num)
+    {
+        const char *pend = src + num*N;
+        for (const char *ptr = src; ptr<pend; ptr+=N, dest+=N)
+            reverse_copy(ptr, ptr+N, dest);
+    }
+
+    bool GetRow(size_t row)
+    {
+        ReadRow(row);
+        if (!good())
+            return good();
+
+        for (Addresses::const_iterator it=fAddresses.begin(); it!=fAddresses.end(); it++)
+        {
+            const Table::Column &c = it->second;
+
+            const char *src = fBufferRow.data() + c.offset;
+            char *dest = reinterpret_cast<char*>(it->first);
+
+            // Let the compiler do some optimization by
+            // knowing the we only have 1, 2, 4 and 8
+            switch (c.size)
+            {
+            case 1: memcpy   (dest, src, c.num*c.size); break;
+            case 2: revcpy<2>(dest, src, c.num);        break;
+            case 4: revcpy<4>(dest, src, c.num);        break;
+            case 8: revcpy<8>(dest, src, c.num);        break;
+            }
+        }
+
+        return good();
+    }
+
+    bool GetNextRow()
+    {
+        return GetRow(fRow+1);
+    }
+
+    bool SkipNextRow()
+    {
+        seekg(fTable.offset+(++fRow)*fTable.bytes_per_row);
+        return good();
+    }
+
+    static bool Compare(const Address &p1, const Address &p2)
+    {
+        return p1.first>p2.first;
+    }
+
+    template<typename T>
+    bool SetPtrAddress(const string &name, T *ptr, size_t cnt)
+    {
+        if (fTable.cols.count(name)==0)
+        {
+            ostringstream str;
+            str <<"SetPtrAddress('" << name << "') - Column not found." << endl;
+#ifdef __EXCEPTIONS
+            throw runtime_error(str.str());
+#else
+            gLog << err << "ERROR - " << str.str() << endl;
+            return false;
+#endif
+        }
+
+        if (sizeof(T)!=fTable.cols[name].size)
+        {
+            ostringstream str;
+            str << "SetPtrAddress('" << name << "') - Element size mismatch: expected "
+                << fTable.cols[name].size << " from header, got " << sizeof(T) << endl;
+#ifdef __EXCEPTIONS
+            throw runtime_error(str.str());
+#else
+            gLog << err << "ERROR - " << str.str() << endl;
+            return false;
+#endif
+        }
+
+        if (cnt!=fTable.cols[name].num)
+        {
+            ostringstream str;
+            str << "SetPtrAddress('" << name << "') - Element count mismatch: expected "
+                << fTable.cols[name].num << " from header, got " << cnt << endl;
+#ifdef __EXCEPTIONS
+            throw runtime_error(str.str());
+#else
+            gLog << err << "ERROR - " << str.str() << endl;
+            return false;
+#endif
+        }
+
+        // if (fAddresses.count(ptr)>0)
+        //     gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
+
+        //fAddresses[ptr] = fTable.cols[name];
+        fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
+        sort(fAddresses.begin(), fAddresses.end(), Compare);
+        return true;
+    }
+
+    template<class T>
+    bool SetRefAddress(const string &name, T &ptr)
+    {
+        return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T));
+    }
+
+    template<typename T>
+    bool SetVecAddress(const string &name, vector<T> &vec)
+    {
+        return SetPtrAddress(name, vec.data(), vec.size());
+    }
+
+    template<typename T>
+        T Get(const string &key) const
+    {
+        return fTable.Get<T>(key);
+    }
+
+    bool SetPtrAddress(const string &name, void *ptr)
+    {
+        if (fTable.cols.count(name)==0)
+        {
+            ostringstream str;
+            str <<"SetPtrAddress('" << name << "') - Column not found." << endl;
+#ifdef __EXCEPTIONS
+            throw runtime_error(str.str());
+#else
+            gLog << err << "ERROR - " << str.str() << endl;
+            return false;
+#endif
+        }
+
+        // if (fAddresses.count(ptr)>0)
+        //     gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
+
+        //fAddresses[ptr] = fTable.cols[name];
+        fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
+        sort(fAddresses.begin(), fAddresses.end(), Compare);
+        return true;
+    }
+
+    bool     HasKey(const string &key) const { return fTable.HasKey(key); }
+    int64_t  GetInt(const string &key) const { return fTable.Get<int64_t>(key); }
+    uint64_t GetUInt(const string &key) const { return fTable.Get<uint64_t>(key); }
+    double   GetFloat(const string &key) const { return fTable.Get<double>(key); }
+    string   GetStr(const string &key) const { return fTable.Get<string>(key); }
+
+    size_t GetN(const string &key) const
+    {
+        return fTable.GetN(key);
+    }
+
+    size_t GetNumRows() const { return fTable.num_rows; }
+    size_t GetRow() const { return fRow; }
+
+    operator bool() const { return fTable && fTable.offset!=0; }
+
+    void PrintKeys() const { fTable.PrintKeys(); }
+    void PrintColumns() const { fTable.PrintColumns(); }
+};
+
+};
+#endif
Index: fact/tools/rootmacros/gainanalysis.C
===================================================================
--- fact/tools/rootmacros/gainanalysis.C	(revision 12166)
+++ fact/tools/rootmacros/gainanalysis.C	(revision 12166)
@@ -0,0 +1,273 @@
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TProfile.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TGraph.h>
+
+#include <stdint.h>
+#include <cstdio>
+
+#include "fits.h"
+#include "FOpenDataFile.c"
+#include "FOpenCalibFile.c"
+#include "FPedestal.c"
+#include "FOscilloscopeAllPx.c"
+
+#define n_peaks 2
+
+int gainanalysis(const char *name, const char *drsname, size_t pixelnr)
+{
+//******************************************************************************
+//Read a datafile and plot the DRS-calibrated data
+//ATTENTION: only works for ROI=1024
+// (array indices of the calibration wrong otherwise)
+//******************************************************************************
+
+	gROOT->SetStyle("Plain");
+	
+//-------------------------------------------
+//Define the data variables
+//-------------------------------------------
+	vector<int16_t> data;
+	vector<int16_t> data_offset;
+	unsigned int data_num;
+	size_t data_n;
+	UInt_t data_px;
+	UInt_t data_roi;
+	
+//-------------------------------------------
+//Get the DRS calibration
+//-------------------------------------------
+	size_t drs_n;
+	vector<float> drs_basemean;
+	vector<float> drs_gainmean;
+	vector<float> drs_triggeroffsetmean;
+	FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
+	
+//-------------------------------------------
+//Check the sizes of the data columns
+//-------------------------------------------
+//	if(drs_n!=data_n)
+//	{
+//		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
+//		return 1;
+//	}
+	
+//-------------------------------------------
+//Open the file
+//-------------------------------------------
+	fits datafile(name);
+	if (!datafile)
+	{
+	cout << "Couldn't properly open the datafile." << endl;
+	return 1;
+	}
+	FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
+	
+//-------------------------------------------
+//Create the canvas, plot and title
+//-------------------------------------------
+	char title_spect[500];
+	char name_spect[50];
+	std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",pixelnr, name,drsname);
+	TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
+	TH1F* spectrum[data_px];
+	for(int i=0; i<data_px; i++) {
+		std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",i,name,drsname);
+		std::sprintf(name_spect,"spect%i",i);
+		spectrum[i] = new TH1F(name_spect,title_spect,300,0,600);
+		spectrum[i]->GetXaxis()->SetTitle("Pulse size (a.u.)");
+		spectrum[i]->GetYaxis()->SetTitle("Counts");
+	}
+	
+	TCanvas *canv_gain= new TCanvas( "canv_gain", "Gain in arbitrary units", 50, 500, 1800, 800 );
+//-------------------------------------------
+//Define the fit
+//-------------------------------------------
+	int fit_start = 20;
+	int fit_end = 200;
+	TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp((-x-[7])/[8])",fit_start,fit_end);
+	Double_t par_init[5+2*n_peaks] = {70,90,100,15,10,15,10,30,30};
+	Double_t parread_tmp[5+2*n_peaks];
+	Double_t parread[data_px][5+2*n_peaks];
+//	TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp(-(x-[12])*(x-[12])/(2*[11]*[11]))",fit_start,fit_end);
+//	Double_t par_init[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
+	f_peaks->SetParameters(par_init);
+	f_peaks->SetLineColor(2);
+	
+	TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end);
+	f_pedestal->SetLineColor(2);
+	f_pedestal->SetLineWidth(1);
+	
+	TF1* f_peak[n_peaks];
+	for(int i=0; i<n_peaks; i++) {
+		f_peak[i] = new TF1("peak","gaus",fit_start,fit_end);
+		f_peak[i]->SetLineColor(2);
+		f_peak[i]->SetLineWidth(1);
+	}
+	
+//-------------------------------------------
+//Find the baseline
+//-------------------------------------------
+	float pedestal_mean, pedestal_rms;
+//	FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms);
+	//Return values
+	//Value of maximal probability: -2.225 +- 4.3144
+	pedestal_mean = -2.225;
+	pedestal_rms = 4.3144;
+	
+//-------------------------------------------
+//Oscilloscope
+//-------------------------------------------
+//	float threshold = pedestal_mean+pedestal_rms;
+	float threshold = 2.5;
+
+//***********************
+	FOscilloscopeAllPx(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
+	
+	canv_spect->cd();
+	spectrum[0]->Draw();
+	canv_spect->Modified();
+	canv_spect->Update();
+//***********************
+//	fits datafile2("20110804_025.fits");
+//	if (!datafile2)
+//	{
+//	cout << "Couldn't properly open the datafile2." << endl;
+//	return 1;
+//	}
+//	FOpenDataFile(datafile2, data, data_offset, data_num, data_n, data_roi, data_px);
+//	FOscilloscopeAllPx(datafile2, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);//	
+//	canv_spect->cd();
+//	spectrum[0]->Draw();
+//	canv_spect->Modified();
+//	canv_spect->Update();
+////***********************
+//	fits datafile3("20110804_026.fits");
+//	if (!datafile3)
+//	{
+//	cout << "Couldn't properly open the datafile3." << endl;
+//	return 1;
+//	}
+//	FOpenDataFile(datafile3, data, data_offset, data_num, data_n, data_roi, data_px);
+//	FOscilloscopeAllPx(datafile3, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
+//	
+//	canv_spect->cd();
+//	spectrum[0]->Draw();
+//	canv_spect->Modified();
+//	canv_spect->Update();
+////***********************
+//	fits datafile4("20110804_027.fits");
+//	if (!datafile4)
+//	{
+//	cout << "Couldn't properly open the datafile4." << endl;
+//	return 1;
+//	}
+//	FOpenDataFile(datafile4, data, data_offset, data_num, data_n, data_roi, data_px);
+//	FOscilloscopeAllPx(datafile4, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
+//	
+//	canv_spect->cd();
+//	spectrum[0]->Draw();
+//	canv_spect->Modified();
+//	canv_spect->Update();
+//***********************
+	Double_t parread_dist[data_px];
+	Double_t parread_x[data_px];
+	for(int i=0; i<data_px; i++) {
+		f_peaks->SetParameters(par_init);
+		if(!(spectrum[i]->Fit(f_peaks,"NR"))) {
+//			std::cout << "Fit successful." << std::endl;
+			f_peaks->GetParameters(&parread_tmp[0]);
+			for(int j=0; j<5+2*n_peaks; j++) {
+				parread[i][j] = parread_tmp[j];
+			}
+			parread_dist[i] = parread_tmp[1];
+			parread_x[i] = i;
+		}
+		else {
+			for(int j=0; j<5+2*n_peaks; j++) {
+				parread[i][j] = 0;
+			}
+			parread_dist[i] = 0;
+			parread_x[i] = i;
+		}
+	}
+	
+	canv_gain->cd();
+	TGraph* plot_gain = new TGraph(data_px,parread_x,parread_dist);
+	plot_gain->GetYaxis()->SetRangeUser(0,150);
+	plot_gain->Draw("A*");
+	canv_gain->Modified();
+	canv_gain->Update();
+	
+	canv_spect->cd();
+	canv_spect->SetLogy();
+	
+//-------------------------------------------
+//Draw the data where the gain is out of limits
+//-------------------------------------------
+	char temp = 'x';
+	std::cout << "Plot the spectra out of limits: Enter for next, 'a' to abort." << std::endl;
+	for(int i=0; i<data_px; i++) {
+		if((parread_dist[i]<80)||(parread_dist[i]>105)) {
+			std::cout << "Gain out of limits: Px " << i << " at " << parread_dist[i] << std::endl;
+			spectrum[i]->Draw();
+			
+			f_peaks->SetParameters(&parread[i][0]);
+			f_peak[0]->SetParameters(parread[i][2],parread[i][0],parread[i][3]);
+			f_peak[1]->SetParameters(parread[i][4],parread[i][0]+parread[i][1],parread[i][5]);
+//			f_peak[2]->SetParameters(parread[i][6],parread[i][0]+2*parread[i][1],parread[i][7]);
+//			f_peak[3]->SetParameters(parread[i][8],parread[i][0]+3*parread[i][1],parread[i][9]);
+//			f_peak[4]->SetParameters(parread[i][10],parread[i][12],parread[i][11]);
+			f_pedestal->SetParameters(parread[i][6],parread[i][7],parread[i][8]);
+			
+			f_peaks->Draw("same");
+			for(int i=0; i<n_peaks; i++) {
+				f_peak[i]->Draw("same");
+			}
+			f_pedestal->Draw("same");
+			
+			canv_spect->Modified();
+			canv_spect->Update();
+			temp='x';
+			while ((temp!='\n') && (temp!='a'))
+			    cin.get(temp);
+			if(temp=='a') break;
+		}
+	}
+//-------------------------------------------
+//Draw the data
+//-------------------------------------------
+	canv_spect->cd();
+	Int_t plotspect_nr;
+	std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
+	std::cin >> plotspect_nr;
+	while(plotspect_nr>=0) {
+		plotspect_nr = plotspect_nr % data_px;
+		spectrum[plotspect_nr]->Draw();
+		
+		f_peaks->SetParameters(&parread[plotspect_nr][0]);
+		f_peak[0]->SetParameters(parread[plotspect_nr][2],parread[plotspect_nr][0],parread[plotspect_nr][3]);
+		f_peak[1]->SetParameters(parread[plotspect_nr][4],parread[plotspect_nr][0]+parread[plotspect_nr][1],parread[plotspect_nr][5]);
+//		f_peak[2]->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][0]+2*parread[plotspect_nr][1],parread[plotspect_nr][7]);
+//		f_peak[3]->SetParameters(parread[plotspect_nr][8],parread[plotspect_nr][0]+3*parread[plotspect_nr][1],parread[plotspect_nr][9]);
+//		f_peak[4]->SetParameters(parread[plotspect_nr][10],parread[plotspect_nr][12],parread[plotspect_nr][11]);
+		f_pedestal->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][7],parread[plotspect_nr][8]);
+		
+		f_peaks->Draw("same");
+		for(int i=0; i<n_peaks; i++) {
+			f_peak[i]->Draw("same");
+		}
+		f_pedestal->Draw("same");
+		
+		canv_spect->Modified();
+		canv_spect->Update();
+		
+		std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
+		std::cin >> plotspect_nr;
+	}
+	
+	return 0;
+}
Index: fact/tools/rootmacros/peakfinderscope.C
===================================================================
--- fact/tools/rootmacros/peakfinderscope.C	(revision 12166)
+++ fact/tools/rootmacros/peakfinderscope.C	(revision 12166)
@@ -0,0 +1,231 @@
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TProfile.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TF1.h>
+
+#include <stdint.h>
+#include <cstdio>
+
+#include "fits.h"
+#include "FOpenDataFile.c"
+#include "FOpenCalibFile.c"
+#include "FPedestal.c"
+#include "FOscilloscope.c"
+
+#define n_peaks 4
+
+int peakfinderscope(const char *name, const char *drsname, size_t eventnr, size_t pixelnr)
+{
+//******************************************************************************
+//Read a datafile and plot the DRS-calibrated data
+//ATTENTION: only works for ROI=1024
+// (array indices of the calibration wrong otherwise)
+//******************************************************************************
+
+	gROOT->SetStyle("Plain");
+	
+//-------------------------------------------
+//Define the data variables
+//-------------------------------------------
+	vector<int16_t> data;
+	vector<int16_t> data_offset;
+	unsigned int data_num;
+	size_t data_n;
+	UInt_t data_px;
+	UInt_t data_roi;
+	
+//-------------------------------------------
+//Get the DRS calibration
+//-------------------------------------------
+	size_t drs_n;
+	vector<float> drs_basemean;
+	vector<float> drs_gainmean;
+	vector<float> drs_triggeroffsetmean;
+	FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
+	
+//-------------------------------------------
+//Check the sizes of the data columns
+//-------------------------------------------
+//	if(drs_n!=data_n)
+//	{
+//		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
+//		return 1;
+//	}
+	
+//-------------------------------------------
+//Create the canvas, plot and title
+//-------------------------------------------
+	char title[500];
+	std::sprintf(title,"Data: %s, DRS: %s, Px %i Ev %i",name,drsname,pixelnr,eventnr);
+	TCanvas *canv = new TCanvas( "canv", "Waveform", 200, 10, 700, 500 );
+//	TProfile *pulseshape = new TProfile("pulseshape", title, 70, -10.5, 59.5);
+	TH2F *pulseshape = new TH2F("pulseshape",title,70,-10.5,59.5,350,-10.5,59.5);
+	pulseshape->GetXaxis()->SetTitle("Time after trigger (bins @2 GHz)");
+	pulseshape->GetYaxis()->SetTitle("Amplitude (mV)");
+	
+	char title_base[500];
+	std::sprintf(title_base,"All samples of Px %i (data: %s, DRS: %s)",pixelnr,name,drsname);
+	
+	char title_spect[500];
+	std::sprintf(title_spect,"Spectrum of all pixels (data: %s, DRS: %s)",name,drsname);
+	TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
+	TH1F *spectrum = new TH1F("spectrum", title_spect,300,0,600);
+	spectrum->GetXaxis()->SetTitle("Pulse size (a.u.)");
+	spectrum->GetYaxis()->SetTitle("Counts");
+	
+//-------------------------------------------
+//Define the fit
+//-------------------------------------------
+	int fit_start = 30;
+	int fit_end = 380;
+	TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp((-x-[11])/[12])",fit_start,fit_end);
+	Double_t par[5+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,10000,30,30};
+//	TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp(-(x-[12])*(x-[12])/(2*[11]*[11]))",fit_start,fit_end);
+//	Double_t par[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
+	f_peaks->SetParameters(par);
+	f_peaks->SetLineColor(2);
+	
+	TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end);
+	f_pedestal->SetLineColor(2);
+	f_pedestal->SetLineWidth(1);
+	
+	TF1* f_peak[n_peaks];
+	for(int i=0; i<n_peaks; i++) {
+		f_peak[i] = new TF1("peak","gaus",fit_start,fit_end);
+		f_peak[i]->SetLineColor(2);
+		f_peak[i]->SetLineWidth(1);
+	}
+	
+//-------------------------------------------
+//Open the file
+//-------------------------------------------
+	fits datafile(name);
+	if (!datafile)
+	{
+	cout << "Couldn't properly open the datafile." << endl;
+	return 1;
+	}
+	FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
+	
+//-------------------------------------------
+//Find the baseline
+//-------------------------------------------
+	float pedestal_mean, pedestal_rms;
+//	FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms);
+	//Return values
+	//Value of maximal probability: -2.225 +- 4.3144
+	pedestal_mean = -2.225;
+	pedestal_rms = 4.3144;
+	
+//-------------------------------------------
+//Oscilloscope
+//-------------------------------------------
+//	float threshold = pedestal_mean+pedestal_rms;
+	float threshold = 2.5;
+
+//***********************
+	FOscilloscope(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
+	canv->cd();
+	pulseshape->Draw();
+	canv->Modified();
+	canv->Update();
+	
+	canv_spect->cd();
+	spectrum->Draw();
+	canv_spect->Modified();
+	canv_spect->Update();
+//***********************
+	fits datafile2("20110804_025.fits");
+	if (!datafile2)
+	{
+	cout << "Couldn't properly open the datafile2." << endl;
+	return 1;
+	}
+	FOpenDataFile(datafile2, data, data_offset, data_num, data_n, data_roi, data_px);
+	FOscilloscope(datafile2, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
+	
+	canv->cd();
+	pulseshape->Draw();
+	canv->Modified();
+	canv->Update();
+	
+	canv_spect->cd();
+	spectrum->Draw();
+	canv_spect->Modified();
+	canv_spect->Update();
+//***********************
+	fits datafile3("20110804_026.fits");
+	if (!datafile3)
+	{
+	cout << "Couldn't properly open the datafile3." << endl;
+	return 1;
+	}
+	FOpenDataFile(datafile3, data, data_offset, data_num, data_n, data_roi, data_px);
+	FOscilloscope(datafile3, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
+	
+	canv->cd();
+	pulseshape->Draw();
+	canv->Modified();
+	canv->Update();
+	
+	canv_spect->cd();
+	spectrum->Draw();
+	canv_spect->Modified();
+	canv_spect->Update();
+//***********************
+	fits datafile4("20110804_027.fits");
+	if (!datafile4)
+	{
+	cout << "Couldn't properly open the datafile4." << endl;
+	return 1;
+	}
+	FOpenDataFile(datafile4, data, data_offset, data_num, data_n, data_roi, data_px);
+	FOscilloscope(datafile4, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
+	
+	canv->cd();
+	pulseshape->Draw();
+	canv->Modified();
+	canv->Update();
+	
+	canv_spect->cd();
+	spectrum->Draw();
+	canv_spect->Modified();
+	canv_spect->Update();
+//***********************
+	spectrum->Fit(f_peaks,"R+");
+	f_peaks->GetParameters(&par[0]);
+	
+	f_peak[0]->SetParameters(par[2],par[0],par[3]);
+	f_peak[1]->SetParameters(par[4],par[0]+par[1],par[5]);
+	f_peak[2]->SetParameters(par[6],par[0]+2*par[1],par[7]);
+	f_peak[3]->SetParameters(par[8],par[0]+3*par[1],par[9]);
+//	f_peak[4]->SetParameters(par[10],par[12],par[11]);
+	f_pedestal->SetParameters(par[10],par[11],par[12]);
+	
+	std::cout << "Crosstalk propability approx. " << ((par[4]*par[5])/(par[2]*par[3]+par[4]*par[5])) << std::endl;
+//-------------------------------------------
+//Draw the data
+//-------------------------------------------
+	canv->cd();
+	pulseshape->Draw();
+	canv->Modified();
+	canv->Update();
+	
+	canv_spect->cd();
+	canv_spect->SetLogy();
+	spectrum->Draw();
+	f_peaks->Draw("same");
+	
+	for(int i=0; i<n_peaks; i++) {
+		f_peak[i]->Draw("same");
+	}
+	
+	f_pedestal->Draw("same");
+	
+	canv_spect->Modified();
+	canv_spect->Update();
+	
+	return 0;
+}
Index: fact/tools/rootmacros/poissonanalysis.C
===================================================================
--- fact/tools/rootmacros/poissonanalysis.C	(revision 12166)
+++ fact/tools/rootmacros/poissonanalysis.C	(revision 12166)
@@ -0,0 +1,237 @@
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TProfile.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TGraph.h>
+
+#include <stdint.h>
+#include <cstdio>
+
+#include "fits.h"
+#include "FOpenDataFile.c"
+#include "FOpenCalibFile.c"
+#include "FPedestal.c"
+#include "FIntFixedPosAllPx.c"
+
+#define n_peaks 2
+
+int poissonanalysis(const char *name, const char *drsname, UInt_t integration_delay, UInt_t integration_size, UInt_t spect_max )
+{
+//******************************************************************************
+//Read a datafile and plot the DRS-calibrated data
+//ATTENTION: only works for ROI=1024
+// (array indices of the calibration wrong otherwise)
+//******************************************************************************
+
+	gROOT->SetStyle("Plain");
+//	TFile f("spectra.root","RECREATE");
+//-------------------------------------------
+//Define the data variables
+//-------------------------------------------
+	vector<int16_t> data;
+	vector<int16_t> data_offset;
+	unsigned int data_num;
+	size_t data_n;
+	UInt_t data_px;
+	UInt_t data_roi;
+	
+//-------------------------------------------
+//Get the DRS calibration
+//-------------------------------------------
+	size_t drs_n;
+	vector<float> drs_basemean;
+	vector<float> drs_gainmean;
+	vector<float> drs_triggeroffsetmean;
+	FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
+	
+//-------------------------------------------
+//Check the sizes of the data columns
+//-------------------------------------------
+//	if(drs_n!=data_n)
+//	{
+//		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
+//		return 1;
+//	}
+	
+//-------------------------------------------
+//Open the file
+//-------------------------------------------
+	fits datafile(name);
+	if (!datafile)
+	{
+	cout << "Couldn't properly open the datafile." << endl;
+	return 1;
+	}
+	FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
+	
+//-------------------------------------------
+//Create the canvas, plot and title
+//-------------------------------------------
+	char title_spect[500];
+	char name_spect[50];
+//	std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",pixelnr, name,drsname);
+	TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
+	TH1F* spectrum[data_px];
+	for(int i=0; i<data_px; i++) {
+		std::sprintf(title_spect,"Spectrum of Px %i, %i-%i (data: %s, DRS: %s)",i,integration_delay,integration_delay+integration_size,name,drsname);
+		std::sprintf(name_spect,"spect%i",i);
+//		spectrum[i] = new TH1F(name_spect,title_spect,400,-150,650);
+		spectrum[i] = new TH1F(name_spect,title_spect,400,-150,spect_max);
+		spectrum[i]->GetXaxis()->SetTitle("Pulse size (a.u.)");
+		spectrum[i]->GetYaxis()->SetTitle("Counts");
+	}
+	
+//	TCanvas *canv_gain= new TCanvas( "canv_gain", "Gain in arbitrary units", 50, 500, 1800, 800 );
+//-------------------------------------------
+//Define the fit
+//-------------------------------------------
+//	int fit_start = 20;
+//	int fit_end = 200;
+//	TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp((-x-[7])/[8])",fit_start,fit_end);
+//	Double_t par_init[5+2*n_peaks] = {70,90,100,15,10,15,10,30,30};
+//	Double_t parread_tmp[5+2*n_peaks];
+//	Double_t parread[data_px][5+2*n_peaks];
+////	TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp(-(x-[12])*(x-[12])/(2*[11]*[11]))",fit_start,fit_end);
+////	Double_t par_init[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
+//	f_peaks->SetParameters(par_init);
+//	f_peaks->SetLineColor(2);
+//	
+//	TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end);
+//	f_pedestal->SetLineColor(2);
+//	f_pedestal->SetLineWidth(1);
+//	
+//	TF1* f_peak[n_peaks];
+//	for(int i=0; i<n_peaks; i++) {
+//		f_peak[i] = new TF1("peak","gaus",fit_start,fit_end);
+//		f_peak[i]->SetLineColor(2);
+//		f_peak[i]->SetLineWidth(1);
+//	}
+	
+//-------------------------------------------
+//Find the baseline
+//-------------------------------------------
+//	float pedestal_mean, pedestal_rms;
+////	FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms);
+//	//Return values
+//	//Value of maximal probability: -2.225 +- 4.3144
+//	pedestal_mean = -2.225;
+//	pedestal_rms = 4.3144;
+	
+//-------------------------------------------
+//Oscilloscope
+//-------------------------------------------
+//	float threshold = pedestal_mean+pedestal_rms;
+//	float threshold = 2.5;
+
+//***********************
+//	UInt_t integration_size = 5;
+//	UInt_t integration_delay = 238;
+
+	FIntFixedPosAllPx(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, integration_size, integration_delay, spectrum);
+	
+	canv_spect->cd();
+	spectrum[0]->Draw();
+	canv_spect->Modified();
+	canv_spect->Update();
+	
+//	Double_t parread_dist[data_px];
+//	Double_t parread_x[data_px];
+//	for(int i=0; i<data_px; i++) {
+//		f_peaks->SetParameters(par_init);
+//		if(!(spectrum[i]->Fit(f_peaks,"NR"))) {
+////			std::cout << "Fit successful." << std::endl;
+//			f_peaks->GetParameters(&parread_tmp[0]);
+//			for(int j=0; j<5+2*n_peaks; j++) {
+//				parread[i][j] = parread_tmp[j];
+//			}
+//			parread_dist[i] = parread_tmp[1];
+//			parread_x[i] = i;
+//		}
+//		else {
+//			for(int j=0; j<5+2*n_peaks; j++) {
+//				parread[i][j] = 0;
+//			}
+//			parread_dist[i] = 0;
+//			parread_x[i] = i;
+//		}
+//	}
+//	
+//	canv_gain->cd();
+//	TGraph* plot_gain = new TGraph(data_px,parread_x,parread_dist);
+//	plot_gain->GetYaxis()->SetRangeUser(0,150);
+//	plot_gain->Draw("A*");
+//	canv_gain->Modified();
+//	canv_gain->Update();
+//	
+//	canv_spect->cd();
+//	canv_spect->SetLogy();
+	
+//-------------------------------------------
+//Draw the data where the gain is out of limits
+//-------------------------------------------
+//	char temp = 'x';
+//	std::cout << "Plot the spectra out of limits: Enter for next, 'a' to abort." << std::endl;
+//	for(int i=0; i<data_px; i++) {
+//		if((parread_dist[i]<80)||(parread_dist[i]>105)) {
+//			std::cout << "Gain out of limits: Px " << i << " at " << parread_dist[i] << std::endl;
+//			spectrum[i]->Draw();
+//			
+//			f_peaks->SetParameters(&parread[i][0]);
+//			f_peak[0]->SetParameters(parread[i][2],parread[i][0],parread[i][3]);
+//			f_peak[1]->SetParameters(parread[i][4],parread[i][0]+parread[i][1],parread[i][5]);
+////			f_peak[2]->SetParameters(parread[i][6],parread[i][0]+2*parread[i][1],parread[i][7]);
+////			f_peak[3]->SetParameters(parread[i][8],parread[i][0]+3*parread[i][1],parread[i][9]);
+////			f_peak[4]->SetParameters(parread[i][10],parread[i][12],parread[i][11]);
+//			f_pedestal->SetParameters(parread[i][6],parread[i][7],parread[i][8]);
+//			
+//			f_peaks->Draw("same");
+//			for(int i=0; i<n_peaks; i++) {
+//				f_peak[i]->Draw("same");
+//			}
+//			f_pedestal->Draw("same");
+//			
+//			canv_spect->Modified();
+//			canv_spect->Update();
+//			temp='x';
+//			while ((temp!='\n') && (temp!='a'))
+//			    cin.get(temp);
+//			if(temp=='a') break;
+//		}
+//	}
+//-------------------------------------------
+//Draw the data
+//-------------------------------------------
+	canv_spect->cd();
+	Int_t plotspect_nr;
+	std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
+	std::cin >> plotspect_nr;
+	while(plotspect_nr>=0) {
+		plotspect_nr = plotspect_nr % data_px;
+		spectrum[plotspect_nr]->Draw();
+		
+//		f_peaks->SetParameters(&parread[plotspect_nr][0]);
+//		f_peak[0]->SetParameters(parread[plotspect_nr][2],parread[plotspect_nr][0],parread[plotspect_nr][3]);
+//		f_peak[1]->SetParameters(parread[plotspect_nr][4],parread[plotspect_nr][0]+parread[plotspect_nr][1],parread[plotspect_nr][5]);
+////		f_peak[2]->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][0]+2*parread[plotspect_nr][1],parread[plotspect_nr][7]);
+////		f_peak[3]->SetParameters(parread[plotspect_nr][8],parread[plotspect_nr][0]+3*parread[plotspect_nr][1],parread[plotspect_nr][9]);
+////		f_peak[4]->SetParameters(parread[plotspect_nr][10],parread[plotspect_nr][12],parread[plotspect_nr][11]);
+//		f_pedestal->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][7],parread[plotspect_nr][8]);
+//		
+//		f_peaks->Draw("same");
+//		for(int i=0; i<n_peaks; i++) {
+//			f_peak[i]->Draw("same");
+//		}
+//		f_pedestal->Draw("same");
+		
+		canv_spect->Modified();
+		canv_spect->Update();
+		
+		std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
+		std::cin >> plotspect_nr;
+	}
+	
+//	f->Write();
+	return 0;
+}
Index: fact/tools/rootmacros/test.C
===================================================================
--- fact/tools/rootmacros/test.C	(revision 12166)
+++ fact/tools/rootmacros/test.C	(revision 12166)
@@ -0,0 +1,47 @@
+#include "fits.h"
+
+#include <stdint.h>
+
+int test(const char *name)
+{
+    fits file(name);
+    if (!file)
+    {
+        cout << "Couldn't properly open the file." << endl;
+        return 1;
+    }
+
+    cout << "-------------------- Header -------------------" << endl;
+    file.PrintKeys();
+    cout << "------------------- Columns -------------------" << endl;
+    file.PrintColumns();
+
+    cout << "--------------------- Data --------------------" << endl;
+
+    unsigned int num;
+    short type;
+
+    const size_t n   = file.GetN("Data");     // Size of column "Data"
+    const UInt_t roi = file.GetUInt("NROI");  // Value from header
+
+    vector<int16_t> data(n);
+
+    file.SetRefAddress("EventNum",    num);
+    file.SetRefAddress("TriggerType", type);
+    file.SetVecAddress("Data",        data);
+
+    for (size_t i=0; i<file.GetNumRows(); i++)
+    {
+        file.GetRow(i);
+        cout << "Event number: " << num << endl;
+
+        for (int j=0; j<1440; j++)
+            for (UInt_t k=0; k<roi; k++)
+            {
+                int16_t sample = data[j*roi+k];
+                sample+=0; // Supress compiler warning
+            }
+    }
+
+    return 0;
+}
Index: fact/tools/rootmacros/testFPedestalAllPx.C
===================================================================
--- fact/tools/rootmacros/testFPedestalAllPx.C	(revision 12166)
+++ fact/tools/rootmacros/testFPedestalAllPx.C	(revision 12166)
@@ -0,0 +1,106 @@
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TH1.h>
+
+#include <stdint.h>
+#include <cstdio>
+
+#include "fits.h"
+#include "FOpenDataFile.c"
+#include "FOpenCalibFile.c"
+#include "FPedestalAllPx.c"
+
+int testFPedestalAllPx(const char *name, const char *drsname)
+{
+//******************************************************************************
+//Read a datafile and plot the DRS-calibrated data
+//ATTENTION: only works for ROI=1024
+// (array indices of the calibration wrong otherwise)
+//Example call in ROOT:
+//root [74] .x calscope.C++("20110804_024.fits","20110804_023.drs.fits",10,1348)
+//T. Krähenbühl, August 2011, tpk@phys.ethz.ch
+//******************************************************************************
+
+	gROOT->SetStyle("Plain");
+	
+//-------------------------------------------
+//Open the file
+//-------------------------------------------
+	fits datafile(name);
+	if (!datafile)
+	{
+	cout << "Couldn't properly open the datafile." << endl;
+	return 1;
+	}
+	
+//-------------------------------------------
+//Get the data
+//-------------------------------------------
+	vector<int16_t> data;
+	vector<int16_t> data_offset;
+	unsigned int data_num;
+	size_t data_n;
+	UInt_t data_px;
+	UInt_t data_roi;
+	FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
+	
+//-------------------------------------------
+//Get the DRS calibration
+//-------------------------------------------
+	size_t drs_n;
+	vector<float> drs_basemean;
+	vector<float> drs_gainmean;
+	vector<float> drs_triggeroffsetmean;
+	FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
+	
+//-------------------------------------------
+//Check the sizes of the data columns
+//-------------------------------------------
+	if(drs_n!=data_n)
+	{
+		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
+		return 1;
+	}
+	
+	
+//-------------------------------------------
+//Create the title
+//-------------------------------------------
+//	char title[500];
+//	std::sprintf(title,"Data: %s, DRS: %s, Px %i Ev %i",name,drsname,pixelnr,eventnr);
+	
+//-------------------------------------------
+//Get the pedestals & RMS
+//-------------------------------------------
+	float pedestal_mean[data_px];
+	float pedestal_rms[data_px];
+	FPedestalAllPx(datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_px, data_roi, pedestal_mean, pedestal_rms);
+//	vector<float> calevent(data_px*data_roi); //Vector for the calibrated event
+//	FGetCalEvent(data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, calevent, data_px, data_roi);
+	
+//-------------------------------------------
+//Draw the data
+//-------------------------------------------
+	TCanvas *canv_mean= new TCanvas( "canv_mean", "mean in mV", 30, 300, 900, 500 );
+	TH1F* histo_mean = new TH1F("histo_mean","Value of maximal probability",2000,-99.5,100.5);
+	TCanvas *canv_rms= new TCanvas( "canv_rms", "RMS in mV", 960, 300, 900, 500 );
+	TH1F* histo_rms = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
+	
+	for (int j=0; j<data_px; j++)
+	{
+		histo_mean->Fill(pedestal_mean[j]);
+		histo_rms->Fill(pedestal_rms[j]);
+	}
+//	histo_mean->GetYaxis()->SetRangeUser(0,150);
+	
+	canv_mean->cd();
+	histo_mean->Draw();
+	canv_mean->Modified();
+	canv_mean->Update();
+	canv_rms->cd();
+	histo_rms->Draw();
+	canv_rms->Modified();
+	canv_rms->Update();
+	
+	return 0;
+}
Index: fact/tools/rootmacros/testfir.C
===================================================================
--- fact/tools/rootmacros/testfir.C	(revision 12166)
+++ fact/tools/rootmacros/testfir.C	(revision 12166)
@@ -0,0 +1,62 @@
+#include "TRandom2.h"
+#include "TH1F.h"
+#include "TCanvas.h"
+
+#include <stdio.h>
+#include <iostream>
+
+
+#define FAD_MAX_SAMPLES 1024
+
+vector<float> Vmeas(FAD_MAX_SAMPLES);
+vector <float> Vfir(FAD_MAX_SAMPLES);
+vector <float> Vpulse(FAD_MAX_SAMPLES);
+#include "factfir.C"
+
+
+#define k 16
+double a[k];
+double b=k;
+
+void testfir(){
+
+	for (int j=0; j<k; j++){
+		a[j]=1;
+	}
+
+	TRandom2 myrandom;
+	myrandom.RndmArray(1024,Vmeas);
+	
+	int pulse_start=350;
+	float pulse_height=0.6;
+	int pulse_len = 60;
+	for (int i = pulse_start; i<pulse_start + pulse_len; i++){
+		Vpulse[i] += pulse_height/pulse_len * (pulse_len-(i-pulse_start));
+		// in cas e you want the pulse in the data
+		Vmeas[i] += Vpulse[i];	 
+	}
+
+
+	factfir(b,a,k, Vfir);
+
+	TH1F *meas = new TH1F("Vmeas","Vmeas test",FAD_MAX_SAMPLES, -0.5 , FAD_MAX_SAMPLES-0.5);
+	TH1F *fir = new TH1F("Vfir","Vfir test",FAD_MAX_SAMPLES, -0.5 , FAD_MAX_SAMPLES-0.5);
+	TH1F *pulse = new TH1F("Vpulse","Vpulse test",FAD_MAX_SAMPLES, -0.5 , FAD_MAX_SAMPLES-0.5);
+
+	for (int i=0; i<FAD_MAX_SAMPLES; i++){
+		meas->SetBinContent(i,Vmeas[i]);
+		fir->SetBinContent(i,Vfir[i]);
+		pulse->SetBinContent(i,Vpulse[i]);
+	}
+
+	TCanvas *c = new TCanvas();
+	c->Divide(1,3);
+	c->cd(1);
+	meas->Draw("HIST");
+	c->cd(2);
+	fir->Draw("HIST");
+	c->cd(3);
+	pulse->Draw("HIST");
+
+
+}
Index: fact/tools/rootmacros/zerosearch.C
===================================================================
--- fact/tools/rootmacros/zerosearch.C	(revision 12166)
+++ fact/tools/rootmacros/zerosearch.C	(revision 12166)
@@ -0,0 +1,44 @@
+#include <vector>
+
+// searches for zero crossings in a given vector of floats
+// for zero crossings in falling edge 
+// give edge = -1
+// 
+// 
+// returns pointer to vetctor of ints
+
+vector<int> *zerosearch(vector<float> &input, int edge =-1 , int pre=10, int post=30 ){
+
+	
+	vector<int> * zeroPositions = new vector<int>;
+	
+	for (int sl =pre ; sl < input.size()-post; sl++){
+//cout << "sl:" << sl << endl;
+		if (input[sl] * input[sl-1] < 0 || input[sl]==0.0 ){ // sign change --> zero crossing OR really zero
+
+			if ( (input[sl-1] - input[sl]) * edge < 0){ // this is the crossing the user wanted 
+				
+				// check if we go lower than a certain limit in the next few slices
+				for ( int lala=0; lala<post; lala++){
+					if (input[sl+lala] < -1.5) {
+						zeroPositions->push_back(sl);
+						sl += lala;
+						break;
+					}
+				}
+
+			} else if ( (input[sl-1] - input[sl]) * edge > 0){ // this is the zero x-ing the user did not want
+
+				// do nothing
+
+			} else { // sl and sl-1 are equal .. don't know waht do to...
+	
+			}		
+
+		} 
+
+	} // end of loop over slices - between pre and post
+
+	
+	return zeroPositions;
+}
