/** * @file SPeakFinder.cpp * @author Dan Saunders, on behalf of the SoLid collaboration. * @date 17 Feb 2016 */ #include "SPeakFinder.h" //============================================================================== //! Constructor setting up default values. SPeakFinder::SPeakFinder(SDetector * dtr, SClipboard * cb) : ISAlgorithm(dtr, cb, "SPeakFinder"), m_nSaveExamples(5), m_nExamplesSaved(0), m_nPeaksFound(0), m_threshold(50), m_nWaveformsScanned(0), m_integralRange(20), m_shortIntegralRange(15), m_nSamplesRisingEdge(5) { m_options.push_back(new SOptionDouble("Threshold", &m_threshold)); m_options.push_back(new SOptionInt("IntegralRange", &m_integralRange)); m_options.push_back(new SOptionInt("ShortIntegralRange", &m_shortIntegralRange)); m_options.push_back(new SOptionInt("nSamplesRisingEdge", &m_nSamplesRisingEdge)); } //============================================================================== void SPeakFinder::initialize() { std::cout<<"[Note]: Threshold set to: "<histosFile()->mkdir((name() + "/Examples").c_str()); } //============================================================================== //! Peak finding. /*! Currently a simple algorithm. Looks for a bin above some threshold, that * should be set far above noise (say 1 PA) and takes a peak as a bin above * its left and right neighbours. Loop is done over channels. Small number of * examples are kept for visual inspection. */ void SPeakFinder::execute() { for (auto c : (*dtr()->channels())) { for (auto w : (*c->waveforms())) findPeaks(w, c); timeSortReconObjects(c->peaks()); } timeSortReconObjects(cb()->peaks()); } //============================================================================== //! Peak finding for a specific waveform and channel - see execute() void SPeakFinder::findPeaks(SWaveform * wf, SChannel * chan) { std::vector peaks; // Temporary container for peaks found in this waveform. auto start = wf->samples()->begin()+1; auto end = wf->samples()->end()-1-m_integralRange; int iSample=0; for (auto iSamp = start; iSamp != end; iSamp++, iSample++) { // See if the sample is above its neighbours and above threshold. // Allow cases where two touching samples have the same value. if (*iSamp > *(iSamp-1) && *iSamp >= *(iSamp+1) && *iSamp > m_threshold+chan->baseline()) { // We have a peak. SPeak * p = makePeak(wf, chan, iSample, iSamp); cb()->peaks()->push_back(p); chan->peaks()->push_back(p); peaks.push_back(p); m_nPeaksFound++; } } if (m_nExamplesSaved < m_nSaveExamples) saveExample(wf, &peaks); m_nWaveformsScanned++; } //============================================================================== //! Creating the peak object from the waveform. /*! Amplitude and time is currently taken as that of the peak bin. * Integral is taken from the specified range. * SPeakFitter can be used to refine this for better estimates. */ SPeak * SPeakFinder::makePeak(SWaveform * wf, SChannel * chan, int iSample, std::vector::iterator iSamp) { uint time = iSample + wf->time(); double htime = wf->htime() + iSample*chan->sampleWidthNs(); float amplitude = float(wf->samples()->at(iSample)) - chan->baseline(); /* Recall integrals to include rising edge. First arg is start position, then * the num samples to integrate over. If the sample doesn't exist, continue. */ float integral = wf->integrate(iSample - m_nSamplesRisingEdge, m_integralRange, chan->baseline()); float shortIntegral = wf->integrate(iSample - m_nSamplesRisingEdge, m_shortIntegralRange, chan->baseline()); SPeak * p = new SPeak(time, htime, chan, wf); p->setAmplitude(amplitude); p->setIntegral(integral); p->setShortIntegral(shortIntegral); // if (amplitude < 2000 && integral > 60000) h_examples.push_back(wf->getHistogram("eg_A" + std::to_string(amplitude) + "_I" + std::to_string(integral))); p->setWaveformPos(iSamp); // Since iSamp is an iterator. return p; } //============================================================================== void SPeakFinder::finalize() { std::cout<<"[Note]: Number of peaks found: "<histosFile()->cd((name() + "/Examples").c_str()); for (auto h : h_examples) h->Write(); } //============================================================================== //! Saves the waveform and corresponding peaks in histograms. void SPeakFinder::saveExample(SWaveform * wf, std::vector * wf_peaks) { std::string wfName = "wf_" + std::to_string(m_nExamplesSaved); std::string hTitle = "Waveform " + std::to_string(m_nExamplesSaved); hTitle += "; DAQ samples (25 ns); Amplitude (ADC)"; h_examples.push_back(wf->getHistogram(wfName, hTitle)); std::string peakName = "peaksForWf_" + std::to_string(m_nExamplesSaved); std::string peakTitle = "Peaks for waveform " + std::to_string(m_nExamplesSaved); peakTitle += "; DAQ samples (25 ns); Amplitude (ADC)"; TH1D *hPeak = new TH1D(peakName.c_str(), peakTitle.c_str(), wf->samples()->size(), 0, wf->samples()->size()); for (auto p : (*wf_peaks)) hPeak->SetBinContent(p->time()-wf->time()+1, p->amplitude()); h_examples.push_back(hPeak); m_nExamplesSaved++; } //==============================================================================