Index: fact/tools/rootmacros/fpeak_cdf.C =================================================================== --- fact/tools/rootmacros/fpeak_cdf.C (revision 12368) +++ fact/tools/rootmacros/fpeak_cdf.C (revision 12369) @@ -31,16 +31,19 @@ #include "factfir.C" - +#include "FOpenDataFile.h" +#include "FOpenDataFile.c" + +#include "DrsCalibration.C" +#include "DrsCalibration.h" + +bool breakout=false; + +int NEvents; vector AllPixelDataVector; vector 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,7 +52,4 @@ vector drs_triggeroffsetmean; -int FOpenDataFile( fits & ); - - vector Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude vector N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors @@ -62,5 +62,5 @@ -float getValue( int, int ); +//float getValue( int, int ); void computeN1mean( int ); void removeSpikes( int ); @@ -79,6 +79,6 @@ 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,18 +92,28 @@ 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,16 +140,5 @@ } -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 a_slide(k_slide, 1); -// double b_slide = k_slide; + // CFD filter settings @@ -150,22 +149,29 @@ a_cfd[k_cfd-1]=1.; -// 2nd slinding window filter -// int ks2 = 16; -// vector 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,24 +180,5 @@ 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,19 +186,5 @@ 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,18 +194,8 @@ } - // 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,7 +223,8 @@ 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; isize(); i++){ @@ -280,4 +244,5 @@ 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,12 +309,4 @@ 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,6 +330,6 @@ 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,5 +342,5 @@ - SaveHistograms( datafilename ); + SaveHistograms( OutRootFileName ); return( 0 ); @@ -401,5 +358,4 @@ 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,12 +367,8 @@ 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,6 +377,4 @@ 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,26 +408,4 @@ } // 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,11 +415,16 @@ 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,9 +463,4 @@ 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,16 +499,7 @@ } -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,33 +509,7 @@ 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 )