Index: fact/tools/rootmacros/tpeak.C =================================================================== --- fact/tools/rootmacros/tpeak.C (revision 12366) +++ fact/tools/rootmacros/tpeak.C (revision 12366) @@ -0,0 +1,847 @@ +/********************** Peak-Template *********************** + +this makro shell search for Peaks in data and form a template +of these, so it can be used for detection of other peaks + +*************************************************************/ + +//root libraries + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#define NPIX 1440 +#define NCELL 1024 +#define FAD_MAX_SAMPLES 1024 + +#define HAVE_ZLIB + +//rootmacros + +#include "fits.h" +#include "FOpenCalibFile.c" + +#include "discriminator.h" +#include "discriminator.C" +#include "zerosearch.h" +#include "zerosearch.C" +#include "factfir.C" + + +//----------------------------------------------------------------------------- +// Decleration of variables +//----------------------------------------------------------------------------- + +vector AllPixelDataVector; +vector StartCellVector; + +unsigned int CurrentEventID; + +bool breakout=false; + +size_t ROIxNOP; +UInt_t NumberOfPixels; +UInt_t RegionOfInterest; +int NEvents; + +size_t drs_n; +vector drs_basemean; +vector drs_gainmean; +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 +const int NumberOfDebugHistoTypes = 7; +const unsigned int + Ameas_ = 0, + N1mean_ = 1, + Vcorr_ = 2, + Vtest_ = 3, + Vslide_ = 4, + Vcfd_ = 5, + Vcfd2_ = 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 + +TObjArray hList; +TObjArray hListBaseline; + +void BookHistos( ); +void SaveHistograms( char * ); + +int searchSinglesPeaks = 0; + +//----------------------------------------------------------------------------- +// 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 + ) +{ + +//----------------------------------------------------------------------------- +// Create (pointer to) Canvases, which are used in every run, +// 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); + } + // 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"); + +//----------------------------------------------------------------------------- +// read FACT raw data +// * remove spikes +// * calculate baseline +// * find peaks (CFD and normal discriminator) +// * compute pulse height and pulse integral spektrum of the peaks +//----------------------------------------------------------------------------- + + +//----------------------------------------------------------------------------- +// 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; + vector 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 as2(ks2, 1); + double bs2 = ks2; + +//----------------------------------------------------------------------------- +// prepare datafile +//----------------------------------------------------------------------------- + +// 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); + if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all! + +//Get the DRS calibration + FOpenCalibFile( drsfilename, + drs_basemean, + drs_gainmean, + drs_triggeroffsetmean, + 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 +//----------------------------------------------------------------------------- + for ( int ev = firstevent; ev < firstevent + nevents; ev++) { + // Get an Event --> consists of 1440 Pixel ...erm....data + datafile->GetRow( ev ); + + //------------------------------------- + // Loop over every Pixel of Event + //------------------------------------- + + 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); + } + + // peaks in Ameas[] are found by searching for zero crossings + // in Vcfd2 + // first Argument 1 means ... *rising* edge + // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well + 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] ) ; + } + + //following Code should produce the Overlay of peaks + vector::iterator it; + + 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++){ + + hPeakOverlay->Fill( pos - (it->maxPos), Vcorr[pos]) ; + } + } + + + + 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); + gPad->SetGrid(); + debugHistos[Vcorr_].Draw(); + + cFiltered->cd(1); + 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; + for (unsigned int i=0; isize(); i++){ + OneBox = new TBox( + zXings->at(i).maxPos -10 , + zXings->at(i).maxVal -0.5, + zXings->at(i).maxPos +10 , + zXings->at(i).maxVal +0.5); + OneBox->SetLineColor(kBlue); + OneBox->SetLineWidth(1); + OneBox->SetFillStyle(0); + OneBox->SetFillColor(kRed); + MyBoxes.push_back(OneBox); + OneBox->Draw(); + } + + cFiltered->cd(2); + gPad->SetGrid(); + debugHistos[Vcfd2_].Draw(); + TLine *zeroline = new TLine(0, 0, 1024, 0); + zeroline->SetLineColor(kBlue); + zeroline->Draw(); + + cRawAndSpikeRemoval->Update(); + cFiltered->Update(); + + //OVERLAY PEAKS + cPeakOverlay->cd(); + hPeakOverlay->Draw("COLZ"); + cPeakOverlay->Modified(); + cPeakOverlay->Update(); + + //Process gui events asynchronously during input + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + timer.TurnOn(); + TString input = Getline("Type 'q' to exit, to go on: "); + timer.TurnOff(); + if (input=="q\n") { + breakout=true; + } + + //TODO!!!!!!!!! + // do some Garbage collection ... + // all the Objects on the heap should be deleted here. + + }// 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; + + } + //------------------------------------- + // 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(); + + //OVERLAY PEAKS + cPeakOverlay->cd(); + hPeakOverlay->Draw("COLZ"); + cPeakOverlay->Modified(); + cPeakOverlay->Update(); + +} + + + if (breakout) + break; + } +//------------------------------------- +// 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 ); + + 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 ]; + for ( int type = 0; type < NumberOfDebugHistoTypes; type++){ + debugHistos[ type ].SetBins(1024, 0, 1024); + debugHistos[ type ].SetLineColor(1); + debugHistos[ type ].SetLineWidth(2); + + // set X axis paras + debugHistos[ type ].GetXaxis()->SetLabelSize(0.1); + debugHistos[ type ].GetXaxis()->SetTitleSize(0.1); + debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2); + debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.)); + + // set Y axis paras + debugHistos[ type ].GetYaxis()->SetLabelSize(0.1); + debugHistos[ type ].GetYaxis()->SetTitleSize(0.1); + debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3); + 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->GetXaxis()->SetTitle( "Timeslices" ); + hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" ); + //hList.Add( hPeakOverlay h); + +} + +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 + + // create the histogram file (replace if already existing) + TFile tf( 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); + +} +*/ +