source: fact/tools/rootmacros/fpeak_cdf.C@ 12363

#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 <TBox.h>
#include <TMath.h>
#include <TFile.h>
#include <TStyle.h>

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

#define NPIX  1440
#define NCELL 1024
#define FAD_MAX_SAMPLES 1024

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

#include "discriminator.h"
#include "discriminator.C"
#include "zerosearch.h"
#include "zerosearch.C"
#include "factfir.C"


vector<int16_t> AllPixelDataVector;
vector<int16_t> StartCellVector;

unsigned int CurrentEventID;

bool breakout=false;

size_t ROIxNOP;
UInt_t NumberOfPixels;
UInt_t RegionOfInterest;
int NEvents;

size_t drs_n;
vector<float> drs_basemean;
vector<float> drs_gainmean;
vector<float> drs_triggeroffsetmean;

int FOpenDataFile( fits & );


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> 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 NumberOfDebugHistoTypes = 7;
const unsigned int
        Ameas_  = 0,
        N1mean_ = 1,
        Vcorr_  = 2,
        Vtest_  = 3,
        Vslide_ = 4,
        Vcfd_   = 5,
        Vcfd2_  = 6;

TH1F* h;
TH1F *debugHistos;
TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
                    // x = pixel id, y = DRS cell id

TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
TH1F *hMeanBsl, *hpltMeanBsl;
TH1F *hRmsBsl, *hpltRmsBsl; 
TH2F * hAmplSpek_cfd;
TH2F * hAmplSpek_discri;
TH2F * hTemp_Array[1441];
TObjArray hList;
TObjArray hListBaseline, hListTemplates;

void BookHistos( );
void SaveHistograms( char * );

int searchSinglesPeaks = 0;

int fpeak(
        char *datafilename              = "../../20111011_055.fits.gz",
        const char *drsfilename = "../../20111011_054.drs.fits.gz",
        int nevents                     = -1,
        int firstevent                  = 0,
        int PixelID                             = -1,
        bool spikeDebug = false,
        int verbosityLevel = 1 // different verbosity levels can be implemented here
        )
{

// Create (pointer to) Canvases, which are used in every run,
// also in 'non-debug' runs
        TCanvas * cSpektrum;
        TCanvas * cStartCell;
        TCanvas *cTemplate;
        cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400);
        cSpektrum->Divide(1,2);
        cStartCell = new TCanvas("cStartCell ", "The Startcells of this run", 10,410,400,400);
        cTemplate = new TCanvas("cTemplate","Template of current Pixel",1,1,1600,1000);

        // Canvases only need if spike Debug, but I want to deklare
        // the pointers anyway ...
        TCanvas *cRawAndSpikeRemoval = NULL;
        TCanvas *cFiltered = NULL;


        if (spikeDebug){
                cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",410,10,400,400);
                cRawAndSpikeRemoval->Divide(1, 2);
                cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
                cFiltered->Divide(1, 2);
        }

gStyle->SetPalette(1,0);
gROOT->SetStyle("Plain");
// read FACT raw data 
//      * remove spikes
//      * calculate baseline
//      * find peaks (CFD and normal discriminator)
//      * compute pulse height and pulse integral spektrum of the peaks

// sliding window filter settings
//      int k_slide = 16;
//      vector<double> a_slide(k_slide, 1);
//      double b_slide = k_slide;

// CFD filter settings
        int k_cfd = 10;
        vector<double> a_cfd(k_cfd, 0);
        double b_cfd = 1.;
        a_cfd[0]=-0.75;
        a_cfd[k_cfd-1]=1.;

// 2nd slinding window filter
//      int ks2 = 16;
//      vector<double> as2(ks2, 1);
//      double bs2 = ks2;
        
// Open the data file
        fits *datafile = new fits( datafilename );
        if (!datafile) {
            printf( "Could not open the file: %s\n", datafilename );
                  return 1;
        }
        
