source: fact/tools/rootmacros/fbsl.C@ 12385

#include <TROOT.h>
#include <TCanvas.h>
#include <TProfile.h>
#include <TTimer.h>
#include <TH1F.h>
#include <TH2F.h>
#include <Getline.h>
#include <TLine.h>
#include <TMath.h>
#include <TFile.h>

#include <stdio.h>
#include <stdint.h>
#include <cstdio>

#define HAVE_ZLIB
#include "fits.h"
#include "FOpenCalibFile.c"

#include "zerosearch.C"
#include "factfir.C"

#define NPIX  1440
#define NCELL 1024

// data access and treatment
#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;
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
);


vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
vector<float> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors
vector<float> Vcorr(FAD_MAX_SAMPLES);  // corrected Values
vector<float> Vdiff(FAD_MAX_SAMPLES);  // numerical derivative

vector<float> Vslide(FAD_MAX_SAMPLES);  // sliding average result
vector<float> Vcfd(FAD_MAX_SAMPLES);    // CDF result
vector<float> Vcfd2(FAD_MAX_SAMPLES);    // CDF result + 2nd sliding average


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 hPixelCellData(
    "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
TH1F *hMeanBsl, *hpltMeanBsl;
TH1F *hRmsBsl, *hpltRmsBsl; 
TObjArray hList;
TObjArray hListBaseline;

void BookHistos( int , int);
void SaveHistograms( const char * );

// Create a canvas
TCanvas* CW;
TCanvas* cFilter;

int spikeDebug = 0;
int searchSinglesPeaks = 0;


int fbsl( 
        const char *datafilename                = "path-to-datafile.fits.gz",
        const char *drsfilename = "path-to-calibfile.drs.fits.gz",
        const char *TextOutFileName = "./appendfile.txt",
        const char *RootOutFileName = "./datafile.root",
        int firstevent                  = 0, 
        int nevents                     = -1, 
        int firstpixel          = 0,
        int npixel                              = -1,
        bool produceGraphic = false
){
        fits * datafile = NULL;
        NEvents =  FOpenDataFile( 
                        datafilename, 
                        &datafile, 
                        Data, 
                        StartCells, 
                        EventID, 
                        RegionOfInterest, 
                        NumberOfPixels, 
                        PXLxROI);

    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;
        
                if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels;

        FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);

  BookHistos( RegionOfInterest, npixel );
    
                // loop over events
    for ( int ev = firstevent; ev < firstevent + nevents; ev++) {

            datafile->GetRow( ev );
                
                if ( ev % 100 == 0 ){
                        cout << "Event ID: " << EventID << endl;
                }
                
        // 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);

            for ( unsigned int sl = 0; sl <RegionOfInterest ; sl++){
               // hPixelCellData.Fill( sl, Vcorr[sl] );
               hBaseline[pix-firstpixel]->Fill( Vslide[sl] ) ;
            }   
            }
    }

    FILE *fp;
    TString fName; 
    fName = TextOutFileName;

    fp = fopen(fName, "a+");
    fprintf( fp, "FILE: %s\n", datafilename );
    fprintf( fp, "NEVENTS: %d\n", nevents);
    NBSLeve = nevents; // this has to be changed when computation on a subset of a run is implemented
    fprintf( fp, "NBSLEVE:  %d\n", NBSLeve );

    for (int pix = firstpixel; pix < firstpixel+npixel; pix++) {

        float vmaxprob = hBaseline[pix]->GetXaxis()->GetBinCenter(
            hBaseline[pix-firstpixel]->GetMaximumBin() );

        fprintf( fp, "%8.3f", vmaxprob );

        hMeanBsl->Fill( vmaxprob );
        hpltMeanBsl->SetBinContent(pix+1, vmaxprob );

        hRmsBsl->Fill(hBaseline[pix-firstpixel]->GetRMS() );
        hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() );   
    }
    fprintf( fp, "\n" );

    fclose( fp );

    SaveHistograms( RootOutFileName );
                if (produceGraphic){
        TCanvas * cMeanBsl = new TCanvas();
        cMeanBsl->cd();
        hMeanBsl->Draw();
        cMeanBsl->Update();

        TCanvas * cRmsBsl = new TCanvas();
        cRmsBsl->cd();
        hRmsBsl->Draw();
        cMeanBsl->Update();
                }
        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 );
}

void BookHistos( int Samples , int NumberOfPixel){
// booking and parameter settings for all histos

    // histograms for baseline extraction
    char hName[500];
    char hTitle[500];

    TH1F *h;
    
    printf("inside BookHistos\n");

    for( int i = 0; i < NumberOfPixel; i++ ) {

        // printf("call sprintf [%d]\n", i );
        sprintf(&hTitle[0],"all events all slices of pixel %d", i);
        sprintf(&hName[0],"base%d", i);
        // printf("call sprintf [%d] done\n", i );

        h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );

        // printf("histo booked\n");
        h->GetXaxis()->SetTitle( "Sample value (mV)" );
        h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
        // printf("histo title set\n");
        hListBaseline.Add( h );
        // printf("histo added to List\n");
        hBaseline[i] = h;
        // printf("histo assigned to array\n");
    }

    printf("made HBaseline * 1440\n");

    hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
    hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
    hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
    hList.Add( hMeanBsl );

    hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
    hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
    hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
    hList.Add( hpltMeanBsl );

    hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
    hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
    hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
    hList.Add( hRmsBsl );

    hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
    hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
    hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
    hList.Add( hpltRmsBsl );
}


void SaveHistograms( const char * loc_fname ){

  TString fName; // name of the histogram file

  /* create the filename for the histogram file */
  fName = loc_fname; // use the name of the tree file
  //fName.Remove(fName.Length() - 5); // remove the extension .root
  //fName += "_histo.root"; // add the new extension
  //fName += ".root"; 

  TFile tf( fName, "RECREATE"); // create the histogram file (replace if already existing)

  hList.Write(); // write the major histograms into the top level directory 
  tf.mkdir("BaselineHisto"); tf.cd("BaselineHisto"); // go to new subdirectory
  hListBaseline.Write(); // write histos into subdirectory
  
  tf.Close(); // close the file

} // 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;
}