Index: /fact/tools/rootmacros/fbsl.C =================================================================== --- /fact/tools/rootmacros/fbsl.C (revision 12513) +++ /fact/tools/rootmacros/fbsl.C (revision 12514) @@ -16,7 +16,15 @@ #define HAVE_ZLIB #include "fits.h" -#include "FOpenCalibFile.c" - -#include "zerosearch.C" + +#include "openFits.h" +#include "openFits.c" + +#include "DrsCalibration.h" +#include "DrsCalibration.C" + +#include "SpikeRemoval.h" +#include "SpikeRemoval.C" + +//#include "zerosearch.C" #include "factfir.C" @@ -27,37 +35,16 @@ #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; - 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 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 drs_n; -vector drs_basemean; -vector drs_gainmean; -vector drs_triggeroffsetmean; - -size_t FOpenDataFile( - const char *datafilename, // path to fits file containing FACT raw data - fits * * datafile, // pointer to pointer, where to return the fits object - 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 - // this can be used, to x-check RegionOfInterest and NumberOfPixels - int VerbosityLevel=1 -); +size_t TriggerOffsetROI, RC; +vector Offset, Gain, TriggerOffset; @@ -73,7 +60,4 @@ -float getValue( int, int ); -void computeN1mean( int ); -void removeSpikes( int ); // histograms @@ -107,5 +91,4 @@ int searchSinglesPeaks = 0; - int fbsl( const char *datafilename = "path-to-datafile.fits.gz", @@ -117,26 +100,38 @@ int firstpixel = 0, int npixel = -1, - bool produceGraphic = false + bool produceGraphic = true ){ fits * datafile = NULL; - NEvents = FOpenDataFile( - datafilename, - &datafile, - Data, - StartCells, - EventID, - RegionOfInterest, - NumberOfPixels, - PXLxROI); + + 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; - - if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels; - - FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n); + if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels; + + RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI); + if (RC == 0){ + cout << "return code of openCalibFits:" << drsfilename << endl; + cout << "is zero -> aborting." << endl; + return 1; + } + BookHistos( RegionOfInterest, npixel ); @@ -151,20 +146,22 @@ } - // loop over pixel - for ( int pix = firstpixel ; pix < npixel+firstpixel ; pix++ ){ - - // loop over DRS slices - for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){ - Ameas[ sl ] = getValue(sl, pix); - } - - computeN1mean( RegionOfInterest ); // prepare spike removal - removeSpikes (RegionOfInterest ); // output in Vcorr - - // filter Vcorr with sliding average using FIR filter function - // 8 is here the HalfWidth of the sliding average window. - sliding_avg(Vcorr, Vslide, 8); - - // factfir(b_slide , a_slide, k_slide, Vcorr, Vslide); + // loop over pixel + for ( int pix = firstpixel ; pix < npixel+firstpixel ; pix++ ){ + + + // 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. + applyDrsCalibration( Ameas,pix,12,12, + Offset, Gain, TriggerOffset, + RegionOfInterest, Data, StartCells); + + // finds spikes in the raw data, and interpolates the value + // spikes are: 1 or 2 slice wide, positive non physical artifacts + removeSpikes (Ameas, Vcorr); + + // filter Vcorr with sliding average using FIR filter function + // 8 is here the HalfWidth of the sliding average window. + sliding_avg(Vcorr, Vslide, 8); for ( unsigned int sl = 0; sl 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 * RegionOfInterest; - slice_pt = pixel_pt + slice; - drs_cal_offset = ( slice + StartCells[ pixel ] )%RegionOfInterest; - 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 ); } @@ -375,59 +288,2 @@ } // end of function: void ana::SaveHistograms( void ) -size_t FOpenDataFile( - const char *datafilename, // path to fits file containing FACT raw data - fits * * datafile, // ptr to pointer, where to return the fits object - 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 - // this can be used, to x-check RegionOfInterest and NumberOfPixels - int VerbosityLevel // -) { - size_t NumberOfEvents; - *datafile = new fits(datafilename); - if (!(*(*datafile))) { - if (VerbosityLevel > 0) - cout << "Couldn't properly open the datafile: " << datafilename << endl; - return 0; - } - - NumberOfEvents = (*datafile)->GetNumRows(); - if (NumberOfEvents < 1){ - if (VerbosityLevel > 0){ - cout << "Warning in FOpenDataFile of file: " << datafilename << endl; - cout << "the file contains no events." << endl; - } - } - - RegionOfInterest = (*datafile)->GetUInt("NROI"); - NumberOfPixels = (*datafile)->GetUInt("NPIX"); - PXLxROI = (*datafile)->GetN("Data"); - - if ( RegionOfInterest * NumberOfPixels != PXLxROI) // something in the file went wrong - { - if (VerbosityLevel > 0){ - cout << "Warning in FOpenDataFile of file: " << datafilename << endl; - cout << "RegionOfInterest * NumberOfPixels != PXLxROI" << endl; - cout << "--> " << RegionOfInterest; - cout << " * " << NumberOfPixels; - cout << " = " << RegionOfInterest * NumberOfPixels; - cout << ", but PXLxROI =" << PXLxROI << endl; - } - return 0; - } - - //Set the sizes of the data vectors - Data.resize(PXLxROI, 0); - StartCells.resize(NumberOfPixels, 0); - - //Link the data to variables - (*datafile)->SetRefAddress("EventNum", EventID); - (*datafile)->SetVecAddress("Data", Data); - (*datafile)->SetVecAddress("StartCellData", StartCells); - - return NumberOfEvents; -} -