// access data
        NEvents = FOpenDataFile( *datafile );
        printf("number of events in file: %d\n", NEvents);
  if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!

//Get the DRS calibration
        FOpenCalibFile( drsfilename,
                                        drs_basemean,
                                        drs_gainmean,
                                        drs_triggeroffsetmean,
                                        drs_n);

//Check the sizes of the data columns
        if (drs_n != 1474560) {
                cerr << "error: DRS calib file has wrong ...erm...size ... drs_n is: "
                        << drs_n << endl;
                cerr << " Aborting." << endl;
                return 1;
        }

        if(ROIxNOP != 1474560)
        {
                cout << "warning: data_n should better be 1440x1024=1474560, but is "
                        << ROIxNOP << endl;
                cout << "this script is not guaranteed to run under these "
                        <<" circumstances....any way ... it is never guaranteed." << endl;
        }
// Book the histograms

        BookHistos( );

        float calibratedVoltage;

        for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
                // Get an Event --> consists of 1440 Pixel ...erm....data
                datafile->GetRow( ev );

                for ( int pix = 0; pix < 1440; pix++ ){

                        hStartCell->Fill( pix, StartCellVector[pix] );

                        // this is a stupid hack ... there is more code at the
                        // end of this loop to complete this hack ...
                        // beginning with if (PixelID != -1) as well.
                        if (PixelID != -1) {
                                pix = PixelID;
                                if (verbosityLevel > 0){
                                        cout << "Processing Event number: " << CurrentEventID << "\t"
                                                << "Pixel number: "<< pix << endl;
                                }
                        }

                        if (verbosityLevel > 0){
                                if (pix % 20 ==0){
                                        cout << "Processing Event number: " << CurrentEventID << "\t"
                                                << "Pixel number: "<< pix << endl;
                                }
                        }

                        // compute the DRs calibrated values and put them into Ameas[]
                        for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
                                calibratedVoltage = getValue( sl, pix);
                                if (verbosityLevel > 10){
                                        printf("calibratedVoltage = %f\n", calibratedVoltage);
                                }
                                Ameas[ sl ] = calibratedVoltage;

                                // in case we want to plot this ... we need to put it into a
                                // Histgramm
                                if (spikeDebug){
                                        debugHistos[Ameas_].SetBinContent(sl, calibratedVoltage);
                                }
                        }
                        // operates on Ameas[] and writes to N1mean[]
                        computeN1mean( RegionOfInterest );

                        // operates on Ameas[] and N1mean[], then writes to Vcorr[]
                        removeSpikes( RegionOfInterest );

                        // filter Vcorr with sliding average using FIR filter function
                        //factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
                        sliding_avg(Vcorr, Vslide, 8);
                        // filter Vslide with CFD using FIR filter function
                        factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
                        // filter Vcfd with sliding average using FIR filter function
                        //factfir(bs2 , as2, ks2, Vcfd, Vcfd2);
                        sliding_avg(Vcfd, Vcfd2, 8);


                        // peaks in Ameas[] are found by searching for zero crossings
                        // in Vcfd2
                        // first Argument 1 means ... *rising* edge
                        // second Argument 1 means ... search with stepsize 1 ... 8 is okay as well
                        vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 8);
                        // zXings means "zero cross ings"
                        //ShiftRegionBy(*zXings, -ks2/2);
                        EnlargeRegion(*zXings, 10, 10);
                        findAbsMaxInRegions(*zXings, Vslide);
                        removeMaximaBelow( *zXings, 11.0, 0);
                        removeMaximaAbove( *zXings, 14.0, 0);

                        if (zXings->size() != 0 ){
                                for (unsigned int i=0; i<zXings->size(); i++){
                                        if (verbosityLevel > 1){
                                                cout << zXings->at(i).maxPos << ":\t"
                                                        << zXings->at(i).maxVal <<endl;
                                        }
                                        hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);

                                        for (int j=-100; j<150; j++){
                                                if (zXings->at(i).maxPos + j >= 0 && zXings->at(i).maxPos + j <= 1023)
                                                        hTemp_Array[pix]->Fill(j+101, Vslide[zXings->at(i).maxPos + j]);
                                        }
                                }
                        }

                        if ( spikeDebug ){
                                for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
                                        debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
                                        debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
                                        debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
                                }

                                cRawAndSpikeRemoval->cd( 1);
                                gPad->SetGrid();
                                debugHistos[Ameas_].Draw();

                                cRawAndSpikeRemoval->cd( 2);
                                gPad->SetGrid();
                                debugHistos[Vcorr_].Draw();

                                cFiltered->cd(1);
                                gPad->SetGrid();
                                debugHistos[Vslide_].Draw();

                                TBox *OneBox;
                                vector<TBox*> MyBoxes;
                                for (unsigned int i=0; i<zXings->size(); i++){
                                        OneBox = new TBox(
                                                zXings->at(i).maxPos -10 ,
                                                zXings->at(i).maxVal -0.5,
                                                zXings->at(i).maxPos +10 ,
                                                zXings->at(i).maxVal +0.5);
                                        OneBox->SetLineColor(kBlue);
                                        OneBox->SetLineWidth(1);
                                        OneBox->SetFillStyle(0);
                                        OneBox->SetFillColor(kRed);
                                        MyBoxes.push_back(OneBox);
                                        OneBox->Draw();
                                }

                                cFiltered->cd(2);
                                gPad->SetGrid();
                                debugHistos[Vcfd2_].Draw();
                                TLine *zeroline = new TLine(0, 0, 1024, 0);
                                zeroline->SetLineColor(kBlue);
                                zeroline->Draw();

                                cRawAndSpikeRemoval->Update();
                                cFiltered->Update();

                                cTemplate->cd();
                                hTemp_Array[pix]->Draw("COLZ");
                                cTemplate->Modified();
                                cTemplate->Update();

                                //Process gui events asynchronously during input
                                TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
                                timer.TurnOn();
                                TString input = Getline("Type 'q' to exit, <return> to go on: ");
                                timer.TurnOff();
                                if (input=="q\n") {
                                        breakout=true;
                                }

                                //TODO!!!!!!!!!
                                // do some Garbage collection ...
                                // all the Objects on the heap should be deleted here.

                        }// end of if(spikeDebug)

                        delete zXings;

                        // this is the 2nd part of the ugly hack in the beginning.
                        // this makes sure, that the loop ends
                        if (PixelID != -1){
                                pix = 1440;
                        }


                        if (breakout)
                                break;

                } // end of loop over pixels

