Commit 93b80247 authored by Andrei Muresan's avatar Andrei Muresan
Browse files

Set channel baseline values from Pulse Finding

This is kind of hacky imo, the better way to do this would be a new
EventProcessor altogther (EventPreprocessor?) that would do this sort of
thing in a more generic way. Consider this a placeholder.
parent c9d5d55a
Pipeline #32645 passed with stages
in 2 minutes and 17 seconds
......@@ -3,8 +3,9 @@
#include <cmath>
#include <string>
std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
std::vector<Pulse> PFS::findPulses(Channel &channel, Json::Value cfg) {
std::vector<Pulse> pulses;
Waveform *wf = channel.getWaveform();
int numBins = cfg["NumBins"].asInt();
double riseTimeSigma = cfg["RiseTimeSigma"].asDouble();
......@@ -28,25 +29,25 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
int pulseDiscardCalcCond3NoiseCoefficient =
cfg["PulseDiscardCalcCond3NoiseCoefficient"].asInt();
double baselineMu = calculateBaselineMu(wf, cfg);
double baselineMu = calculateBaselineMu(channel, cfg);
if (!baselineMu) {
// skip noisy waveform
return pulses;
}
numBins = numBins ? numBins : wf.GetNbinsX();
double binWidth = wf.GetBinWidth(1);
numBins = numBins ? numBins : wf->GetNbinsX();
double binWidth = wf->GetBinWidth(1);
int binRiseTime = round(riseTimeSigma * riseTimeSigmaCoefficient / binWidth);
int binFallTime =
round((fallTimeTau * fallTimeTauCoefficient + fallTime2Tau) / binWidth);
for (int bin = binRiseTime + 1; bin <= numBins; bin++) {
double binAmp = wf.GetBinContent(bin);
double binAmp = wf->GetBinContent(bin);
// TODO: rename these condition variables and extract that 5.0 coefficient
bool cond1 = polarity * (binAmp - baselineMu) > (noise * 6.0);
bool cond2 = polarity * (binAmp - wf.GetBinContent(bin - binRiseTime)) >
bool cond2 = polarity * (binAmp - wf->GetBinContent(bin - binRiseTime)) >
(noise * 5.0);
if (!cond1 || !cond2) {
......@@ -58,7 +59,7 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
// double pulseAbsAmp = polarity * binAmp;
// for (int tmpBin = bin + 1; tmpBin <= numBins; tmpBin++) {
// double tmpAbsAmp = polarity * wf.GetBinContent(tmpBin);
// double tmpAbsAmp = polarity * wf->GetBinContent(tmpBin);
// if (pulseAbsAmp < tmpAbsAmp) {
// pulseAbsAmp = tmpAbsAmp;
......@@ -67,9 +68,9 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
double pulseAbsAmp = polarity * binAmp;
bin++;
while (bin <= wf.GetNbinsX() &&
pulseAbsAmp <= polarity * wf.GetBinContent(bin)) {
pulseAbsAmp = polarity * wf.GetBinContent(bin);
while (bin <= wf->GetNbinsX() &&
pulseAbsAmp <= polarity * wf->GetBinContent(bin)) {
pulseAbsAmp = polarity * wf->GetBinContent(bin);
bin++;
}
......@@ -81,7 +82,7 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
for (int tmpBin = bin - binRiseTime - pulseBaselineCalcBeginOffset;
tmpBin < (bin - binRiseTime - pulseBaselineCalcEndOffset); tmpBin++) {
double tmpBinAmp = wf.GetBinContent(tmpBin);
double tmpBinAmp = wf->GetBinContent(tmpBin);
double absTmpBinAmp = polarity * tmpBinAmp;
pulseBaseline += tmpBinAmp;
......@@ -97,13 +98,13 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
double charge{0.0};
for (int tmpBin = (bin - binRiseTime - pulseChargeCalcOffset);
tmpBin < (bin + binFallTime); tmpBin++) {
charge += wf.GetBinContent(tmpBin);
charge += wf->GetBinContent(tmpBin);
}
charge -= baselineMu * (binRiseTime + binFallTime + pulseChargeCalcOffset);
charge *= polarity;
charge *= wf.GetBinWidth(1);
charge *= wf->GetBinWidth(1);
// find the pulse's starting point
int pulseStartBin;
......@@ -111,7 +112,7 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
pulseStartBin < (bin + binFallTime); pulseStartBin++) {
double noiseThreshold = pulseStartCalcNoiseCoefficient * noise;
double absAmpDelta =
polarity * (wf.GetBinContent(pulseStartBin) - pulseBaseline);
polarity * (wf->GetBinContent(pulseStartBin) - pulseBaseline);
if (absAmpDelta >= noiseThreshold) {
break;
}
......@@ -125,7 +126,7 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
pulseEndBin--) {
double noiseThreshold = pulseEndCalcNoiseCoefficient * noise;
double absAmpDelta =
polarity * (wf.GetBinContent(pulseStartBin) - pulseBaseline);
polarity * (wf->GetBinContent(pulseStartBin) - pulseBaseline);
if (absAmpDelta >= noiseThreshold) {
break;
......@@ -135,7 +136,7 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
int pulseWidth = pulseEndBin - pulseStartBin;
// TODO: figure out what's going on here and give these better names
bool discardCond1 = pulseStartBin < wf.GetNbinsX();
bool discardCond1 = pulseStartBin < wf->GetNbinsX();
bool discardCond2 = charge > pulseDiscardCalcCond2NoiseCoefficient * noise *
std::sqrt(binRiseTime + binFallTime);
bool discardCond3 = (-1 * prePulseMin - baselineMu) <
......@@ -145,7 +146,7 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
bin--;
const Pulse pulse{
wf.GetBinCenter(bin), // time
wf->GetBinCenter(bin), // time
pulseAbsAmp, // absAmp
polarity * (pulseAbsAmp - pulseBaseline), // amp
pulseBaseline, // baseline
......@@ -162,15 +163,17 @@ std::vector<Pulse> PFS::findPulses(const Waveform &wf, Json::Value cfg) {
return pulses;
}
double PFS::calculateBaselineMu(const Waveform &wf, Json::Value cfg) {
double PFS::calculateBaselineMu(Channel &channel, Json::Value cfg) {
double ampMax = cfg["MaxAmplitude"].asDouble();
double ampMin = cfg["MinAmplitude"].asDouble();
double wfAmpBinning = cfg["WaveformAmplitudeBinning"].asDouble();
double noise = cfg["Noise"].asDouble();
double noiseCoefficient = cfg["NoiseCoefficient"].asDouble();
Waveform *wf = channel.getWaveform();
const std::string baselineHistName =
std::string{wf.GetName()} + std::string{"Baseline"};
std::string{wf->GetName()} + std::string{"Baseline"};
TH1F *baselineHist;
......@@ -187,14 +190,17 @@ double PFS::calculateBaselineMu(const Waveform &wf, Json::Value cfg) {
// baseline histogram represents the distribution of amplitude values
// across the waveform
for (int bin = 1; bin <= wf.GetNbinsX(); bin++) {
baselineHist->Fill(wf.GetBinContent(bin));
for (int bin = 1; bin <= wf->GetNbinsX(); bin++) {
baselineHist->Fill(wf->GetBinContent(bin));
}
// calculate baseline mu value
double baselineMu =
baselineHist->GetBinLowEdge(baselineHist->GetMaximumBin());
channel.setBaselineMu(baselineMu);
channel.setBaselineRMS(baselineHist->GetRMS());
if (baselineHist->GetRMS() >= noise * noiseCoefficient) {
return 0.0;
}
......
......@@ -5,11 +5,11 @@
class PFS : public IPulseFindingStrategy {
public:
std::vector<Pulse> findPulses(const Waveform &wf, Json::Value cfg) override;
std::vector<Pulse> findPulses(Channel &channel, Json::Value cfg) override;
bool checkParams(Json::Value params) override;
private:
double calculateBaselineMu(const Waveform &wf, Json::Value cfg);
double calculateBaselineMu(Channel &channel, Json::Value cfg);
};
#endif // __PFS_H_
......@@ -28,9 +28,8 @@ int PulseFindingProcessor::processEvent(Event *event) {
auto channels = event->getChannels();
for (int i = 0; i < channels->size(); i++) {
auto channel = channels[i];
auto pulses =
strategies[i]->findPulses(*(channels->at(i).getWaveform()), vars[i]);
auto channel = channels->at(i);
auto pulses = strategies[i]->findPulses(channels->at(i), vars[i]);
channels->at(i).setPulses(pulses);
pulseCounter += pulses.size();
......
#ifndef __PULSEFINDINGSTRATEGY_H
#define __PULSEFINDINGSTRATEGY_H
#include "../Pulse.h"
#include "../Waveform.h"
#include "../Channel.h"
#include "json/json.h"
#include <vector>
class IPulseFindingStrategy {
public:
virtual std::vector<Pulse> findPulses(const Waveform &wf,
Json::Value cfg) = 0;
virtual std::vector<Pulse> findPulses(Channel &channel, Json::Value cfg) = 0;
virtual bool checkParams(Json::Value params) = 0;
virtual ~IPulseFindingStrategy() {}
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment