Index: /fact/tools/rootmacros/zippedfitstest.C =================================================================== --- /fact/tools/rootmacros/zippedfitstest.C (revision 12193) +++ /fact/tools/rootmacros/zippedfitstest.C (revision 12193) @@ -0,0 +1,437 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +#define HAVE_ZLIB +#include "izstream.h" +#include "fits.h" +//#include "TPKplotevent.c" +//#include "FOpenDataFile.c" +#include "FOpenCalibFile.c" + +#include "discriminator.h" +#include "discriminator.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.); + +void BookHistos( int ); + +// Create a canvas +TCanvas* CW; +TCanvas* cFilter; + +int spikeDebug = 0; + + +int fana( + char *datafilename = "../raw/20111011_026.fits.gz", + const char *drsfilename = "../raw/20111011_022.drs.fits.gz", + int pixelnr = 4, + 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; + + +//------------------------------------------- +//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; + 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 ]); + } + + // 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); + + vector *my = discriminator( Vslide, 3.5, 100 ); + for (int p=0; psize(); p++ ){ + sp->Fill(my->at(p).maxVal); + } + delete my; + + + if ( spikeDebug ){ + for ( unsigned int sl = 0; sl < data_roi; sl++){ + h[tVslide].SetBinContent( sl, Vslide[sl] ); + h[tVcfd].SetBinContent( sl, Vcfd[sl] ); + h[tVcfd2].SetBinContent( sl, Vcfd2[sl] ); + } + } +/* + vector * zeros = zerosearch( Vcfd2, -1, 10, 20 ); + if (zeros->size() == 0 ){ + continue; + } + // check value of Vside at zero position + for ( int i=0; isize(); i++){ + cout << zeros->at(i) << ":\t" << Vslide[ zeros->at(i) ]<Fill(Vslide[zeros->at(i)]); + } +*/ + + 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( tVtest + 1); + // gPad->SetGrid(); + // h[tVtest].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; + } + + + } + + 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, 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); + +} +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() ; +} +