if (ev % 10 == 0){
                cSpektrum->cd(1);
                hAmplSpek_cfd->Draw("COLZ");
                cSpektrum->cd(2);
                hAmplSpek_discri->Draw("COLZ");
                cSpektrum->Modified();
                cSpektrum->Update();

                // updating seems not to work ..
                // debug cout
                cStartCell->cd();
                hStartCell->Draw();
                cStartCell->Modified();
                cStartCell->Update();

                cTemplate->cd();
                hTemp_Array[PixelID]->GetYaxis()->SetRangeUser(-10,40);
                hTemp_Array[PixelID]->Draw("COLZ");
                cTemplate->Modified();
                cTemplate->Update();
}


                if (breakout)
                        break;
        }       // end of loop over pixels


        SaveHistograms( datafilename );

        return( 0 );
}

void removeSpikes(int Samples){

    const float fract = 0.8;
    float x, xp, xpp, x3p;

    // 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 = 0; i < Samples; i++) {

    // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ 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];
            x3p = Ameas[i+3] - N1mean[i+3];

            // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p);

            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.;
                // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]);
                // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]);
                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
        for ( int i = 0; i < Samples; i++ ) debugHistos[ Vcorr_ ].SetBinContent( i, Vcorr[i] );
}

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 + StartCellVector[ pixel ] )%RegionOfInterest;
        cal_pt    = pixel_pt + drs_cal_offset;

        vraw = AllPixelDataVector[ slice_pt ] * dconv;
        vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;

        return( vcal );
}

