Index: fact/tools/rootmacros/fbsl.C
===================================================================
--- fact/tools/rootmacros/fbsl.C	(revision 12206)
+++ fact/tools/rootmacros/fbsl.C	(revision 12206)
@@ -0,0 +1,416 @@
+#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 <TFile.h>
+
+#include <stdio.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> Vcfd2(FAD_MAX_SAMPLES);    // CDF result + 2nd sliding average
+
+#include "factfir.C"
+
+float getValue( int, int );
+void  computeN1mean( int );
+void  removeSpikes( int );
+
+// histograms
+const int Ntypes = 7;
+const unsigned int  // arranged by Dominik
+	tAmeas 	= 0, 
+	tN1mean = 1, 
+	tVcorr  = 2,
+	tVtest  = 3,
+	tVslide = 4,
+	tVcfd   = 5,
+	tVcfd2  = 6;
+
+TH1F* h;
+TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
+                    // x = pixel id, y = DRS cell id
+TH2F hPixelCellData(
+    "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
+TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
+TH1F *hMeanBsl, *hpltMeanBsl;
+TH1F *hRmsBsl, *hpltRmsBsl; 
+TObjArray hList;
+TObjArray hListBaseline;
+
+void BookHistos( int );
+void SaveHistograms( char * );
+
+// Create a canvas
+TCanvas* CW;
+TCanvas* cFilter;
+
+int spikeDebug = 0;
+int searchSinglesPeaks = 0;
+
+
+int fbsl( 
+	char *datafilename 		= "../raw/20110916_025.fits",
+	const char *drsfilename = "../raw/20110916_024.drs.fits",
+	int pixelnr 			= 0,
+	int firstevent 			= 0, 
+	int nevents 			= -1 ){
+// read and analyze FACT raw data
+
+// 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;
+
+
+    //  FILE *pedFile;
+	// pedFile = fopen("t.txt","u");
+    // fprintf(pedFile,"# Pedestal Data");
+    // fclose( pedFile );
+//-------------------------------------------
+//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
+    printf("call BookHistos\n");
+    BookHistos( data_roi );
+    printf("returned from BookHistos\n");
+    // Loop over events
+	cout << "--------------------- Data --------------------" << endl;
+
+    float value;
+    
+    // loop over events
+//    for ( int ev = firstevent; ev < nevents; ev++) {
+    for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
+
+	    datafile.GetRow( ev );
+		
+		if ( ev % 100 == 0 ){
+			cout << "Event number: " << data_num << endl;
+		}
+		
+        // loop over pixel
+        for ( int pix = 0; pix < 1440; pix++ ){
+
+            // hStartCell->Fill( pix, data_offset[pix] );
+        
+            // loop over DRS slices
+            for ( unsigned int sl = 0; sl < data_roi; sl++){
+
+                Ameas[ sl ] = getValue(sl, pix);
+                
+            }
+            // printf("Ameas[100] = %f\n", Ameas[100] );
+
+            computeN1mean( data_roi ); // prepare spike removal
+            removeSpikes( data_roi );  // output in Vcorr
+
+            // filter Vcorr with sliding average using FIR filter function
+            factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
+
+            for ( unsigned int sl = 0; sl < data_roi; sl++){
+               // hPixelCellData.Fill( sl, Vcorr[sl] );
+               hBaseline[pix]->Fill( Vslide[sl] ) ;
+            }   
+	    }
+    }
+
+    for (int pix = 0; pix < NPIX; pix++) {
+        //printf("Maximum bin pix[%d] %f\n", pix ,hBaseline[pix]->GetMaximumBin() );
+
+        float vmaxprob = hBaseline[pix]->GetXaxis()->GetBinCenter(
+            hBaseline[pix]->GetMaximumBin() );
+        hMeanBsl->Fill( vmaxprob );
+        hpltMeanBsl->SetBinContent(pix+1, vmaxprob );
+
+        hRmsBsl->Fill(hBaseline[pix]->GetRMS() );
+        hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() );   
+    }
+
+    SaveHistograms( datafilename );
+
+    TCanvas * cMeanBsl = new TCanvas();
+    cMeanBsl->cd();
+    hMeanBsl->Draw();
+    //canv_mean->Modified();
+    cMeanBsl->Update();
+
+    TCanvas * cRmsBsl = new TCanvas();
+    cRmsBsl->cd();
+    hRmsBsl->Draw();
+    //canv_rms->Modified();
+    cMeanBsl->Update();
+
+	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 = 2; i < Samples-2 ; 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];
+
+            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.;
+                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
+}
+
+void computeN1mean( int Samples ){
+// compute the mean of the left and right neighbors of a channel
+
+    for( int i = 2; i < Samples - 2; i++){
+/*        if (i == 0){ // use right sample es mean
+            N1mean[i] = Ameas[i+1];
+        }
+        else if ( i == Samples-1 ){ //use left sample as mean
+            N1mean[i] = Ameas[i-1];
+        }
+        else{
+            N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
+        }
+*/
+        N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
+    }
+} // end of computeN1mean computation
+
+float getValue( int slice, int pixel ){
+
+    const float dconv = 2000/4096.0;
+
+    float vraw, vcal;
+
+    unsigned int pixel_pt;
+    unsigned int slice_pt;
+    unsigned int cal_pt;
+    unsigned int drs_cal_offset;
+
+    // printf("pixel = %d, slice = %d\n", slice, pixel);
+
+    pixel_pt = pixel * 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
+
+    // histograms for baseline extraction
+    char hName[500];
+    char hTitle[500];
+
+    TH1F *h;
+    
+    printf("inside BookHistos\n");
+
+    for( int i = 0; i < NPIX; i++ ) {
+
+        // printf("call sprintf [%d]\n", i );
+        sprintf(&hTitle[0],"all events all slices of pixel %d", i);
+        sprintf(&hName[0],"base%d", i);
+        // printf("call sprintf [%d] done\n", i );
+
+        h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );
+
+        // printf("histo booked\n");
+        h->GetXaxis()->SetTitle( "Sample value (mV)" );
+        h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
+        // printf("histo title set\n");
+        hListBaseline.Add( h );
+        // printf("histo added to List\n");
+        hBaseline[i] = h;
+        // printf("histo assigned to array\n");
+    }
+
+    printf("made HBaseline * 1440\n");
+
+    hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
+    hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
+    hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
+    hList.Add( hMeanBsl );
+
+    hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
+    hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
+    hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
+    hList.Add( hpltMeanBsl );
+
+    hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
+    hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
+    hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
+    hList.Add( hRmsBsl );
+
+    hpltRmsBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
+    hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
+    hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
+    hList.Add( hpltRmsBsl );
+}
+
+
+void SaveHistograms( char * loc_fname ){
+
+  TString fName; // name of the histogram file
+
+  /* create the filename for the histogram file */
+  fName = loc_fname; // use the name of the tree file
+  fName.Remove(fName.Length() - 5); // remove the extension .root
+  fName += "_histo.root"; // add the new extension
+
+  TFile tf( fName, "RECREATE"); // create the histogram file (replace if already existing)
+
+  hList.Write(); // write the major histograms into the top level directory 
+  tf.mkdir("BaselineHisto"); tf.cd("BaselineHisto"); // go to new subdirectory
+  hListBaseline.Write(); // write histos into subdirectory
+  
+  tf.Close(); // close the file
+
+} // end of function: void ana::SaveHistograms( void )
+
+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() ;
+}
+
