Index: fact/tools/rootmacros/fana2.C =================================================================== --- fact/tools/rootmacros/fana2.C (revision 12604) +++ fact/tools/rootmacros/fana2.C (revision 12604) @@ -0,0 +1,296 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#define HAVE_ZLIB +#include "fits.h" + +#include "openFits.h" +#include "openFits.c" + +#include "DrsCalibration.h" +#include "DrsCalibration.C" + +#include "SpikeRemoval.h" +#include "SpikeRemoval.C" + +#define NPIX 1440 +#define NCELL 1024 + +// data access and treatment +#define FAD_MAX_SAMPLES 1024 + +int NEvents; +vector Data; // vector, which will be filled with raw data +vector StartCells; // vector, which will be filled with DRS start positions +unsigned int EventID; // index of the current event +UInt_t RegionOfInterest; // Width of the Region, read out of the DRS +UInt_t NumberOfPixels; // Total number of pixel, read out of the camera +size_t PXLxROI; // Size of column "Data" = #Pixel x ROI + +int NBSLeve = 1000; + +size_t TriggerOffsetROI, RC; +vector Offset, Gain, TriggerOffset; + + +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 SpikeEst(FAD_MAX_SAMPLES); // corrected Values +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" + +void computeSpikeEstimator( int ); + +// histograms +const int Ntypes = 7; +const unsigned int // arranged by Dominik + tAmeas = 0, + tSpikeEst = 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.); + +void BookHistos( int ); + +// Create a canvas +TCanvas* CW; +TCanvas* cFilter; + +int spikeDebug = 1; + + +int fana2( + 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 = NULL; + NEvents = openDataFits( + datafilename, + &datafile, + Data, + StartCells, + EventID, + RegionOfInterest, + NumberOfPixels, + PXLxROI ); +- + printf( "number of events in file: %d\n", NEvents ); + if (NEvents == 0){ + cout << "return code of openDataFits:" << datafilename<< endl; + cout << "is zero -> aborting." << endl; + return 1; + } + + // 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 +//------------------------------------------- + + RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI, 3); + if (RC == 0){ + cout << "return code of openCalibFits:" << drsfilename << endl; + cout << "is zero -> aborting." << endl; + return 1; + } + +// Book the histograms + BookHistos( RegionOfInterest ); + +// Loop over events + cout << "--------------------- Data --------------------" << endl; + + // float value; + + // TH1F * sp = new TH1F("spektrum", "test of Stepktrum", 256, -0.5, 63.5); + + size_t calib_RC = 1; + for ( int ev = firstevent; ev < firstevent + nevents; ev++) { + + datafile->GetRow( ev ); + if (ev % 50 ==0){ + cout << "Event number: " << EventID << endl; + } + + // get the data of this pixel from the Data vector + // apply the Drs Calibration and cut off 12 slices at the beginning + // and at the end. + calib_RC = applyDrsCalibration( Ameas, pixelnr,12,12, Offset, Gain, TriggerOffset, + RegionOfInterest, Data, StartCells); + if (calib_RC == 0){ + break; + } + + computeSpikeEstimator( RegionOfInterest ); + removeSpikes( Ameas, Vcorr ); + sliding_avg( Vcorr, Vslide, 8 ); + + // for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){ + // hPixelCellData.Fill(sl, Vcorr[ sl ]); + //} + + // filter Vcorr with sliding average using FIR filter function + factfir(b_slide , a_slide, k_slide, Vcorr, Vslide); + + // filter Vslide with CFD using FIR filter function + factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd); + // filter Vcfd with sliding average using FIR filter function + factfir(b_slide , a_slide, k_slide, Vcfd, Vcfd2); + + if ( spikeDebug ){ + for ( unsigned int sl = 0; sl < RegionOfInterest; sl++ ){ + + h[tAmeas].SetBinContent( sl, Ameas[ sl ] ); + h[tVcorr].SetBinContent( sl, Vcorr[ sl ] ); + + h[tVslide].SetBinContent( sl, Vslide[ sl ] ); + h[tVcfd].SetBinContent( sl, Vcfd[ sl ] ); + h[tVcfd2].SetBinContent( sl, Vcfd2[ sl ] ); + } + } + + if ( spikeDebug ){ + CW->cd( tAmeas + 1); + gPad->SetGrid(); + h[tAmeas].Draw(); + + CW->cd( tSpikeEst + 1); + gPad->SetGrid(); + h[tSpikeEst].Draw(); + + CW->cd( tVcorr + 1); + gPad->SetGrid(); + h[tVcorr].Draw(); + + cFilter->cd( Ntypes - tVslide ); + cFilter->cd(1); + gPad->SetGrid(); + h[tVslide].Draw(); + cFilter->cd( Ntypes - tVcfd ); + cFilter->cd(2); + gPad->SetGrid(); + h[tVcfd].Draw(); + TLine zeroline(0, 0, 1024, 0); + zeroline.SetLineColor(kBlue); + zeroline.Draw(); + + cFilter->cd( Ntypes - tVcfd2 ); + cFilter->cd(3); + gPad->SetGrid(); + h[tVcfd2].Draw(); + + zeroline.Draw(); + + CW->Update(); + cFilter->Update(); + + //Process gui events asynchronously during input + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + timer.TurnOn(); + TString input = Getline("Type 'q' to exit, to go on: "); + timer.TurnOff(); + if (input=="q\n") break; + } + + + } + + return( 0 ); +} + +void computeSpikeEstimator( int Samples ){ + +// compute the mean of the left and right neighbors of a channel + + for( int i = 1; i < Samples-1; i++){ + + N1mean[ i ] = ( Ameas[i-1] + Ameas[i+1] ) / 2.; + SpikeEst[ i ] = Ameas[ i ] - N1mean[ i ]; + h[tSpikeEst].SetBinContent( i, SpikeEst[ i ] ); + } +} // end of computeN1mean computation + + +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, 3); + cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600); + cFilter->Divide(1, 3); + + // hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024); + +}