Changeset 12514

Timestamp:

11/14/11 19:14:49 (13 years ago)

Author:

neise

Message:

try with new calibration

File:

• fact/tools/rootmacros/fbsl.C (modified) (8 diffs)

fact/tools/rootmacros/fbsl.C

@@ -16 +16 @@
 #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 +35 @@
 #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;
-        vector<int16_t> Data;                           // vector, which will be filled with raw data
-        vector<int16_t> 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<int16_t> Data;                           // vector, which will be filled with raw data
+vector<int16_t> 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<float> drs_basemean;
-vector<float> drs_gainmean;
-vector<float> 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<int16_t> &Data,                  // vector, which will be filled with raw data
-        vector<int16_t> &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<float> Offset, Gain, TriggerOffset;
 
 
@@ -73 +60 @@
 
 
-float getValue( int, int );
-void  computeN1mean( int );
-void  removeSpikes( int );
 
 // histograms
@@ -107 +91 @@
 int searchSinglesPeaks = 0;
 
-
 int fbsl( 
         const char *datafilename                = "path-to-datafile.fits.gz",
@@ -117 +100 @@
         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 +146 @@
                 }
                 
-        // 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 <RegionOfInterest ; sl++){
@@ -215 +212 @@
                 }
         return( 0 );
-}
-
-void removeSpikes(int Samples){
-
-    const float fract = 0.8;
-    float x, xp, xpp;
-
-    // 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 * 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 +288 @@
 } // 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<int16_t> &Data,                  // vector, which will be filled with raw data
-        vector<int16_t> &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;
-}
-