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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#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 data; +vector 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 drs_basemean; +vector drs_gainmean; +vector drs_triggeroffsetmean; + +int FOpenDataFile( fits & ); + + +vector Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude +vector N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors +vector N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors +vector Vcorr(FAD_MAX_SAMPLES); // corrected Values +vector Vdiff(FAD_MAX_SAMPLES); // numerical derivative + +vector Vslide(FAD_MAX_SAMPLES); // sliding average result +vector Vcfd(FAD_MAX_SAMPLES); // CDF result +vector 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 a_slide(k_slide, 1); + double b_slide = k_slide; + +// CFD filter settings + int k_cfd = 10; + vector 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 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() ; +} +