Index: /fact/tools/rootmacros/README.txt =================================================================== --- /fact/tools/rootmacros/README.txt (revision 12391) +++ /fact/tools/rootmacros/README.txt (revision 12392) @@ -105,2 +105,49 @@ +----------------------------------------------------------------------------------------------- +tpeak.C ROOT macros to produce overlays of data around a found peak in a given +window. + +int tpeak( + char *datafilename = "data/20111016_013.fits.gz", + const char *drsfilename = "../../20111016_011.drs.fits.gz", + const char *OutRootFileName = "../analysis/fpeak_cdf.Coutput.root", + int firstevent = 0, + int nevents = -1, + int firstpixel = 0, + int npixel = -1, + bool spikeDebug = false, + int avg1 = 14, + int avg2 = 8, + int OverlayWindowLeft = 50, + int OverlayWindowRight = 150, + int verbosityLevel = 1, // different verbosity levels can be implemented +here + bool ProduceGraphic = true + ) + +call it like this, if you want to overlay a certain number of peaks for a +single channel +tpeak("data/20111029_017.fits.gz","data/20111029_013.drs.fits.gz","test.root",0,1000,333,1,true,0,0,50,150,1,true) + +set the first bool 'spikeDebug' to false in order to overlay quicker ... every +50th event the canvas is updated. +depending on what kind of peaks you like to detect the sliding average filters +should be set to e.g. 14,8 for singles in a quiet G-APD pedestal run +or 0,0 if you want to see just every thing that might be a signal... + +the window size 50,150 means .. 50 slices to the left of the maximum and 150 +to the right. + +there is a bit of cleaning included in the file, maybe you want to switch it +off, then just search for the calls of these methods... + removeMaximaBelow( *zXings, 3.0); + removeRegionWithMaxOnEdge( *zXings, 2); + removeRegionOnFallingEdge( *zXings, 100); + +and comment them out or play with the settings. + +they are defined in zerosearch.C, which is maybe a bad choice... + + + Index: /fact/tools/rootmacros/tpeak.C =================================================================== --- /fact/tools/rootmacros/tpeak.C (revision 12391) +++ /fact/tools/rootmacros/tpeak.C (revision 12392) @@ -42,20 +42,25 @@ #include "factfir.C" +#include "FOpenDataFile.h" +#include "FOpenDataFile.c" + +#include "DrsCalibration.C" +#include "DrsCalibration.h" + +#include "SpikeRemoval.h" +#include "SpikeRemoval.C" //----------------------------------------------------------------------------- // Decleration of variables //----------------------------------------------------------------------------- - +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; @@ -64,21 +69,9 @@ 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 vector Vcorr(FAD_MAX_SAMPLES); // corrected Values -vector Vdiff(FAD_MAX_SAMPLES); // numerical derivative - vector Vslide(FAD_MAX_SAMPLES); // sliding average result vector Vcfd(FAD_MAX_SAMPLES); // CDF result vector Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average - - -float getValue( int, int ); -float correctDrsOffset( int slice, int pixel ); -void computeN1mean( int ); -void removeSpikes( int ); // histograms @@ -95,14 +88,6 @@ TH1F* h; TH1F *debugHistos; -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 -TH1F *hMeanBsl, *hpltMeanBsl; -TH1F *hRmsBsl, *hpltRmsBsl; -TH2F * hAmplSpek_cfd; -TH2F * hAmplSpek_discri; - TH2F *hPeakOverlay; //histogrammm for overlay of detected Peaks +int hPeakOverlayXaxisLeft,hPeakOverlayXaxisRight; TObjArray hList; @@ -110,22 +95,32 @@ void BookHistos( ); -void SaveHistograms( char * ); - -int searchSinglesPeaks = 0; +void SaveHistograms( const char * ); //----------------------------------------------------------------------------- // Main //----------------------------------------------------------------------------- - int tpeak( - char *datafilename = "../../20111011_055.fits.gz", - const char *drsfilename = "../../20111011_054.drs.fits.gz", - int PixelID = -1, - int firstevent = 0, - int nevents = -1, - bool spikeDebug = false, - int verbosityLevel = 1 // different verbosity levels can be implemented here - ) + char *datafilename = "data/20111016_013.fits.gz", + const char *drsfilename = "../../20111016_011.drs.fits.gz", + const char *OutRootFileName = "../analysis/fpeak_cdf.Coutput.root", + int firstevent = 0, + int nevents = -1, + int firstpixel = 0, + int npixel = -1, + bool spikeDebug = false, + int avg1 = 14, + int avg2 = 8, + int OverlayWindowLeft = 50, + int OverlayWindowRight = 150, + int verbosityLevel = 1, // different verbosity levels can be implemented here + bool ProduceGraphic = true + ) + { +hPeakOverlayXaxisLeft = OverlayWindowLeft; +hPeakOverlayXaxisRight = OverlayWindowRight; + + gStyle->SetPalette(1,0); + gROOT->SetStyle("Plain"); //----------------------------------------------------------------------------- @@ -133,27 +128,15 @@ // also in 'non-debug' runs //----------------------------------------------------------------------------- - - TCanvas * cSpektrum; - TCanvas * cStartCell; - cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400); - cSpektrum->Divide(1,2); - cStartCell = new TCanvas("cStartCell ", "The Startcells of this run", 10,410,400,400); - // Canvases only need if spike Debug, but I want to deklare // the pointers anyway ... - TCanvas *cRawAndSpikeRemoval = NULL; TCanvas *cFiltered = NULL; if (spikeDebug){ - cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",410,10,400,400); - cRawAndSpikeRemoval->Divide(1, 2); - cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400); - cFiltered->Divide(1, 2); + cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",1,10,400,300); + cFiltered->Divide(1, 3); } // Canvases to show the peak template TCanvas *cPeakOverlay = NULL; - cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 410, 10, 600, 600); - -gStyle->SetPalette(1,0); -gROOT->SetStyle("Plain"); + cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 1, 310, 400, 400); + //----------------------------------------------------------------------------- @@ -169,10 +152,4 @@ // Filter-Settings //----------------------------------------------------------------------------- - -// sliding window filter settings - int k_slide = 16; - vector a_slide(k_slide, 1); //erzeugt vector der länge k_slide, und füllt die elemente mit 1 - double b_slide = k_slide; - // CFD filter settings int k_cfd = 10; @@ -182,9 +159,4 @@ a_cfd[k_cfd-1]=1.; -// 2nd slinding window filter - int ks2 = 16; - vector as2(ks2, 1); - double bs2 = ks2; - //----------------------------------------------------------------------------- // prepare datafile @@ -192,14 +164,29 @@ // Open the data file - fits *datafile = new fits( datafilename ); - if (!datafile) { - printf( "Could not open the file: %s\n", datafilename ); - return 1; - } - -// access data - NEvents = FOpenDataFile( *datafile ); - printf("number of events in file: %d\n", NEvents); + + fits * datafile; + // Opens the raw data file and 'binds' the variables given as + // Parameters to the data file. So they are filled with + // raw data as soon as datafile->GetRow(int) is called. + NEvents = OpenDataFile( datafilename, &datafile, + AllPixelDataVector, StartCellVector, CurrentEventID, + RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel); + if (NEvents == 0){ + cout << "return code of OpenDataFile:" << datafilename<< endl; + cout << "is zero -> aborting." << endl; + return 1; + } + + 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 <<"Total # of Pixel: "<< NumberOfPixels << "\t"; + if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all! + if (verbosityLevel > 0) + cout <<"of, which "<< npixel << "will be processed"<< endl; //Get the DRS calibration @@ -210,26 +197,7 @@ 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 myTHR = 3.5; - float calibratedVoltage; //----------------------------------------------------------------------------- // Loops over Every Event and Pixel @@ -243,58 +211,26 @@ //------------------------------------- - 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; - } - } - - if (verbosityLevel > 0){ - if (pix % 20 ==0){ - cout << "Processing Event number: " << CurrentEventID << "\t" - << "Pixel number: "<< pix << endl; - } - } - - // 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); - } - } - // operates on Ameas[] and writes to N1mean[] - computeN1mean( RegionOfInterest ); - - // operates on Ameas[] and N1mean[], then writes to Vcorr[] - removeSpikes( RegionOfInterest ); - - // filter Vcorr with sliding average using FIR filter function - factfir(b_slide , a_slide, k_slide, Vcorr, Vslide); - // filter Vslide with CFD using FIR filter function - factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd); - // filter Vcfd with sliding average using FIR filter function - factfir(bs2 , as2, ks2, Vcfd, Vcfd2); - - - vector *Regions = discriminator( Vslide, myTHR,true , 120); - for (unsigned int p=0; psize(); p++ ){ - hAmplSpek_discri->Fill(pix , Regions->at(p).maxVal); - } + for ( int pix = firstpixel; pix < firstpixel + npixel; pix++ ){ + + if (verbosityLevel > 0){ + if (pix == firstpixel){ + cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl; + } + } + + applyDrsCalibration( Ameas,pix,12,12, + drs_basemean, drs_gainmean, drs_triggeroffsetmean, + RegionOfInterest, AllPixelDataVector, StartCellVector); + + // finds spikes in the raw data, and interpolates the value + // spikes are: 1 or 2 slice wide, positive non physical artifacts + removeSpikes (Ameas, Vcorr); + + // filter Vcorr with sliding average using FIR filter function + sliding_avg(Vcorr, Vslide, avg1); + // filter Vslide with CFD using FIR filter function + factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd); + // filter Vcfd with sliding average using FIR filter function + sliding_avg(Vcfd, Vcfd2, avg2); // peaks in Ameas[] are found by searching for zero crossings @@ -304,12 +240,9 @@ vector * zXings = zerosearch( Vcfd2 , 1 , 1); // zXings means "zero cross ings" - ShiftRegionBy(*zXings, -ks2/2); EnlargeRegion(*zXings, 10, 10); findAbsMaxInRegions(*zXings, Vslide); - - for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){ - // hPixelCellData.Fill( sl, Vcorr[sl] ); - hBaseline[pix]->Fill( Vslide[sl] ) ; - } + removeMaximaBelow( *zXings, 3.0); + removeRegionWithMaxOnEdge( *zXings, 2); + removeRegionOnFallingEdge( *zXings, 100); //following Code should produce the Overlay of peaks @@ -317,11 +250,11 @@ for (it = zXings->begin() ; it < zXings->end() ; it++){ - int a = it->maxPos-100; - int b = it->maxPos+100; - if (a < 0) - a =0; - if (b>1023) - b=1023; - for ( int pos = a; pos <= b; pos++){ + int Left = it->maxPos - OverlayWindowLeft; + int Right = it->maxPos + OverlayWindowRight; + if (Left < 0) + Left =0; + if (Right > (int)Vcorr.size() ) + Right=Vcorr.size() ; + for ( int pos = Left; pos < Right; pos++){ hPeakOverlay->Fill( pos - (it->maxPos), Vcorr[pos]) ; @@ -330,50 +263,31 @@ - - if (zXings->size() != 0 ){ - for (unsigned int i=0; isize(); i++){ - if (verbosityLevel > 1){ - cout << zXings->at(i).maxPos << ":\t" - << zXings->at(i).maxVal <Fill(pix, zXings->at(i).maxVal); - } - } - if ( spikeDebug ){ - for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){ - debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] ); - debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] ); - debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] ); - } - - cRawAndSpikeRemoval->cd( 1); - gPad->SetGrid(); - debugHistos[Ameas_].Draw(); - - cRawAndSpikeRemoval->cd( 2); + + + // TODO do this correct. The vectors should be the rigt ones... this is just luck + debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + + for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){ + debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]); + debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]); + debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] ); + debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] ); + debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] ); + } + + + cFiltered->cd(1); gPad->SetGrid(); debugHistos[Vcorr_].Draw(); - cFiltered->cd(1); + cFiltered->cd(2); gPad->SetGrid(); debugHistos[Vslide_].Draw(); - TLine * thrLine = new TLine(0, myTHR, 1024, myTHR); - thrLine->SetLineColor(kRed); - thrLine->SetLineWidth(2); - thrLine->Draw(); - TLine * OneLine; - vector MyLines; - for (unsigned int p=0; psize(); p++ ){ - OneLine = new TLine( - Regions->at(p).begin , - Regions->at(p).maxVal, - Regions->at(p).end, - Regions->at(p).maxVal); - OneLine->SetLineColor(kRed); - OneLine->SetLineWidth(2); - MyLines.push_back(OneLine); - OneLine->Draw(); - } + TBox *OneBox; vector MyBoxes; @@ -392,5 +306,5 @@ } - cFiltered->cd(2); + cFiltered->cd(3); gPad->SetGrid(); debugHistos[Vcfd2_].Draw(); @@ -399,5 +313,4 @@ zeroline->Draw(); - cRawAndSpikeRemoval->Update(); cFiltered->Update(); @@ -423,13 +336,5 @@ }// end of if(spikeDebug) - delete Regions; 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 = 2000; - } - if (breakout) break; @@ -439,19 +344,6 @@ // end of loop over pixels //------------------------------------- -if (ev % 50 == 0){ - cSpektrum->cd(1); - hAmplSpek_cfd->Draw("COLZ"); - cSpektrum->cd(2); - hAmplSpek_discri->Draw("COLZ"); - cSpektrum->Modified(); - cSpektrum->Update(); - - // updating seems not to work .. - // debug cout - cStartCell->cd(); - hStartCell->Draw(); - cStartCell->Modified(); - cStartCell->Update(); - + if (ProduceGraphic){ + if (ev % 50 == 0){ //OVERLAY PEAKS cPeakOverlay->cd(); @@ -459,8 +351,6 @@ cPeakOverlay->Modified(); cPeakOverlay->Update(); - -} - - + } + } if (breakout) break; @@ -469,271 +359,19 @@ // end of loop over Events //------------------------------------- - - - for (int pix = 0; pix < NumberOfPixels; pix++) { - //printf("Maximum bin pix[%d] %f\n", pix ,hBaseline[pix]->GetMaximumBin() ); - - float vmaxprob = hBaseline[pix]->GetXaxis()->GetBinCenter( - hBaseline[pix]->GetMaximumBin() ); - - printf("%8.3f", vmaxprob ); - - hMeanBsl->Fill( vmaxprob ); - hpltMeanBsl->SetBinContent(pix+1, vmaxprob ); - - hRmsBsl->Fill(hBaseline[pix]->GetRMS() ); - hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() ); - } - - - SaveHistograms( datafilename ); + if (ProduceGraphic){ + //OVERLAY PEAKS + cPeakOverlay->cd(); + hPeakOverlay->Draw("COLZ"); + cPeakOverlay->Modified(); + cPeakOverlay->Update(); + } + + SaveHistograms( OutRootFileName ); return( 0 ); } -void removeSpikes(int Samples){ - - const float fract = 0.8; - float x, xp, xpp, x3p; - - // assume that there are no spikes - for ( int i = 0; i < Samples; i++) Vcorr[i] = Ameas[i]; - -// find the spike and replace it by mean value of neighbors - 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]; - - if ( x < -5. ){ // a spike candidate - // check consistency with a single channel spike - xp = Ameas[i+1] - N1mean[i+1]; - xpp = Ameas[i+2] - N1mean[i+2]; - 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; - } - else{ - - 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 - } - // treatment for the end of the pipeline must be added !!! - } - } - else{ - // do nothing - } - } // end of spike search and correction - for ( int i = 0; i < Samples; i++ ) debugHistos[ Vcorr_ ].SetBinContent( i, Vcorr[i] ); -} -/* -void computeN1mean( int Samples ){ -cout << "In compute N1mean" << endl; -// compute the mean of the left and right neighbors of a channel - - for( int i = 0; i < Samples; i++){ - if (i == 0){ // use right sample es mean - N1mean[i] = Ameas[i+1]; - } - else if ( i == Samples-1 ){ //use left sample as mean - N1mean[i] = Ameas[i-1]; - } - else{ - N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.; - } - h[tN1mean].SetBinContent(i, Ameas[i] - N1mean[i]); - } -} // end of computeN1mean computation -*/ - -void computeN1mean( int Samples ){ -// compute the mean of the left and right neighbors of a channel - - for( int i = 2; i < Samples - 2; i++){ -/* if (i == 0){ // use right sample es mean - N1mean[i] = Ameas[i+1]; - } - else if ( i == Samples-1 ){ //use left sample as mean - N1mean[i] = Ameas[i-1]; - } - else{ - N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.; - } -*/ - N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.; - } -} // end of computeN1mean computation - - - -/* -// this function takes a vector and tries to apply the DRS calibration to it. -// therfor it uses a DRS calibration file, which must be given as a parameter. -// in case the calibration file can not be opened the function generates a warning on cerr -// and returns an empty vector -// -// the function assumes the first parameter conatins the data of all Pixel -// the RegionOfInterest is deduced from its size(). - - // - -vector * calibrateDrsData( - vector & AllPixelData, - vector & AllStartCells, - const char * DrsCalibFilePath, - bool ApplyCalibToSourceVector = false) -{ - vector * CalibratedDRSData; - - size_t DrsLengthOfVectors; - vector DrsOffset; - vector DrsGain; - vector DrsOffsetSpecialTrigger; - - // the Total Number of Pixel, we have data of - // is deduced from the size() of AllPixelStartCell - unsigned int NumberOfPixel = AllStartCells.size(); - //sanity Check: - if (NumberOfPixel < 1) { - cerr << "calibrateDrsData - Error: NumberOfPixel=AllStartCells.size() ia:" - << AllStartCells.size() << " must be >= 1" << endl; - return CalibratedDRSData - } - - // The RegionOfInterest is deduced from the size() of AllPixel Data - unsigned int RegionOfInterest = AllPixelData.size() / NumberOfPixel; - // sanity Check - if ( RegionOfInterest < 1){ - cerr << "calibrateDrsData - Error: RegionOfInterest = AllPixelData.size() / NumberOfPixel is:" - << RegionOfInterest << " must be >= 1" << endl; - return CalibratedDRSData - } - - // now lets see if we can find the drs calib file - FOpenCalibFile( DrsCalibFilePath, - DrsOffset, - DrsGain, - DrsOffsetSpecialTrigger, - DrsLengthOfVectors); - //sanity check is done in FOpenCalibFile, I guess ... check and maybe TODO - // maybe at least the return code of FOpenCalibFile should be checked here ... TODO - - - - -} - -*/ - -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 ); -} - -// this is a faster version of getValue, that does not do a full calibration -// was used for a performace test... -float correctDrsOffset( int slice, int pixel ){ - - const float dconv = 2000/4096.0; - - // here 1024 is not the RegionOfInterest, but really the lenth of the pipeline - unsigned int physical_slice = ( slice + StartCellVector[ pixel ] ) % 1024; - - unsigned int slice_pt; - unsigned int physical_slice_pt; - slice_pt = pixel * RegionOfInterest + slice; - physical_slice_pt = pixel * RegionOfInterest + physical_slice; - - float vcal = AllPixelDataVector[ slice_pt ] * - dconv - drs_basemean[ physical_slice_pt ]; - return( vcal ); -} - // booking and parameter settings for all histos void BookHistos( ){ - // histograms for baseline extraction - char hName[500]; - char hTitle[500]; - - 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 ); - - h->GetXaxis()->SetTitle( "Sample value (mV)" ); - h->GetYaxis()->SetTitle( "Entries / 0.5 mV" ); - hListBaseline.Add( h ); - hBaseline[i] = h; - } - - hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5); - hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" ); - hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" ); - hList.Add( hMeanBsl ); - - hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5); - hpltMeanBsl->GetXaxis()->SetTitle( "pixel" ); - hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" ); - hList.Add( hpltMeanBsl ); - - hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5); - hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" ); - hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" ); - hList.Add( hRmsBsl ); - - hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5); - hpltRmsBsl->GetXaxis()->SetTitle( "pixel" ); - hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" ); - hList.Add( hpltRmsBsl ); - - hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",1440,-0.5,1439.5, 256, -27.5, 100.5); - hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" ); - hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" ); - hList.Add( hAmplSpek_cfd ); - - hAmplSpek_discri = new TH2F("hAmplSpek_discri","amplitude spektrum - std discriminator",1440,-0.5,1439.5, 256, -27.5, 100.5); - hAmplSpek_discri->GetXaxis()->SetTitle( "pixel" ); - 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 ]; @@ -755,93 +393,20 @@ debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)"); } - hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks", 201, -100.5, 100.5, 1024, -99.5, 100.5); + hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks", + hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 , + 4096, -48.5, 1999.5 ); hPeakOverlay->GetXaxis()->SetTitle( "Timeslices" ); hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" ); - //hList.Add( hPeakOverlay h); + hList.Add( hPeakOverlay ); } -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 - tf.mkdir("BaselineHisto"); - tf.cd("BaselineHisto"); // go to new subdirectory - hListBaseline.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() ; -} - - -///////////////////////////////////////////////////////////////////////////// -/// old BookHistos -/* -void BookHistos( int Samples ){ -// booking and parameter settings for all histos - - h = new TH1F[ Ntypes ]; - - for ( int type = 0; type < Ntypes; type++){ - - h[ type ].SetBins(Samples, 0, Samples); - h[ type ].SetLineColor(1); - h[ type ].SetLineWidth(2); - - // set X axis paras - h[ type ].GetXaxis()->SetLabelSize(0.1); - h[ type ].GetXaxis()->SetTitleSize(0.1); - h[ type ].GetXaxis()->SetTitleOffset(1.2); - h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.)); - - // set Y axis paras - h[ type ].GetYaxis()->SetLabelSize(0.1); - h[ type ].GetYaxis()->SetTitleSize(0.1); - h[ type ].GetYaxis()->SetTitleOffset(0.3); - h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)"); - } - cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",10,10,800,600); - cRawAndSpikeRemoval->Divide(1, 3); - cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600); - cFilter->Divide(1, 3); - - hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024); - -} -*/ -