Changeset 12369 for fact/tools/rootmacros

Timestamp:

11/04/11 12:45:30 (13 years ago)

Author:

neise

Message:

cleaned up a bit ...still needs testing

File:

• fact/tools/rootmacros/fpeak_cdf.C (modified) (23 diffs)

fact/tools/rootmacros/fpeak_cdf.C

@@ -31 +31 @@
 #include "factfir.C"
 
-
+#include "FOpenDataFile.h"
+#include "FOpenDataFile.c"
+
+#include "DrsCalibration.C"
+#include "DrsCalibration.h"
+
+bool breakout=false;
+
+int NEvents;
 vector<int16_t> AllPixelDataVector;
 vector<int16_t> StartCellVector;
-
 unsigned int CurrentEventID;
-
-bool breakout=false;
-
-size_t ROIxNOP;
+size_t PXLxROI;
+UInt_t RegionOfInterest;
 UInt_t NumberOfPixels;
-UInt_t RegionOfInterest;
-int NEvents;
 
 size_t drs_n;
@@ -49 +52 @@
 vector<float> drs_triggeroffsetmean;
 
-int FOpenDataFile( fits & );
-
-
 vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
 vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
@@ -62 +62 @@
 
 
-float getValue( int, int );
+//float getValue( int, int );
 void computeN1mean( int );
 void removeSpikes( int );
@@ -79 +79 @@
 TH1F* h;
 TH1F *debugHistos;
-TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
-                    // x = pixel id, y = DRS cell id
+TH2F* hStartCell;               // id of the DRS physical pipeline cell where readout starts
+                                                // x = pixel id, y = DRS cell id
 
 TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
@@ -92 +92 @@
 
 void BookHistos( );
-void SaveHistograms( char * );
-
-int searchSinglesPeaks = 0;
-
+void SaveHistograms(const char * );
+
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// 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 fpeak(
         char *datafilename              = "../../20111011_055.fits.gz",
         const char *drsfilename = "../../20111011_054.drs.fits.gz",
+        const char *OutRootFileName = "./fpeak_cdf.Coutput.root",
+        int firstevent                  = 0,
         int nevents                     = -1,
-        int firstevent                  = 0,
-        int PixelID                             = -1,
+        int firstpixel                  = 0,
+        int npixel                              = -1,
         bool spikeDebug = false,
         int verbosityLevel = 1 // different verbosity levels can be implemented here
         )
 {
+gStyle->SetPalette(1,0);
+gROOT->SetStyle("Plain");
 
 // Create (pointer to) Canvases, which are used in every run,
@@ -130 +140 @@
         }
 
-gStyle->SetPalette(1,0);
-gROOT->SetStyle("Plain");
-// read FACT raw data 
-//      * remove spikes
-//      * calculate baseline
-//      * find peaks (CFD and normal discriminator)
-//      * compute pulse height and pulse integral spektrum of the peaks
-
-// sliding window filter settings
-//      int k_slide = 16;
-//      vector<double> a_slide(k_slide, 1);
-//      double b_slide = k_slide;
+
 
 // CFD filter settings
@@ -150 +149 @@
         a_cfd[k_cfd-1]=1.;
 
