Ignore:
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.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • fact/tools/rootmacros/fbsl.C

    r12307 r12356  
    1616#define HAVE_ZLIB
    1717#include "fits.h"
    18 //#include "TPKplotevent.c"
    19 //#include "FOpenDataFile.c"
    2018#include "FOpenCalibFile.c"
    2119
    2220#include "zerosearch.C"
     21#include "factfir.C"
    2322
    2423#define NPIX  1440
     
    2827#define FAD_MAX_SAMPLES 1024
    2928
    30 vector<int16_t> data;
    31 vector<int16_t> data_offset;
    32 
    33 unsigned int data_num;
    34 
    35 size_t data_n;
    36 UInt_t data_px;
    37 UInt_t data_roi;
     29//vector<int16_t> data;
     30//vector<int16_t> data_offset;
     31//unsigned int data_num;
     32//size_t data_n;
     33//UInt_t data_px;
     34//UInt_t data_roi;
     35        vector<int16_t> Data;                           // vector, which will be filled with raw data
     36        vector<int16_t> StartCells;     // vector, which will be filled with DRS start positions
     37        unsigned int EventID;                           // index of the current event
     38        UInt_t RegionOfInterest;                // Width of the Region, read out of the DRS
     39        UInt_t NumberOfPixels;                  // Total number of pixel, read out of the camera
     40        size_t PXLxROI;                                                 // Size of column "Data" = #Pixel x ROI
     41
    3842int NEvents;
    3943int NBSLeve = 1000;
     
    4448vector<float> drs_triggeroffsetmean;
    4549
    46 int FOpenDataFile( fits & );
     50size_t FOpenDataFile(
     51        const char *datafilename,               // path to fits file containing FACT raw data
     52        fits * * datafile,                              // pointer to pointer, where to return the fits object
     53        vector<int16_t> &Data,                  // vector, which will be filled with raw data
     54        vector<int16_t> &StartCells,    // vector, which will be filled with DRS start positions
     55        unsigned int &EventID,                  // index of the current event
     56        UInt_t &RegionOfInterest,               // Width of the Region, read out of the DRS
     57        UInt_t &NumberOfPixels,                 // Total number of pixel, read out of the camera
     58        size_t &PXLxROI,                                // Size of column "Data" = #Pixel x ROI
     59                        // this can be used, to x-check RegionOfInterest and NumberOfPixels
     60        int VerbosityLevel=1
     61);
    4762
    4863
     
    5772vector<float> Vcfd2(FAD_MAX_SAMPLES);    // CDF result + 2nd sliding average
    5873
    59 #include "factfir.C"
    6074
    6175float getValue( int, int );
     
    7589
    7690TH1F* h;
    77 TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
    78                     // x = pixel id, y = DRS cell id
    7991TH2F hPixelCellData(
    8092    "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
     
    8597TObjArray hListBaseline;
    8698
    87 void BookHistos( int );
    88 void SaveHistograms( char * );
     99void BookHistos( int , int);
     100void SaveHistograms( const char * );
    89101
    90102// Create a canvas
     
    97109
    98110int fbsl(
    99         char *datafilename              = "../../20110804_024.fits.gz",
    100         const char *drsfilename = "../../20110804_023.drs.fits.gz",
    101         int pixelnr                     = 0,
     111        const char *datafilename                = "path-to-datafile.fits.gz",
     112        const char *drsfilename = "path-to-calibfile.drs.fits.gz",
     113        const char *TextOutFileName = "./appendfile.txt",
     114        const char *RootOutFileName = "./datafile.root",
    102115        int firstevent                  = 0,
    103         int nevents                     = -1 ){
    104 // read and analyze FACT raw data
    105 
    106 // sliding window filter settings
    107         int k_slide = 16;
    108         vector<double> a_slide(k_slide, 1);
    109         double b_slide = k_slide;
    110 
    111 // CFD filter settings
    112         int k_cfd = 10;
    113         vector<double> a_cfd(k_cfd, 0);
    114         double b_cfd = 1.;
    115         a_cfd[0]=0.75;
    116         a_cfd[k_cfd-1]=-1.;
    117 
    118 // 2nd slinding window filter
    119         int ks2 = 16;
    120         vector<double> as2(ks2, 1);
    121         double bs2 = ks2;
    122         gROOT->SetStyle("Plain");
    123        
    124 //-------------------------------------------
    125 // Open the file
    126 //-------------------------------------------
    127         fits datafile( datafilename );
    128         if (!datafile) {
    129             printf( "Could not open the file: %s\n", datafilename );
    130             return 1;
    131         }
    132        
    133     // access data
    134         NEvents = FOpenDataFile( datafile );
     116        int nevents                     = -1,
     117        int firstpixel          = 0,
     118        int npixel                              = -1,
     119        bool produceGraphic = false
     120){
     121        fits * datafile = NULL;
     122        NEvents =  FOpenDataFile(
     123                        datafilename,
     124                        &datafile,
     125                        Data,
     126                        StartCells,
     127                        EventID,
     128                        RegionOfInterest,
     129                        NumberOfPixels,
     130                        PXLxROI);
    135131
    136132    printf("number of events in file: %d\n", NEvents);
     
    139135    // the user would like to read. nevents = -1 means: read all
    140136    if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;
    141 
    142 
    143     //  FILE *pedFile;
    144         // pedFile = fopen("t.txt","u");
    145     // fprintf(pedFile,"# Pedestal Data");
    146     // fclose( pedFile );
    147 //-------------------------------------------
    148 //Get the DRS calibration
    149 //-------------------------------------------
     137       
     138                if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels;
    150139
    151140        FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
    152        
    153 
    154     //Check the sizes of the data columns
    155         if(drs_n!=data_n)
    156         {
    157                 cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
    158                 return 1;
    159         }
    160     // Book the histograms
    161     // printf("call BookHistos\n");
    162     BookHistos( data_roi );
    163     // printf("returned from BookHistos\n");
    164     // Loop over events
    165         cout << "--------------------- Data --------------------" << endl;
    166 
    167     float value;
     141
     142  BookHistos( RegionOfInterest, npixel );
    168143   
    169     // loop over events
    170 //    for ( int ev = firstevent; ev < nevents; ev++) {
     144                // loop over events
    171145    for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
    172146
    173             datafile.GetRow( ev );
     147            datafile->GetRow( ev );
    174148               
    175149                if ( ev % 100 == 0 ){
    176                         cout << "Event number: " << data_num << endl;
     150                        cout << "Event ID: " << EventID << endl;
    177151                }
    178152               
    179153        // loop over pixel
    180         for ( int pix = 0; pix < 1440; pix++ ){
    181 
    182             // hStartCell->Fill( pix, data_offset[pix] );
    183        
     154        for ( int pix = firstpixel ; pix < npixel+firstpixel ; pix++ ){
     155
    184156            // loop over DRS slices
    185             for ( unsigned int sl = 0; sl < data_roi; sl++){
    186 
     157            for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
    187158                Ameas[ sl ] = getValue(sl, pix);
    188                
    189159            }
    190             // printf("Ameas[100] = %f\n", Ameas[100] );
    191 
    192             computeN1mean( data_roi ); // prepare spike removal
    193             removeSpikes( data_roi );  // output in Vcorr
     160
     161            computeN1mean( RegionOfInterest ); // prepare spike removal
     162            removeSpikes (RegionOfInterest );  // output in Vcorr
    194163
    195164            // filter Vcorr with sliding average using FIR filter function
    196             factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
    197 
    198             for ( unsigned int sl = 0; sl < data_roi; sl++){
     165                                                // 8 is here the HalfWidth of the sliding average window. 
     166                                                sliding_avg(Vcorr, Vslide, 8);
     167
     168                                                //            factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
     169
     170            for ( unsigned int sl = 0; sl <RegionOfInterest ; sl++){
    199171               // hPixelCellData.Fill( sl, Vcorr[sl] );
    200                hBaseline[pix]->Fill( Vslide[sl] ) ;
     172               hBaseline[pix-firstpixel]->Fill( Vslide[sl] ) ;
    201173            }   
    202174            }
     
    204176
    205177    FILE *fp;
    206     TString fName; // name of the histogram file
    207     /* create the filename for the histogram file */
    208     fName = datafilename; // use the name of the tree file
    209     fName.Remove(fName.Length() - 5); // remove the extension .root
    210     fName += "_bsl.txt"; // add the new extension
     178    TString fName;
     179    fName = TextOutFileName;
    211180
    212181    fp = fopen(fName, "a+");
    213182    fprintf( fp, "FILE: %s\n", datafilename );
    214     fprintf( fp, "NEVENTS: %d\n", NEvents);
    215     NBSLeve = NEvents; // this has to be changed when computation on a subset of a run is implemented
     183    fprintf( fp, "NEVENTS: %d\n", nevents);
     184    NBSLeve = nevents; // this has to be changed when computation on a subset of a run is implemented
    216185    fprintf( fp, "NBSLEVE:  %d\n", NBSLeve );
    217186
    218     for (int pix = 0; pix < NPIX; pix++) {
    219         //printf("Maximum bin pix[%d] %f\n", pix ,hBaseline[pix]->GetMaximumBin() );
     187    for (int pix = firstpixel; pix < firstpixel+npixel; pix++) {
    220188
    221189        float vmaxprob = hBaseline[pix]->GetXaxis()->GetBinCenter(
    222             hBaseline[pix]->GetMaximumBin() );
     190            hBaseline[pix-firstpixel]->GetMaximumBin() );
    223191
    224192        fprintf( fp, "%8.3f", vmaxprob );
     
    227195        hpltMeanBsl->SetBinContent(pix+1, vmaxprob );
    228196
    229         hRmsBsl->Fill(hBaseline[pix]->GetRMS() );
     197        hRmsBsl->Fill(hBaseline[pix-firstpixel]->GetRMS() );
    230198        hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() );   
    231199    }
     200    fprintf( fp, "\n" );
    232201
    233202    fclose( fp );
    234203
    235     SaveHistograms( datafilename );
    236 
    237     TCanvas * cMeanBsl = new TCanvas();
    238     cMeanBsl->cd();
    239     hMeanBsl->Draw();
    240     //canv_mean->Modified();
    241     cMeanBsl->Update();
    242 
    243     TCanvas * cRmsBsl = new TCanvas();
    244     cRmsBsl->cd();
    245     hRmsBsl->Draw();
    246     //canv_rms->Modified();
    247     cMeanBsl->Update();
    248 
     204    SaveHistograms( RootOutFileName );
     205                if (produceGraphic){
     206        TCanvas * cMeanBsl = new TCanvas();
     207        cMeanBsl->cd();
     208        hMeanBsl->Draw();
     209        cMeanBsl->Update();
     210
     211        TCanvas * cRmsBsl = new TCanvas();
     212        cRmsBsl->cd();
     213        hRmsBsl->Draw();
     214        cMeanBsl->Update();
     215                }
    249216        return( 0 );
    250217}
     
    253220
    254221    const float fract = 0.8;
    255     float x, xp, xpp, x3p;
     222    float x, xp, xpp;
    256223
    257224    // assume that there are no spikes
     
    323290    // printf("pixel = %d, slice = %d\n", slice, pixel);
    324291
    325     pixel_pt = pixel * data_roi;
     292    pixel_pt = pixel * RegionOfInterest;
    326293    slice_pt = pixel_pt + slice;
    327     drs_cal_offset = ( slice + data_offset[ pixel ] )%data_roi;
     294    drs_cal_offset = ( slice + StartCells[ pixel ] )%RegionOfInterest;
    328295    cal_pt    = pixel_pt + drs_cal_offset;
    329296
    330     vraw = data[ slice_pt ] * dconv;
     297    vraw = Data[ slice_pt ] * dconv;
    331298    vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
    332299
     
    334301}
    335302
    336 void BookHistos( int Samples ){
     303void BookHistos( int Samples , int NumberOfPixel){
    337304// booking and parameter settings for all histos
    338305
     
    345312    printf("inside BookHistos\n");
    346313
    347     for( int i = 0; i < NPIX; i++ ) {
     314    for( int i = 0; i < NumberOfPixel; i++ ) {
    348315
    349316        // printf("call sprintf [%d]\n", i );
     
    388355
    389356
    390 void SaveHistograms( char * loc_fname ){
     357void SaveHistograms( const char * loc_fname ){
    391358
    392359  TString fName; // name of the histogram file
     
    394361  /* create the filename for the histogram file */
    395362  fName = loc_fname; // use the name of the tree file
    396   fName.Remove(fName.Length() - 5); // remove the extension .root
    397   fName += "_histo.root"; // add the new extension
     363  //fName.Remove(fName.Length() - 5); // remove the extension .root
     364  //fName += "_histo.root"; // add the new extension
     365  //fName += ".root";
    398366
    399367  TFile tf( fName, "RECREATE"); // create the histogram file (replace if already existing)
     
    407375} // end of function: void ana::SaveHistograms( void )
    408376
    409 int FOpenDataFile(fits &datafile){
    410 
    411 /*  cout << "-------------------- Data Header -------------------" << endl;
    412     datafile.PrintKeys();
    413     cout << "------------------- Data Columns -------------------" << endl;
    414     datafile.PrintColumns();
    415     */
    416 
    417     //Get the size of the data column
    418     data_roi = datafile.GetUInt("NROI");  // Value from header
    419     data_px = datafile.GetUInt("NPIX");
    420     data_n = datafile.GetN("Data");         //Size of column "Data" = #Pixel x ROI
    421    
    422     //Set the sizes of the data vectors
    423     data.resize(data_n,0);
    424     data_offset.resize(data_px,0);
    425    
    426     //Link the data to variables
    427     datafile.SetRefAddress("EventNum", data_num);
    428     datafile.SetVecAddress("Data", data);
    429     datafile.SetVecAddress("StartCellData", data_offset);
    430     datafile.GetRow(0);
    431    
    432     cout << "Opening data file successful..." << endl;
    433 
    434     return datafile.GetNumRows() ;
     377size_t FOpenDataFile(
     378        const char *datafilename,               // path to fits file containing FACT raw data
     379        fits * * datafile,                              // ptr to pointer, where to return the fits object
     380        vector<int16_t> &Data,                  // vector, which will be filled with raw data
     381        vector<int16_t> &StartCells,    // vector, which will be filled with DRS start positions
     382        unsigned int &EventID,                  // index of the current event
     383        UInt_t &RegionOfInterest,               // Width of the Region, read out of the DRS
     384        UInt_t &NumberOfPixels,                 // Total number of pixel, read out of the camera
     385        size_t &PXLxROI,                                // Size of column "Data" = #Pixel x ROI
     386                        // this can be used, to x-check RegionOfInterest and NumberOfPixels
     387        int VerbosityLevel                              //
     388) {
     389        size_t NumberOfEvents;
     390        *datafile = new fits(datafilename);
     391        if (!(*(*datafile))) {
     392                if (VerbosityLevel > 0)
     393                        cout << "Couldn't properly open the datafile: " << datafilename << endl;
     394                return 0;
     395        }
     396
     397        NumberOfEvents = (*datafile)->GetNumRows();
     398        if (NumberOfEvents < 1){
     399                if (VerbosityLevel > 0){
     400                        cout << "Warning in FOpenDataFile of file: " << datafilename << endl;
     401                        cout << "the file contains no events." << endl;
     402                }
     403        }
     404
     405        RegionOfInterest = (*datafile)->GetUInt("NROI");
     406        NumberOfPixels = (*datafile)->GetUInt("NPIX");
     407        PXLxROI = (*datafile)->GetN("Data");
     408
     409        if ( RegionOfInterest * NumberOfPixels != PXLxROI) // something in the file went wrong
     410        {
     411                if (VerbosityLevel > 0){
     412                        cout << "Warning in FOpenDataFile of file: " << datafilename << endl;
     413                        cout << "RegionOfInterest * NumberOfPixels != PXLxROI" << endl;
     414                        cout << "--> " << RegionOfInterest;
     415                        cout << " * " << NumberOfPixels;
     416                        cout << " = " << RegionOfInterest * NumberOfPixels;
     417                        cout << ", but PXLxROI =" << PXLxROI << endl;
     418                }
     419                return 0;
     420        }
     421
     422        //Set the sizes of the data vectors
     423        Data.resize(PXLxROI, 0);
     424        StartCells.resize(NumberOfPixels, 0);
     425
     426        //Link the data to variables
     427        (*datafile)->SetRefAddress("EventNum", EventID);
     428        (*datafile)->SetVecAddress("Data", Data);
     429        (*datafile)->SetVecAddress("StartCellData", StartCells);
     430
     431        return NumberOfEvents;
    435432}
    436433
Note: See TracChangeset for help on using the changeset viewer.