Index: /fact/tools/rootmacros/README.txt =================================================================== --- /fact/tools/rootmacros/README.txt (revision 12355) +++ /fact/tools/rootmacros/README.txt (revision 12356) @@ -1,3 +1,34 @@ README of FACT svntools/rootmacros + fbsl.C ROOT Macro computing the baseline for each pixel + the function is declared as: + 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 + ) + + the baseline and its rms is calculated for each pixel based on the data, + given in the datafile. + the results are appended to a textfile + and the histograms ,which were used for caluculation, as well as some +overviews are stored in a root file. + the 4 ints: firstevent, nevents, firstpixel, npixel can be used to calculate +only for a subset + the last bool can be set to true, which will open 2 Canvases with overview +histograms. + + for automatic production of baseline analysis one can call this macro like +this: +e.g. + root -l -q fbsl.C++'("/data00/fact-construction/raw/2011/10/26/20111026_036.fits.gz", "/data00/fact-construction/raw/2011/10/26/20111026_031.drs.fits.gz", "./fbsl.txt", "./20111026_036.root")' + + + on the ISDC cluster, I used /opt/root5.18.x86_64/bin/root for testing. Index: /fact/tools/rootmacros/factfir.C =================================================================== --- /fact/tools/rootmacros/factfir.C (revision 12355) +++ /fact/tools/rootmacros/factfir.C (revision 12356) @@ -2,20 +2,13 @@ #include - -#ifndef FAD_MAX_SAMPLES - #warning FAD_MAX_SAMPLES not defined. defining: FAD_MAX_SAMPLES to 1024 - #define FAD_MAX_SAMPLES 1024 -#endif - // source vector is -// float Ameas[FAD_MAX_SAMPLES]; void factfir(double b, vector &a, int k, vector &source, vector &dest){ //dest.clear(); - for (int slice =0; slice < FAD_MAX_SAMPLES; slice++) { + for (int slice =0; slice < (int)source.size(); slice++) { float currentval = 0; for (int i=0; i < k ; i++){ - currentval += a[i] * source[ (slice - i + FAD_MAX_SAMPLES) % FAD_MAX_SAMPLES ]; + currentval += a[i] * source[ (slice - i + source.size()) % source.size() ]; } dest[slice] = currentval / b; Index: /fact/tools/rootmacros/fbsl.C =================================================================== --- /fact/tools/rootmacros/fbsl.C (revision 12355) +++ /fact/tools/rootmacros/fbsl.C (revision 12356) @@ -16,9 +16,8 @@ #define HAVE_ZLIB #include "fits.h" -//#include "TPKplotevent.c" -//#include "FOpenDataFile.c" #include "FOpenCalibFile.c" #include "zerosearch.C" +#include "factfir.C" #define NPIX 1440 @@ -28,12 +27,17 @@ #define FAD_MAX_SAMPLES 1024 -vector data; -vector data_offset; - -unsigned int data_num; - -size_t data_n; -UInt_t data_px; -UInt_t data_roi; +//vector data; +//vector data_offset; +//unsigned int data_num; +//size_t data_n; +//UInt_t data_px; +//UInt_t data_roi; + vector Data; // vector, which will be filled with raw data + vector 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; @@ -44,5 +48,16 @@ vector drs_triggeroffsetmean; -int FOpenDataFile( fits & ); +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 &Data, // vector, which will be filled with raw data + vector &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 +); @@ -57,5 +72,4 @@ vector Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average -#include "factfir.C" float getValue( int, int ); @@ -75,6 +89,4 @@ TH1F* h; -TH2F* hStartCell; // id of the DRS physical pipeline cell where readout starts - // x = pixel id, y = DRS cell id TH2F hPixelCellData( "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.); @@ -85,6 +97,6 @@ TObjArray hListBaseline; -void BookHistos( int ); -void SaveHistograms( char * ); +void BookHistos( int , int); +void SaveHistograms( const char * ); // Create a canvas @@ -97,40 +109,24 @@ int fbsl( - char *datafilename = "../../20110804_024.fits.gz", - const char *drsfilename = "../../20110804_023.drs.fits.gz", - int pixelnr = 0, + 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 ){ -// read and analyze FACT raw data - -// sliding window filter settings - int k_slide = 16; - vector a_slide(k_slide, 1); - double b_slide = k_slide; - -// CFD filter settings - int k_cfd = 10; - vector 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 as2(ks2, 1); - double bs2 = ks2; - gROOT->SetStyle("Plain"); - -//------------------------------------------- -// Open the file -//------------------------------------------- - fits datafile( datafilename ); - if (!datafile) { - printf( "Could not open the file: %s\n", datafilename ); - return 1; - } - - // access data - NEvents = FOpenDataFile( datafile ); + 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); @@ -139,64 +135,40 @@ // the user would like to read. nevents = -1 means: read all if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; - - - // FILE *pedFile; - // pedFile = fopen("t.txt","u"); - // fprintf(pedFile,"# Pedestal Data"); - // fclose( pedFile ); -//------------------------------------------- -//Get the DRS calibration -//------------------------------------------- + + if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels; FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n); - - - //Check the sizes of the data columns - if(drs_n!=data_n) - { - cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl; - return 1; - } - // Book the histograms - // printf("call BookHistos\n"); - BookHistos( data_roi ); - // printf("returned from BookHistos\n"); - // Loop over events - cout << "--------------------- Data --------------------" << endl; - - float value; + + BookHistos( RegionOfInterest, npixel ); - // loop over events -// for ( int ev = firstevent; ev < nevents; ev++) { + // loop over events for ( int ev = firstevent; ev < firstevent + nevents; ev++) { - datafile.GetRow( ev ); + datafile->GetRow( ev ); if ( ev % 100 == 0 ){ - cout << "Event number: " << data_num << endl; + cout << "Event ID: " << EventID << endl; } // loop over pixel - for ( int pix = 0; pix < 1440; pix++ ){ - - // hStartCell->Fill( pix, data_offset[pix] ); - + for ( int pix = firstpixel ; pix < npixel+firstpixel ; pix++ ){ + // loop over DRS slices - for ( unsigned int sl = 0; sl < data_roi; sl++){ - + for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){ Ameas[ sl ] = getValue(sl, pix); - } - // printf("Ameas[100] = %f\n", Ameas[100] ); - - computeN1mean( data_roi ); // prepare spike removal - removeSpikes( data_roi ); // output in Vcorr + + computeN1mean( RegionOfInterest ); // prepare spike removal + removeSpikes (RegionOfInterest ); // output in Vcorr // filter Vcorr with sliding average using FIR filter function - factfir(b_slide , a_slide, k_slide, Vcorr, Vslide); - - for ( unsigned int sl = 0; sl < data_roi; sl++){ + // 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 Fill( Vslide[sl] ) ; + hBaseline[pix-firstpixel]->Fill( Vslide[sl] ) ; } } @@ -204,21 +176,17 @@ FILE *fp; - TString fName; // name of the histogram file - /* create the filename for the histogram file */ - fName = datafilename; // use the name of the tree file - fName.Remove(fName.Length() - 5); // remove the extension .root - fName += "_bsl.txt"; // add the new extension + 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, "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 = 0; pix < NPIX; pix++) { - //printf("Maximum bin pix[%d] %f\n", pix ,hBaseline[pix]->GetMaximumBin() ); + for (int pix = firstpixel; pix < firstpixel+npixel; pix++) { float vmaxprob = hBaseline[pix]->GetXaxis()->GetBinCenter( - hBaseline[pix]->GetMaximumBin() ); + hBaseline[pix-firstpixel]->GetMaximumBin() ); fprintf( fp, "%8.3f", vmaxprob ); @@ -227,24 +195,23 @@ hpltMeanBsl->SetBinContent(pix+1, vmaxprob ); - hRmsBsl->Fill(hBaseline[pix]->GetRMS() ); + hRmsBsl->Fill(hBaseline[pix-firstpixel]->GetRMS() ); hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() ); } + fprintf( fp, "\n" ); fclose( fp ); - SaveHistograms( datafilename ); - - TCanvas * cMeanBsl = new TCanvas(); - cMeanBsl->cd(); - hMeanBsl->Draw(); - //canv_mean->Modified(); - cMeanBsl->Update(); - - TCanvas * cRmsBsl = new TCanvas(); - cRmsBsl->cd(); - hRmsBsl->Draw(); - //canv_rms->Modified(); - cMeanBsl->Update(); - + 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 ); } @@ -253,5 +220,5 @@ const float fract = 0.8; - float x, xp, xpp, x3p; + float x, xp, xpp; // assume that there are no spikes @@ -323,10 +290,10 @@ // printf("pixel = %d, slice = %d\n", slice, pixel); - pixel_pt = pixel * data_roi; + pixel_pt = pixel * RegionOfInterest; slice_pt = pixel_pt + slice; - drs_cal_offset = ( slice + data_offset[ pixel ] )%data_roi; + drs_cal_offset = ( slice + StartCells[ pixel ] )%RegionOfInterest; cal_pt = pixel_pt + drs_cal_offset; - vraw = data[ slice_pt ] * dconv; + vraw = Data[ slice_pt ] * dconv; vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35; @@ -334,5 +301,5 @@ } -void BookHistos( int Samples ){ +void BookHistos( int Samples , int NumberOfPixel){ // booking and parameter settings for all histos @@ -345,5 +312,5 @@ printf("inside BookHistos\n"); - for( int i = 0; i < NPIX; i++ ) { + for( int i = 0; i < NumberOfPixel; i++ ) { // printf("call sprintf [%d]\n", i ); @@ -388,5 +355,5 @@ -void SaveHistograms( char * loc_fname ){ +void SaveHistograms( const char * loc_fname ){ TString fName; // name of the histogram file @@ -394,6 +361,7 @@ /* 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.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) @@ -407,30 +375,59 @@ } // end of function: void ana::SaveHistograms( void ) -int FOpenDataFile(fits &datafile){ - -/* cout << "-------------------- Data Header -------------------" << endl; - datafile.PrintKeys(); - cout << "------------------- Data Columns -------------------" << endl; - datafile.PrintColumns(); - */ - - //Get the size of the data column - data_roi = datafile.GetUInt("NROI"); // Value from header - data_px = datafile.GetUInt("NPIX"); - data_n = datafile.GetN("Data"); //Size of column "Data" = #Pixel x ROI - - //Set the sizes of the data vectors - data.resize(data_n,0); - data_offset.resize(data_px,0); - - //Link the data to variables - datafile.SetRefAddress("EventNum", data_num); - datafile.SetVecAddress("Data", data); - datafile.SetVecAddress("StartCellData", data_offset); - datafile.GetRow(0); - - cout << "Opening data file successful..." << endl; - - return datafile.GetNumRows() ; +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 &Data, // vector, which will be filled with raw data + vector &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; }