source: fact/tools/rootmacros/example.C@ 17963

#warning I AM NOT YET TESTED ... DO NOT USE ME

// example.C
// ROOT macro to show, how FACT raw data can be read into ROOT
// -> be DRS calibrated using a dedicated .drs.fits file
// -> also the famous DRS spikes will be removed.

#include <TROOT.h>      // everybody needs the ROOT of everything
#include <TStyle.h>     // in order to set gStyle->SetPallette(1,0);
#include <TCanvas.h>    // well ..we are going to show some plots, right?
#include <TTimer.h>     // dunno
#include <Getline.h>    // well ... for Getline()
#include <TFile.h>      // for root file handling
#include <TH1F.h>       // we will use these Histograms

#define HAVE_ZLIB       // -> Info for fits.h
#include "fits.h"       // fits.h was written by Thomas Bretz

// provides 2 functions to get easily the raw data
// out of our fits files... there is more data in the fits files,
// but it is not needed in this example.
// For further Info, .. ask an expert, or wait for the next Example :-)
#include "openFits.h"   //
#include "openFits.c"   //

// DrsCalibration provides a function, which can be used to retrieve
// DRS-calibrated data of a certain pixel.
#include "DrsCalibration.h"
#include "DrsCalibration.C"
#include "SpikeRemoval.h"
#include "SpikeRemoval.C"

// contains implementation of std FIR filters for vector<float> and
// a simple sldig average filter, which is not FIR.
#include "factfir.C"

// These variables will be linked to the raw data file.
// So when ever the Memberfunction "GetRow(int i)" of a fits file is called
// These variables will contain the data of the i'th Event:
int NEvents;                        // total number of events in file
vector<int16_t> AllPixelDataVector; // All raw data of all pixel of the current Event
vector<int16_t> StartCellVector;    // All StartCells of all pixel of the current Event
unsigned int CurrentEventID;        // The id = upcounting number of the current Event
size_t PXLxROI;                     // the width of data read out of the DRS chips x the number of Pixels
                                    //      so this will be the length of AllPixelDataVector
UInt_t RegionOfInterest;            // The region of Interest of all pixel in the current file
UInt_t NumberOfPixels;              // the number of pixels in the current file
// please note:
// Most of this information is not actually needed to plot or work with
// the raw data. It is just needed to apply the DRS calibration.
// So in principle the user of this script should not even see, that these
// data exists somewhere, but I was not yet able to hide it...



// these variables will contain the calibration constants read out of the
// dedicated DRS calibration files.
vector<float>   Offset;             // the offsets for all 1024 slices for all 1440 channels
vector<float>   Gain;               // the gains for all 1024 slices for all 1440 channels
vector<float>   TriggerOffset;      // a minor effect calibration constant. for details see DrsCalibration.h and .c
size_t          TriggerOffsetROI;   // the 'width' of the TriggerOffset per channel
size_t          RC;                 // RC = return code of openCalibFits() -> should be number of rows = 1 !!!

vector<float> Wavelet;              // here we store the data of one single event of one single pixel
TH1F* hWavelet = NULL;              // this we use for plotting the Wavelet
TCanvas *cExample = NULL;
TObjArray hList;
// The following two functions were introduced by Werner
// BookHistos() creates histograms and takes care of settings like
// Binning
// Naming
// Colors
// SaveHistos() takes care of Saving all Histos created in BookHistos to a root file
// this is a nice idea... makes 'sure' nothing is plotted, but forgotten to save...
void BookHistos();
void SaveHistograms(const char * );

int example(
    const char *datafilename        = "path-to-data/20111016_013.fits.gz",
    const char *drsfilename         = "path-to-data/20111016_011.drs.fits.gz",
    const char *OutRootFileName     = "path-to-analysis/20111016_013_example.root",
    int firstevent                  = 0,
    int nevents                     = -1,
    int firstpixel                  = 0,
    int npixel                      = -1,
    bool Debug                      = true,
    int verbosityLevel              = 1,
    bool ProduceGraphic             = true
        )
{
bool breakout = false;

cout << "I AM NOT YET TESTED ... DO NOT USE ME" << endl;


gStyle->SetPalette(1,0);    // nice color scheme
gROOT->SetStyle("Plain");   // don't know

    if ( Debug ){
        cExample = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
        cExample->Divide(1, 4);
    }

    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 ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
    if ( npixel == -1 || npixel > (int)NumberOfPixels  ) npixel = (int)NumberOfPixels; // -1 means all!
    if (verbosityLevel > 0) {
        cout <<"number of events in file: "<< NEvents << "\t";
        cout <<"of, which "<< nevents << "will be processed"<< endl;
        cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
        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;
    }

    BookHistos();

//-----------------------------------------------------------------------------
// Loops over Every Event and Pixel
//-----------------------------------------------------------------------------
    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++ ){
            if (verbosityLevel > 0){
                if (pix == firstpixel){
                    cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl;
                }
            }

            // this is the interesting point of this script...
            // well actually it is not interesting, but beginning from this point on
            // the user of this script has access to the pure data of:
            // event number ev
            // pixel number pix
            // The vector Wavewlet contains now the data, while the first 12 slices
            // and the last 12 slices were cut out of the raw data...
            // this was decided once on a FACT meeting.
            applyDrsCalibration( Wavelet, pix, 12, 12,
                                Offset, Gain, TriggerOffset,
                                RegionOfInterest, AllPixelDataVector, StartCellVector);

            // HERE we could plot the Wavelet, but lets first remove the spikes.

                        // finds spikes in the raw data, and interpolates the value
                        // spikes are: 1 or 2 slice wide, positive non physical artifacts
            // the data containing the spike is taken from 'Wavelet'
            // and the cleaned result is stored back to 'Wavelet'
            // synopsis: removeSpikes (input-vector, output-vector);
            removeSpikes (Wavelet, Wavelet);

            // Apply a sliding average Filter to the raw data ... it looks nicer for plotting
            // :-) This filter is defined in fir.h and .c, even though this particular Filter is not
            // a FIR filter, since it filters using past and future of a certain sample.
            // The filter width is 2*4+1=9 slices.
            sliding_avg(Wavelet, Wavelet, 4);


            if ( Debug ){
                hWavelet->GetXaxis()->Set(Wavelet.size() , -0.5 , Wavelet.size()-0.5);

                for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
                    hWavelet->SetBinContent(sl, Wavelet[sl]);
                }

                cExample->cd(1);
                gPad->SetGrid();
                hWavelet->Draw();
                cExample->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;
                }
            }// end of if(spikeDebug)

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

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

if (ProduceGraphic){
    cExample = new TCanvas("cExample ","An Example TCanvas",10,10,400,400);
    hWavelet->Draw();
    cExample->Modified();
    cExample->Update();
}
        SaveHistograms( OutRootFileName );
        return( 0 );
}

// booking and parameter settings for all histos
void BookHistos(){

    hWavelet = new TH1F("hWavelet", "to be changed", Wavelet.size(), -0.5, Wavelet.size()-0.5 );
    hWavelet->GetXaxis()->SetTitle( "time in slices" );
    hWavelet->GetYaxis()->SetTitle( "amplitude in mV" );
    hList.Add(hWavelet);
}

void SaveHistograms( const char * loc_fname){
        // create the histogram file (replace if already existing)
        TFile tf( loc_fname, "RECREATE");

        hList.Write(); // write the major histograms into the top level directory

        tf.Close(); // close the file
}