source: fact/tools/rootmacros/example_new.C@ 13782

#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 "openFits.h"
#include "openFits.c"

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

#include "FOpenDataFile.h"
#include "FOpenDataFile.c"

#include "DrsCalibration.C"
#include "DrsCalibration.h"

#include "SpikeRemoval.h"
#include "SpikeRemoval.C"

bool breakout=false;
int verbosityLevel = 1;

int NEvents;
vector<int16_t> AllPixelDataVector;
vector<int16_t> StartCellVector;
unsigned int CurrentEventID;
size_t PXLxROI;
UInt_t RegionOfInterest;
UInt_t NumberOfPixels;

size_t TriggerOffsetROI, RC;
vector<float> Offset, Gain, TriggerOffset;

vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
vector<float> Vcorr(FAD_MAX_SAMPLES);  // corrected Values
vector<float> Vslide(FAD_MAX_SAMPLES);  // sliding average result


TH1F *waveform = NULL;
TCanvas *funnycanvas = NULL;
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// read FACT raw data
//      * remove spikes
//      * calculate baseline
//      * find peaks (CFD and normal discriminator)
//      * compute pulse height and pulse integral spektrum of the peaks
int example_new(
        char *datafilename              = "datafile.fits.gz",
        const char *drsfilename = "calibfile.drs.fits.gz",
        int firstevent                  = 0,
        int nevents                     = -1,
        int firstpixel                  = 0,
        int npixel                      = -1
)
{
        gStyle->SetPalette(1,0);
        gROOT->SetStyle("Plain");
    
    // create a new TH1F for displaying the wave form ... the binning will be set again, just before we use it
    // so the binning is here just a dummy setting...
    waveform = new TH1F("waveform", "content of current pixel in current Event", 1024,-0.5,1023.5);
    funnycanvas = new TCanvas("funnycanvas","funny title", 0 , 0 ,300,300);

        fits * datafile;
        // Opens the raw data file and 'binds' the variables given as
        // Parameters to the data file. So they are filled with
        // raw data as soon as datafile->GetRow(int) is called.
        NEvents = openDataFits( datafilename, &datafile,
                AllPixelDataVector, StartCellVector, CurrentEventID,
                RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
        if (NEvents == 0){
                cout << "return code of OpenDataFile:" << datafilename<< endl;
                cout << "is zero -> aborting." << endl;
                return 1;
        }

        if (verbosityLevel > 0)
                cout <<"number of events in file: "<< NEvents << "\t";
        if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
        if (verbosityLevel > 0)
                cout <<"of, which "<< nevents << "will be processed"<< endl;

        if (verbosityLevel > 0)
                cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
        if ( npixel == -1 || npixel > (int)NumberOfPixels  ) npixel = (int)NumberOfPixels; // -1 means all!
        if (verbosityLevel > 0)
                cout <<"of, which "<< npixel << "will be processed"<< endl;

        RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI);
        if (RC == 0){
                cout << "return code of openCalibFits:" << drsfilename << endl;
                cout << "is zero -> aborting." << endl;
                return 1;
        }

    //-----------------------------------------------------------------------------
    // Loops over Every Event and Pixel
    //-----------------------------------------------------------------------------
        cout << "Processing Events in total: " << nevents << endl;

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

                for ( int pix = firstpixel; pix < firstpixel+npixel; pix++ ){
        
                        applyDrsCalibration( Ameas,pix,12,12,
                                Offset, Gain, TriggerOffset,
                                RegionOfInterest, AllPixelDataVector, StartCellVector);

                        // 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
                        sliding_avg(Vcorr, Vslide, 8);

                                        
                        waveform->GetXaxis()->Set(Vslide.size() , -0.5 , Vslide.size()-0.5);

                        for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
                                waveform->SetBinContent( sl, Vslide[sl] );
                        }
            
            
            waveform->Draw();
            funnycanvas->Modified();
            funnycanvas->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;
            }

            if (breakout) break;
                } // end of loop over pixels
                if (breakout) break;
        }       // end of loop over pixels
        
        cout << "Event processing done." << endl;  

        return 0;
}