Index: fact/tools/rootmacros/example.C =================================================================== --- fact/tools/rootmacros/example.C (revision 12718) +++ fact/tools/rootmacros/example.C (revision 12718) @@ -0,0 +1,235 @@ +#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 // everybody needs the ROOT of everything +#include // in order to set gStyle->SetPallette(1,0); +#include // well ..we are going to show some plots, right? +#include // 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 and +// a simple sldig average filter, which is not FIR. +#include "factfir.h" +#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 AllPixelDataVector; // All raw data of all pixel of the current Event +vector 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 Offset; // the offsets for all 1024 slices for all 1440 channels +vector Gain; // the gains for all 1024 slices for all 1440 channels +vector 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 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; +// 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 = false, + int verbosityLevel = 1, + bool ProduceGraphic = true + ) +{ + + +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 (spikeDebug){ + cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400); + cFiltered->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, 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); + + hDom = new TH1F("hDom", "test", Wavelet.size(), -0.5, Wavelet.size()-0.5 ); + hDom->GetXaxis()->SetTitle( "time in slices" ); + hDom->GetYaxis()->SetTitle( "amplitude in mV" ); + hList.Add(hDom); + +} + +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 +}