Ignore:
Timestamp:
11/14/11 19:14:49 (13 years ago)
Author:
neise
Message:
try with new calibration
File:
1 edited

Legend:

Unmodified
Added
Removed
  • fact/tools/rootmacros/fbsl.C

    r12356 r12514  
    1616#define HAVE_ZLIB
    1717#include "fits.h"
    18 #include "FOpenCalibFile.c"
    19 
    20 #include "zerosearch.C"
     18
     19#include "openFits.h"
     20#include "openFits.c"
     21
     22#include "DrsCalibration.h"
     23#include "DrsCalibration.C"
     24
     25#include "SpikeRemoval.h"
     26#include "SpikeRemoval.C"
     27
     28//#include "zerosearch.C"
    2129#include "factfir.C"
    2230
     
    2735#define FAD_MAX_SAMPLES 1024
    2836
    29 //vector<int16_t> data;
    30 //vector<int16_t> data_offset;
    31 //unsigned int data_num;
    32 //size_t data_n;
    33 //UInt_t data_px;
    34 //UInt_t data_roi;
    35         vector<int16_t> Data;                           // vector, which will be filled with raw data
    36         vector<int16_t> StartCells;     // vector, which will be filled with DRS start positions
    37         unsigned int EventID;                           // index of the current event
    38         UInt_t RegionOfInterest;                // Width of the Region, read out of the DRS
    39         UInt_t NumberOfPixels;                  // Total number of pixel, read out of the camera
    40         size_t PXLxROI;                                                 // Size of column "Data" = #Pixel x ROI
    41 
    4237int NEvents;
     38vector<int16_t> Data;                           // vector, which will be filled with raw data
     39vector<int16_t> StartCells;     // vector, which will be filled with DRS start positions
     40unsigned int EventID;                           // index of the current event
     41UInt_t RegionOfInterest;                // Width of the Region, read out of the DRS
     42UInt_t NumberOfPixels;                  // Total number of pixel, read out of the camera
     43size_t PXLxROI;                                                 // Size of column "Data" = #Pixel x ROI
     44
    4345int NBSLeve = 1000;
    4446
    45 size_t drs_n;
    46 vector<float> drs_basemean;
    47 vector<float> drs_gainmean;
    48 vector<float> drs_triggeroffsetmean;
    49 
    50 size_t FOpenDataFile(
    51         const char *datafilename,               // path to fits file containing FACT raw data
    52         fits * * datafile,                              // pointer to pointer, where to return the fits object
    53         vector<int16_t> &Data,                  // vector, which will be filled with raw data
    54         vector<int16_t> &StartCells,    // vector, which will be filled with DRS start positions
    55         unsigned int &EventID,                  // index of the current event
    56         UInt_t &RegionOfInterest,               // Width of the Region, read out of the DRS
    57         UInt_t &NumberOfPixels,                 // Total number of pixel, read out of the camera
    58         size_t &PXLxROI,                                // Size of column "Data" = #Pixel x ROI
    59                         // this can be used, to x-check RegionOfInterest and NumberOfPixels
    60         int VerbosityLevel=1
    61 );
     47size_t TriggerOffsetROI, RC;
     48vector<float> Offset, Gain, TriggerOffset;
    6249
    6350
     
    7360
    7461
    75 float getValue( int, int );
    76 void  computeN1mean( int );
    77 void  removeSpikes( int );
    7862
    7963// histograms
     
    10791int searchSinglesPeaks = 0;
    10892
    109 
    11093int fbsl(
    11194        const char *datafilename                = "path-to-datafile.fits.gz",
     
    117100        int firstpixel          = 0,
    118101        int npixel                              = -1,
    119         bool produceGraphic = false
     102        bool produceGraphic = true
    120103){
    121104        fits * datafile = NULL;
    122         NEvents =  FOpenDataFile(
    123                         datafilename,
    124                         &datafile,
    125                         Data,
    126                         StartCells,
    127                         EventID,
    128                         RegionOfInterest,
    129                         NumberOfPixels,
    130                         PXLxROI);
     105
     106        NEvents = openDataFits(
     107                datafilename,
     108                &datafile,
     109                Data,
     110                StartCells,
     111                EventID,
     112                RegionOfInterest,
     113                NumberOfPixels,
     114                PXLxROI);
    131115
    132116    printf("number of events in file: %d\n", NEvents);
     117        if (NEvents == 0){
     118                cout << "return code of openDataFits:" << datafilename<< endl;
     119                cout << "is zero -> aborting." << endl;
     120                return 1;
     121        }
     122
    133123
    134124    // compare the number of events in the data file with the nevents the
    135125    // the user would like to read. nevents = -1 means: read all
    136126    if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;
    137        
    138                 if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels;
    139 
    140         FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
     127        if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels;
     128
     129        RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI);
     130        if (RC == 0){
     131                cout << "return code of openCalibFits:" << drsfilename << endl;
     132                cout << "is zero -> aborting." << endl;
     133                return 1;
     134        }
     135
    141136
    142137  BookHistos( RegionOfInterest, npixel );
     
    151146                }
    152147               
    153         // loop over pixel
    154         for ( int pix = firstpixel ; pix < npixel+firstpixel ; pix++ ){
    155 
    156             // loop over DRS slices
    157             for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
    158                 Ameas[ sl ] = getValue(sl, pix);
    159             }
    160 
    161             computeN1mean( RegionOfInterest ); // prepare spike removal
    162             removeSpikes (RegionOfInterest );  // output in Vcorr
    163 
    164             // filter Vcorr with sliding average using FIR filter function
    165                                                 // 8 is here the HalfWidth of the sliding average window. 
    166                                                 sliding_avg(Vcorr, Vslide, 8);
    167 
    168                                                 //            factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
     148                // loop over pixel
     149                for ( int pix = firstpixel ; pix < npixel+firstpixel ; pix++ ){
     150
     151
     152                        // get the data of this pixel from the  Data  vector
     153                        // apply the Drs Calibration and cut off 12 slices at the beginning
     154                        // and at the end.
     155                        applyDrsCalibration( Ameas,pix,12,12,
     156                                Offset, Gain, TriggerOffset,
     157                                RegionOfInterest, Data, StartCells);
     158
     159                        // finds spikes in the raw data, and interpolates the value
     160                        // spikes are: 1 or 2 slice wide, positive non physical artifacts
     161                        removeSpikes (Ameas, Vcorr);
     162
     163                        // filter Vcorr with sliding average using FIR filter function
     164                        // 8 is here the HalfWidth of the sliding average window.
     165                        sliding_avg(Vcorr, Vslide, 8);
    169166
    170167            for ( unsigned int sl = 0; sl <RegionOfInterest ; sl++){
     
    215212                }
    216213        return( 0 );
    217 }
    218 
    219 void removeSpikes(int Samples){
    220 
    221     const float fract = 0.8;
    222     float x, xp, xpp;
    223 
    224     // assume that there are no spikes
    225     for ( int i = 0; i <  Samples; i++) Vcorr[i] = Ameas[i];
    226 
    227     // find the spike and replace it by mean value of neighbors
    228     for ( int i = 2; i < Samples-2 ; i++) {
    229 
    230     x = Ameas[i] - N1mean[i];
    231 
    232         if ( x < -5. ){ // a spike candidate
    233             // check consistency with a single channel spike
    234             xp = Ameas[i+1] - N1mean[i+1];
    235             xpp = Ameas[i+2] - N1mean[i+2];
    236 
    237             if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
    238                 // printf("double spike candidate\n");
    239                 Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
    240                 Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
    241                 i = i + 3;
    242             }
    243             else{
    244    
    245                 if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
    246                     Vcorr[i+1] = N1mean[i+1];
    247                     // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
    248                     // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
    249                     N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
    250                     i = i + 2;//do not care about the next sample it was the spike
    251                 }
    252                 // treatment for the end of the pipeline must be added !!!
    253             }
    254         }
    255         else{
    256              // do nothing
    257         }
    258     } // end of spike search and correction
    259 }
    260 
    261 void computeN1mean( int Samples ){
    262 // compute the mean of the left and right neighbors of a channel
    263 
    264     for( int i = 2; i < Samples - 2; i++){
    265 /*        if (i == 0){ // use right sample es mean
    266             N1mean[i] = Ameas[i+1];
    267         }
    268         else if ( i == Samples-1 ){ //use left sample as mean
    269             N1mean[i] = Ameas[i-1];
    270         }
    271         else{
    272             N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
    273         }
    274 */
    275         N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
    276     }
    277 } // end of computeN1mean computation
    278 
    279 float getValue( int slice, int pixel ){
    280 
    281     const float dconv = 2000/4096.0;
    282 
    283     float vraw, vcal;
    284 
    285     unsigned int pixel_pt;
    286     unsigned int slice_pt;
    287     unsigned int cal_pt;
    288     unsigned int drs_cal_offset;
    289 
    290     // printf("pixel = %d, slice = %d\n", slice, pixel);
    291 
    292     pixel_pt = pixel * RegionOfInterest;
    293     slice_pt = pixel_pt + slice;
    294     drs_cal_offset = ( slice + StartCells[ pixel ] )%RegionOfInterest;
    295     cal_pt    = pixel_pt + drs_cal_offset;
    296 
    297     vraw = Data[ slice_pt ] * dconv;
    298     vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
    299 
    300     return( vcal );
    301214}
    302215
     
    375288} // end of function: void ana::SaveHistograms( void )
    376289
    377 size_t FOpenDataFile(
    378         const char *datafilename,               // path to fits file containing FACT raw data
    379         fits * * datafile,                              // ptr to pointer, where to return the fits object
    380         vector<int16_t> &Data,                  // vector, which will be filled with raw data
    381         vector<int16_t> &StartCells,    // vector, which will be filled with DRS start positions
    382         unsigned int &EventID,                  // index of the current event
    383         UInt_t &RegionOfInterest,               // Width of the Region, read out of the DRS
    384         UInt_t &NumberOfPixels,                 // Total number of pixel, read out of the camera
    385         size_t &PXLxROI,                                // Size of column "Data" = #Pixel x ROI
    386                         // this can be used, to x-check RegionOfInterest and NumberOfPixels
    387         int VerbosityLevel                              //
    388 ) {
    389         size_t NumberOfEvents;
    390         *datafile = new fits(datafilename);
    391         if (!(*(*datafile))) {
    392                 if (VerbosityLevel > 0)
    393                         cout << "Couldn't properly open the datafile: " << datafilename << endl;
    394                 return 0;
    395         }
    396 
    397         NumberOfEvents = (*datafile)->GetNumRows();
    398         if (NumberOfEvents < 1){
    399                 if (VerbosityLevel > 0){
    400                         cout << "Warning in FOpenDataFile of file: " << datafilename << endl;
    401                         cout << "the file contains no events." << endl;
    402                 }
    403         }
    404 
    405         RegionOfInterest = (*datafile)->GetUInt("NROI");
    406         NumberOfPixels = (*datafile)->GetUInt("NPIX");
    407         PXLxROI = (*datafile)->GetN("Data");
    408 
    409         if ( RegionOfInterest * NumberOfPixels != PXLxROI) // something in the file went wrong
    410         {
    411                 if (VerbosityLevel > 0){
    412                         cout << "Warning in FOpenDataFile of file: " << datafilename << endl;
    413                         cout << "RegionOfInterest * NumberOfPixels != PXLxROI" << endl;
    414                         cout << "--> " << RegionOfInterest;
    415                         cout << " * " << NumberOfPixels;
    416                         cout << " = " << RegionOfInterest * NumberOfPixels;
    417                         cout << ", but PXLxROI =" << PXLxROI << endl;
    418                 }
    419                 return 0;
    420         }
    421 
    422         //Set the sizes of the data vectors
    423         Data.resize(PXLxROI, 0);
    424         StartCells.resize(NumberOfPixels, 0);
    425 
    426         //Link the data to variables
    427         (*datafile)->SetRefAddress("EventNum", EventID);
    428         (*datafile)->SetVecAddress("Data", Data);
    429         (*datafile)->SetVecAddress("StartCellData", StartCells);
    430 
    431         return NumberOfEvents;
    432 }
    433 
Note: See TracChangeset for help on using the changeset viewer.