Changeset 12356

Timestamp:

11/01/11 14:59:46 (13 years ago)

Author:

neise

Message:

adjusted fbsl for ISDC cluster .. I basically renamed a few variables and generally tidied up a bit ... not sure what actually helped... maybe its still not solved.

Location:

fact/tools/rootmacros

Files:

• README.txt (modified) (1 diff)
• factfir.C (modified) (1 diff)
• fbsl.C (modified) (17 diffs)

fact/tools/rootmacros/README.txt

@@ -1 +1 @@
 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.

fact/tools/rootmacros/factfir.C

@@ -2 +2 @@
 #include <deque>
 
-
-#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<double> &a, int k, vector<float> &source, vector<float> &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;

fact/tools/rootmacros/fbsl.C

@@ -16 +16 @@
 #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 +27 @@
 #define FAD_MAX_SAMPLES 1024
 
-vector<int16_t> data;
-vector<int16_t> data_offset;
-
-unsigned int data_num;
-
-size_t data_n;
-UInt_t data_px;
-UInt_t data_roi;
+//vector<int16_t> data;
+//vector<int16_t> data_offset;
+//unsigned int data_num;
+//size_t data_n;
+//UInt_t data_px;
+//UInt_t data_roi;
+        vector<int16_t> Data;                           // vector, which will be filled with raw data
+        vector<int16_t> 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 +48 @@
 vector<float> 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<int16_t> &Data,                  // vector, which will be filled with raw data
+        vector<int16_t> &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 +72 @@
 vector<float> Vcfd2(FAD_MAX_SAMPLES);    // CDF result + 2nd sliding average
 
-#include "factfir.C"
 
 float getValue( int, int );
@@ -75 +89 @@
 
 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 +97 @@
 TObjArray hListBaseline;
 
-void BookHistos( int );
-void SaveHistograms( char * );
+void BookHistos( int , int);
+void SaveHistograms( const char * );
 
 // Create a canvas
@@ -97 +109 @@
 
 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<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;
-        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 +135 @@
     // 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 <RegionOfInterest ; sl++){
                // hPixelCellData.Fill( sl, Vcorr[sl] );
-               hBaseline[pix]->Fill( Vslide[sl] ) ;
+               hBaseline[pix-firstpixel]->Fill( Vslide[sl] ) ;
             }   
             }
@@ -204 +176 @@
 
     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 +195 @@
         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 +220 @@
 
     const float fract = 0.8;
-    float x, xp, xpp, x3p;
+    float x, xp, xpp;
 
     // assume that there are no spikes
@@ -323 +290 @@
     // 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 +301 @@
 }
 
-void BookHistos( int Samples ){
+void BookHistos( int Samples , int NumberOfPixel){
 // booking and parameter settings for all histos
 
@@ -345 +312 @@
     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 +355 @@
 
 
-void SaveHistograms( char * loc_fname ){
+void SaveHistograms( const char * loc_fname ){
 
   TString fName; // name of the histogram file
@@ -394 +361 @@
   /* 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 +375 @@
 } // 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<int16_t> &Data,                  // vector, which will be filled with raw data
+        vector<int16_t> &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;
 }