-// 2nd slinding window filter
-//      int ks2 = 16;
-//      vector<double> as2(ks2, 1);
-//      double bs2 = ks2;
-        
-// Open the data file
-        fits *datafile = new fits( datafilename );
-        if (!datafile) {
-            printf( "Could not open the file: %s\n", datafilename );
-                  return 1;
+        fits * datafile;
+        NEvents = OpenDataFile( datafilename, &datafile,
+                AllPixelDataVector, StartCellVector, CurrentEventID,
+                RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
+        if (NEvents == 0){
+                cout << "return code of FOpenDataFile:" << datafilename<< endl;
+                cout << "is zero -> aborting." << endl;
+                return 1;
         }
-        
-// access data
-        NEvents = FOpenDataFile( *datafile );
-        printf("number of events in file: %d\n", NEvents);
-  if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
-
-//Get the DRS calibration
+
+        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 <<"Totel # of Pixel: "<< NumberOfPixels << "\t";
+        if ( npixel == -1 || npixel > (int)NumberOfPixels  ) npixel = (int)NumberOfPixels; // -1 means all!
+        if (verbosityLevel > 0)
+                cout <<"of, which "<< nevents << "will be processed"<< endl;
+
+
+
+        //Get the DRS calibration
         FOpenCalibFile( drsfilename,
                                         drs_basemean,
@@ -174 +180 @@
                                         drs_n);
 
-//Check the sizes of the data columns
-        if (drs_n != 1474560) {
-                cerr << "error: DRS calib file has wrong ...erm...size ... drs_n is: "
-                        << drs_n << endl;
-                cerr << " Aborting." << endl;
-                return 1;
-        }
-
-        if(ROIxNOP != 1474560)
-        {
-                cout << "warning: data_n should better be 1440x1024=1474560, but is "
-                        << ROIxNOP << endl;
-                cout << "this script is not guaranteed to run under these "
-                        <<" circumstances....any way ... it is never guaranteed." << endl;
-        }
-// Book the histograms
-
         BookHistos( );
-
-        float calibratedVoltage;
 
         for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
@@ -199 +186 @@
                 datafile->GetRow( ev );
 
-                for ( int pix = 0; pix < 1440; pix++ ){
-
-                        hStartCell->Fill( pix, StartCellVector[pix] );
-
-                        // this is a stupid hack ... there is more code at the
-                        // end of this loop to complete this hack ...
-                        // beginning with if (PixelID != -1) as well.
-                        if (PixelID != -1) {
-                                pix = PixelID;
-                                if (verbosityLevel > 0){
-                                        cout << "Processing Event number: " << CurrentEventID << "\t"
-                                                << "Pixel number: "<< pix << endl;
-                                }
-                        }
-
+                for ( int pix = firstpixel; pix < firstpixel+npixel; pix++ ){
                         if (verbosityLevel > 0){
                                 if (pix % 20 ==0){
@@ -221 +194 @@
                         }
 
-                        // compute the DRs calibrated values and put them into Ameas[]
-                        for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
-                                calibratedVoltage = getValue( sl, pix);
-                                if (verbosityLevel > 10){
-                                        printf("calibratedVoltage = %f\n", calibratedVoltage);
-                                }
-                                Ameas[ sl ] = calibratedVoltage;
-
-                                // in case we want to plot this ... we need to put it into a
-                                // Histgramm
-                                if (spikeDebug){
-                                        debugHistos[Ameas_].SetBinContent(sl, calibratedVoltage);
-                                }
-                        }
+                        applyDrsCalibration( Ameas,pix,
+                                drs_basemean, drs_gainmean, drs_triggeroffsetmean,
+                                RegionOfInterest, AllPixelDataVector, StartCellVector);
+
                         // operates on Ameas[] and writes to N1mean[]
                         computeN1mean( RegionOfInterest );
@@ -260 +223 @@
                         EnlargeRegion(*zXings, 10, 10);
                         findAbsMaxInRegions(*zXings, Vslide);
-                        removeMaximaBelow( *zXings, 11.0, 0);
-                        removeMaximaAbove( *zXings, 14.0, 0);
-
+                        //removeMaximaBelow( *zXings, 11.0, 0);
+                        //removeMaximaAbove( *zXings, 14.0, 0);
+
+                        // fill maxima in Histogram
                         if (zXings->size() != 0 ){
                                 for (unsigned int i=0; i<zXings->size(); i++){
@@ -280 +244 @@
                         if ( spikeDebug ){
                                 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
+                                        debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
                                         debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
                                         debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
@@ -344 +309 @@
 
                         delete zXings;
-
-                        // this is the 2nd part of the ugly hack in the beginning.
-                        // this makes sure, that the loop ends
-                        if (PixelID != -1){
-                                pix = 1440;
-                        }
-
-
                         if (breakout)
                                 break;
@@ -373 +330 @@
 
                 cTemplate->cd();
-                hTemp_Array[PixelID]->GetYaxis()->SetRangeUser(-10,40);
-                hTemp_Array[PixelID]->Draw("COLZ");
+                hTemp_Array[firstpixel]->GetYaxis()->SetRangeUser(-10,40);
+                hTemp_Array[firstpixel]->Draw("COLZ");
                 cTemplate->Modified();
                 cTemplate->Update();
@@ -385 +342 @@
 
 
-        SaveHistograms( datafilename );
+        SaveHistograms( OutRootFileName );
 
         return( 0 );
@@ -401 +358 @@
     for ( int i = 0; i < Samples; i++) {
 
-    // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ i ] );
 
     x = Ameas[i] - N1mean[i];
@@ -411 +367 @@
             x3p = Ameas[i+3] - N1mean[i+3];
 
-            // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p);
 
             if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
-                // printf("double spike candidate\n");
                 Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
                 Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
-                // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]);
-                // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]);
                 i = i + 3; 
             }
@@ -425 +377 @@
                 if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
                     Vcorr[i+1] = N1mean[i+1];
-                    // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
-                    // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
                     N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
                     i = i + 2;//do not care about the next sample it was the spike
@@ -458 +408 @@
 } // end of computeN1mean computation
 
-float getValue( int slice, int pixel ){
-        const float dconv = 2000/4096.0;
-
-        float vraw, vcal;
-
-        unsigned int pixel_pt;
-        unsigned int slice_pt;
-        unsigned int cal_pt;
-        unsigned int drs_cal_offset;
-
-        // printf("pixel = %d, slice = %d\n", slice, pixel);
-
-        pixel_pt = pixel * RegionOfInterest;
-        slice_pt = pixel_pt + slice;
-        drs_cal_offset = ( slice + StartCellVector[ pixel ] )%RegionOfInterest;
-        cal_pt    = pixel_pt + drs_cal_offset;
-
-        vraw = AllPixelDataVector[ slice_pt ] * dconv;
-        vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
-
-        return( vcal );
-}
 
 // booking and parameter settings for all histos
@@ -487 +415 @@
         char hTitle[500];
 
+        TString name,title;
+        title = "all events all slices of pixel ";
+        name = "base";
         TH1F *h;
 
         for( int i = 0; i < NPIX; i++ ) {
-                sprintf(&hTitle[0],"all events all slices of pixel %d", i);
-                sprintf(&hName[0],"base%d", i);
-
-                h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );
+                title += i;
+                name += i;
+//              sprintf(&hTitle[0],"all events all slices of pixel %d", i);
+//              sprintf(&hName[0],"base%d", i);
+
+                h = new TH1F( name, title, 400, -99.5 ,100.5 );
 
                 h->GetXaxis()->SetTitle( "Sample value (mV)" );
@@ -530 +463 @@
         hAmplSpek_discri->GetXaxis()->SetTitle( "amplitude in mV" );
         hList.Add( hAmplSpek_discri );
-
-        hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
-        hStartCell->GetXaxis()->SetTitle( "pixel" );
-        hStartCell->GetXaxis()->SetTitle( "slice" );
-        hList.Add( hStartCell );
 
         debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
@@ -571 +499 @@
 }
 
-void SaveHistograms( char * loc_fname ){
-
-        TString fName; // name of the histogram file
-
-        // create the filename for the histogram file
-        fName = loc_fname; // use the name of the tree file
-        // TODO ... next statement doesn't work for ".fits.gz"
-        fName.Remove(fName.Length() - 5); // remove the extension .fits
-        fName += "_discri.root"; // add the new extension
-
+void SaveHistograms( const char * loc_fname ){
         // create the histogram file (replace if already existing)
-        TFile tf( fName, "RECREATE");
+        TFile tf( loc_fname, "RECREATE");
 
         hList.Write(); // write the major histograms into the top level directory
@@ -590 +509 @@
         tf.cd("..");
         tf.mkdir("PulseTempplates");
-        tf.cd("PulseTempplates");
+        tf.cd("PulseTemplates");
         hListTemplates.Write(); // write histos into subdirectory
 
         tf.Close(); // close the file
-} // end of SaveHistograms( char * loc_fname )
-
-int FOpenDataFile(fits &datafile){
-
-//      cout << "-------------------- Data Header -------------------" << endl;
-//      datafile.PrintKeys();
-//      cout << "------------------- Data Columns -------------------" << endl;
-//      datafile.PrintColumns();
-
-        //Get the size of the data column
-        RegionOfInterest        = datafile.GetUInt("NROI");
-        NumberOfPixels          = datafile.GetUInt("NPIX");
-        //Size of column "Data" = #Pixel x ROI
-        ROIxNOP                         = datafile.GetN("Data");
-
-        //Set the sizes of the data vectors
-        AllPixelDataVector.resize(ROIxNOP,0);
-        StartCellVector.resize(NumberOfPixels,0);
-
-        //Link the data to variables
-        datafile.SetRefAddress("EventNum",  CurrentEventID);
-        datafile.SetVecAddress("Data", AllPixelDataVector);
-        datafile.SetVecAddress("StartCellData", StartCellVector);
-        datafile.GetRow(0);
-
-        return datafile.GetNumRows() ;
-}
+} // end of SaveHistograms(const char * loc_fname )