Ignore:
Timestamp:
11/04/11 20:18:57 (13 years ago)
Author:
neise
Message:
moved spike removal out of this file
File:
1 edited

Legend:

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

    r12371 r12379  
    3737#include "DrsCalibration.h"
    3838
     39#include "SpikeRemoval.h"
     40#include "SpikeRemoval.C"
     41
    3942bool breakout=false;
    4043
     
    5356
    5457vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
    55 vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
    5658vector<float> Vcorr(FAD_MAX_SAMPLES);  // corrected Values
    57 vector<float> Vdiff(FAD_MAX_SAMPLES);  // numerical derivative
    58 
    5959vector<float> Vslide(FAD_MAX_SAMPLES);  // sliding average result
    6060vector<float> Vcfd(FAD_MAX_SAMPLES);    // CDF result
    6161vector<float> Vcfd2(FAD_MAX_SAMPLES);    // CDF result + 2nd sliding average
    62 
    63 
    64 //float getValue( int, int );
    65 void computeN1mean( int );
    66 void removeSpikes( int );
    6762
    6863// histograms
     
    130125
    131126        fits * datafile;
     127        // Opens the raw data file and 'binds' the variables given as
     128        // Parameters to the data file. So they are filled with
     129        // raw data as soon as datafile->GetRow(int) is called.
    132130        NEvents = OpenDataFile( datafilename, &datafile,
    133131                AllPixelDataVector, StartCellVector, CurrentEventID,
     
    177175                                RegionOfInterest, AllPixelDataVector, StartCellVector);
    178176
    179                         // operates on Ameas[] and writes to N1mean[]
    180                         computeN1mean( RegionOfInterest );
    181 
    182                         // operates on Ameas[] and N1mean[], then writes to Vcorr[]
    183                         removeSpikes( RegionOfInterest );
     177                        // finds spikes in the raw data, and interpolates the value
     178                        // spikes are: 1 or 2 slice wide, positive non physical artifacts
     179                        removeSpikes (Ameas, Vcorr);
    184180
    185181                        // filter Vcorr with sliding average using FIR filter function
     
    288284}
    289285
    290 void removeSpikes(int Samples){
    291 
    292     const float fract = 0.8;
    293     float x, xp, xpp, x3p;
    294 
    295     // assume that there are no spikes
    296     for ( int i = 0; i <  Samples; i++) Vcorr[i] = Ameas[i];
    297 
    298 // find the spike and replace it by mean value of neighbors
    299     for ( int i = 0; i < Samples; i++) {
    300 
    301 
    302     x = Ameas[i] - N1mean[i];
    303 
    304         if ( x < -5. ){ // a spike candidate
    305             // check consistency with a single channel spike
    306             xp = Ameas[i+1] - N1mean[i+1];
    307             xpp = Ameas[i+2] - N1mean[i+2];
    308             x3p = Ameas[i+3] - N1mean[i+3];
    309 
    310 
    311             if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
    312                 Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
    313                 Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
    314                 i = i + 3;
    315             }
    316             else{
    317    
    318                 if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
    319                     Vcorr[i+1] = N1mean[i+1];
    320                     N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
    321                     i = i + 2;//do not care about the next sample it was the spike
    322                 }
    323                 // treatment for the end of the pipeline must be added !!!
    324             }
    325         }
    326         else{
    327              // do nothing
    328         }
    329     } // end of spike search and correction
    330         for ( int i = 0; i < Samples; i++ ) debugHistos[ Vcorr_ ].SetBinContent( i, Vcorr[i] );
    331 }
    332 
    333 void computeN1mean( int Samples ){
    334 // compute the mean of the left and right neighbors of a channel
    335 
    336     for( int i = 2; i < Samples - 2; i++){
    337 /*        if (i == 0){ // use right sample es mean
    338             N1mean[i] = Ameas[i+1];
    339         }
    340         else if ( i == Samples-1 ){ //use left sample as mean
    341             N1mean[i] = Ameas[i-1];
    342         }
    343         else{
    344             N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
    345         }
    346 */
    347         N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
    348     }
    349 } // end of computeN1mean computation
    350 
    351 
    352286// booking and parameter settings for all histos
    353287void BookHistos( ){
Note: See TracChangeset for help on using the changeset viewer.