// booking and parameter settings for all histos
void BookHistos( ){
        // histograms for baseline extraction
        char hName[500];
        char hTitle[500];

        TH1F *h;

        for( int i = 0; i < NPIX; i++ ) {
                sprintf(&hTitle[0],"all events all slices of pixel %d", i);
                sprintf(&hName[0],"base%d", i);

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

                h->GetXaxis()->SetTitle( "Sample value (mV)" );
                h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
                hListBaseline.Add( h );
                hBaseline[i] = h;
        }

    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 );

        hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",1440,-0.5,1439.5, 256, -27.5, 100.5);
        hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" );
        hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" );
        hList.Add( hAmplSpek_cfd );

        hAmplSpek_discri = new TH2F("hAmplSpek_discri","amplitude spektrum - std discriminator",1440,-0.5,1439.5, 256, -27.5, 100.5);
        hAmplSpek_discri->GetXaxis()->SetTitle( "pixel" );
        hAmplSpek_discri->GetXaxis()->SetTitle( "amplitude in mV" );
        hList.Add( hAmplSpek_discri );

        hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
        hStartCell->GetXaxis()->SetTitle( "pixel" );
        hStartCell->GetXaxis()->SetTitle( "slice" );
        hList.Add( hStartCell );

        debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
        for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
                debugHistos[ type ].SetBins(1024, 0, 1024);
                debugHistos[ type ].SetLineColor(1);
                debugHistos[ type ].SetLineWidth(2);

                // set X axis paras
                debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
                debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
                debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
                debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));

                // set Y axis paras
                debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
                debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
                debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
                debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
        }

        TH2F *temp;
        for ( int type = 0; type < 1441; type++){
                sprintf(&hTitle[0],"pulse template of pixel %d", type );
                sprintf(&hName[0],"template_%d", type );

                temp = new TH2F( hName, hTitle, 256, 0, 255, 256, -10.5, 117.5);

                temp->GetXaxis()->SetTitle( "Time slice (%.1f ns/slice)" );
                temp->GetYaxis()->SetTitle( "Amplitude (about mV)" );
                hListTemplates.Add( temp );
                hTemp_Array[ type ] = temp;
        }


}

void SaveHistograms( 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
        // TODO ... next statement doesn't work for ".fits.gz"
        fName.Remove(fName.Length() - 5); // remove the extension .fits
        fName += "_discri.root"; // add the new extension

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

        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.cd("..");
        tf.mkdir("PulseTempplates");
        tf.cd("PulseTempplates");
        hListTemplates.Write(); // write histos into subdirectory

        tf.Close(); // close the file
} // end of SaveHistograms( char * loc_fname )

int FOpenDataFile(fits &datafile){

//      cout << "-------------------- Data Header -------------------" << endl;
//      datafile.PrintKeys();
//      cout << "------------------- Data Columns -------------------" << endl;
//      datafile.PrintColumns();

        //Get the size of the data column
        RegionOfInterest        = datafile.GetUInt("NROI");
        NumberOfPixels          = datafile.GetUInt("NPIX");
        //Size of column "Data" = #Pixel x ROI
        ROIxNOP                         = datafile.GetN("Data");

        //Set the sizes of the data vectors
        AllPixelDataVector.resize(ROIxNOP,0);
        StartCellVector.resize(NumberOfPixels,0);

        //Link the data to variables
        datafile.SetRefAddress("EventNum",  CurrentEventID);
        datafile.SetVecAddress("Data", AllPixelDataVector);
        datafile.SetVecAddress("StartCellData", StartCellVector);
        datafile.GetRow(0);

        return datafile.GetNumRows() ;
}