From 18c5debfdc93e3a4d645a9ef70a08aeeee7ed666 Mon Sep 17 00:00:00 2001
From: Nicholas Styles <nicholas.styles@desy.de>
Date: Thu, 29 Oct 2020 10:01:06 +0000
Subject: [PATCH] Merge branch 'improveVMMSpeed' into '21.9'
Improve VMM CPU performance
See merge request atlas/athena!37699
---
.../MM_Digitization/VMM_Shaper.h | 24 +--
.../src/MM_ElectronicsResponseSimulation.cxx | 2 +-
.../MM_Digitization/src/VMM_Shaper.cxx | 163 ++++++++++++------
3 files changed, 125 insertions(+), 64 deletions(-)
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/VMM_Shaper.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/VMM_Shaper.h
index 3798409d2d77..4def653f55ca 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/VMM_Shaper.h
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/VMM_Shaper.h
@@ -5,40 +5,40 @@
#ifndef MM_DIGITIZATION_VMM_SHAPER_H
#define MM_DIGITIZATION_VMM_SHAPER_H
-#include <TH1F.h>
#include <vector>
class VMM_Shaper{
public:
- VMM_Shaper(float peakTime);
+ VMM_Shaper(const float peakTime, const float lowerTimeWindow, const float upperTimeWindow);
virtual ~VMM_Shaper()=default;
void initialize();
- void vmmPeakResponse(const std::vector<float> effectiveCharge, const std::vector<float> electronsTime, const double electronicsThreshold, double &litudeFirstPeak, double &timeFirstPeak) const;
- void vmmThresholdResponse(const std::vector<float> effectiveCharge, const std::vector<float> electronsTime, const double electronicsThreshold, double &litudeAtFirstPeak, double &timeAtThreshold) const;
+ void vmmPeakResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold, double &litudeFirstPeak, double &timeFirstPeak) const;
+ void vmmThresholdResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold, double &litudeAtFirstPeak, double &timeAtThreshold) const;
- bool hasChargeAboveThreshold(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double &electronicsThreshold) const;
+ bool hasChargeAboveThreshold(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold) const;
private:
double m_peakTime;
- double m_timeStep;
- double m_maxTime;
-
+ double m_lowerTimeWindow, m_upperTimeWindow;
+
+ double m_timeStep, m_inverseTimeStep;
+ double m_preCalculationVMMShaper;
// shaper parameters
double m_a, m_pole0, m_re_pole1, m_im_pole1, m_pole1_square, m_k1_abs, m_argK1;
int m_nBins;
double m_peakTimeChargeScaling;
- void vmmResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, TH1F &response) const;
- int findPeak(const TH1F &response, const double &electronicsThreshold) const;
-
+ double vmmResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, double time) const;
+ double findPeak(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold) const;
+ bool aboveThresholdSimple(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold) const;
};
-#endif // MM_DIGITIZATION_VMM_SHAPER_H
\ No newline at end of file
+#endif // MM_DIGITIZATION_VMM_SHAPER_H
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_ElectronicsResponseSimulation.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_ElectronicsResponseSimulation.cxx
index e227749c3431..c47bbbbcf8a1 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_ElectronicsResponseSimulation.cxx
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_ElectronicsResponseSimulation.cxx
@@ -29,7 +29,7 @@ MM_ElectronicsResponseSimulation::MM_ElectronicsResponseSimulation():
/*******************************************************************************/
void MM_ElectronicsResponseSimulation::initialize()
{
- m_vmmShaper = std::make_unique<VMM_Shaper>(m_peakTime);
+ m_vmmShaper = std::make_unique<VMM_Shaper>(m_peakTime, m_timeWindowLowerOffset, m_timeWindowUpperOffset);
m_vmmShaper->initialize();
}
/*******************************************************************************/
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/VMM_Shaper.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/VMM_Shaper.cxx
index d27d7e2dfdf5..39d0763c2ab7 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/VMM_Shaper.cxx
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/VMM_Shaper.cxx
@@ -4,6 +4,7 @@
#include "MM_Digitization/VMM_Shaper.h"
#include <cmath>
+#include <algorithm>
namespace {
// VMM shaper parameters provided by G. Iakovidis
@@ -16,10 +17,12 @@ namespace {
static constexpr double mmIonFlowTime = 150.; // ns
}
-VMM_Shaper::VMM_Shaper(float peakTime):m_peakTime(peakTime),
-m_timeStep(0.1),
-m_maxTime(700.0)
+VMM_Shaper::VMM_Shaper(const float peakTime,const float lowerTimeWindow,const float upperTimeWindow):m_peakTime(peakTime),
+m_lowerTimeWindow(lowerTimeWindow),
+m_upperTimeWindow(upperTimeWindow),
+m_timeStep(0.1)
{
+ m_inverseTimeStep = 1./m_timeStep;
initialize();
}
@@ -33,80 +36,138 @@ void VMM_Shaper::initialize() {
m_k1_abs = std::sqrt(Re_K1*Re_K1 + Im_K1*Im_K1);
m_argK1 = std::atan2(Im_K1, Re_K1);
- m_nBins = m_maxTime/m_timeStep;
// scale factor for charge taking into account the mm ion flow time of ~150ns
// if the peaking time is lower then that, only a fration of the total charge is integrated
m_peakTimeChargeScaling = (m_peakTime < mmIonFlowTime ? 1.0*m_peakTime/mmIonFlowTime : 1.0);
+
+ // preCalculate factor to avoid recalculating for each electron
+ m_preCalculationVMMShaper = chargeScaleFactor*m_peakTimeChargeScaling*std::pow(m_a, 3)*m_pole0*m_pole1_square;
}
-void VMM_Shaper::vmmResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, TH1F &response) const {
- for (uint i_electron = 0; i_electron < effectiveCharge.size(); i_electron++) {
- for (double ti = electronsTime[i_electron]; ti < m_maxTime; ti += m_timeStep) {
- double t = (ti-electronsTime.at(i_electron))*(10^-9);
+double VMM_Shaper::vmmResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, double time) const{
+ double response = 0;
+ for (unsigned int i_electron = 0; i_electron < effectiveCharge.size(); i_electron++) {
+ if (time < electronsTime.at(i_electron)) continue;
+ double t = (time-electronsTime.at(i_electron))*(10^-9);
// now follows the vmm shaper response function provided by G. Iakovidis
// It is described in section 7.1.3 of https://cds.cern.ch/record/1955475
- double st = effectiveCharge.at(i_electron)*chargeScaleFactor*m_peakTimeChargeScaling*std::pow(m_a, 3)*m_pole0*m_pole1_square*((K0*std::exp(-t*m_pole0))+(2.*m_k1_abs*std::exp(-t*m_re_pole1)*std::cos(-t*m_im_pole1+m_argK1)));
- response.Fill(ti, st);
- }
+ double st = effectiveCharge.at(i_electron)*m_preCalculationVMMShaper*((K0*std::exp(-t*m_pole0))+(2.*m_k1_abs*std::exp(-t*m_re_pole1)*std::cos(-t*m_im_pole1+m_argK1)));
+ response += st;
}
+ return response;
}
-void VMM_Shaper::vmmPeakResponse(const std::vector<float> effectiveCharge, const std::vector<float> electronsTime, const double electronicsThreshold, double &litudeFirstPeak, double &timeFirstPeak) const {
- TH1F response("response", "response", m_nBins, 0, m_maxTime);
+void VMM_Shaper::vmmPeakResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold, double &litudeFirstPeak, double &timeFirstPeak) const{
+ double t_peak = findPeak(effectiveCharge, electronsTime, electronicsThreshold);
+
+ if (t_peak == -9999 ) return; // no peak found
- vmmResponse(effectiveCharge, electronsTime, response);
- int i_bin_peak = findPeak(response, electronicsThreshold);
+ amplitudeFirstPeak = vmmResponse(effectiveCharge, electronsTime, t_peak);
+ timeFirstPeak = t_peak;
+}
- if (i_bin_peak == -1) return; // no peak found
- timeFirstPeak = response.GetXaxis()->GetBinCenter(i_bin_peak);
- amplitudeFirstPeak = response.GetBinContent(i_bin_peak);
-}
+void VMM_Shaper::vmmThresholdResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold, double &litudeAtFirstPeak, double &timeAtThreshold) const {
+ if (!aboveThresholdSimple(effectiveCharge, electronsTime, electronicsThreshold)) return;
+ if (effectiveCharge.size() == 0) return; // protect min_element
+ double startTime = m_lowerTimeWindow;
+ double minElectronTime = *std::min_element(electronsTime.begin(), electronsTime.end());
+ if (startTime < minElectronTime) startTime = minElectronTime; // if smallest strip times are higher then the lower time window, just start the loop from the smallest electron time
+
+ double tmpTimeAtThreshold = -9999;
+ for (double time = startTime; time < m_upperTimeWindow; time += m_timeStep) {
+ if (vmmResponse(effectiveCharge, electronsTime, time) >= electronicsThreshold) {
+ tmpTimeAtThreshold = time;
+ break;
+ }
+ }
-void VMM_Shaper::vmmThresholdResponse(const std::vector<float> effectiveCharge, const std::vector<float> electronsTime, const double electronicsThreshold, double &litudeAtFirstPeak, double &timeAtThreshold) const {
- TH1F response("response", "response", m_nBins, 0, m_maxTime);
+ if (tmpTimeAtThreshold == -9999) return;
- vmmResponse(effectiveCharge, electronsTime, response);
- int i_bin_peak = findPeak(response, electronicsThreshold);
+ double t_peak = findPeak(effectiveCharge, electronsTime, electronicsThreshold);
+ if (t_peak == -9999) return;
- int binAboveThreshold = response.FindFirstBinAbove(electronicsThreshold);
- if (binAboveThreshold == -1) return;
- timeAtThreshold = response.GetXaxis()->GetBinCenter(binAboveThreshold);
- amplitudeAtFirstPeak = response.GetBinContent(i_bin_peak);
+ timeAtThreshold = tmpTimeAtThreshold;
+ amplitudeAtFirstPeak = vmmResponse(effectiveCharge, electronsTime, t_peak);
}
-int VMM_Shaper::findPeak(const TH1F &response, const double &electronicsThreshold) const{
- TH1F derivative("derivative", "derivative", m_nBins, 0, m_maxTime);
+double VMM_Shaper::findPeak(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold) const{
+ if(effectiveCharge.size()==0) return -9999; // protect min_element
+ double startTime = m_lowerTimeWindow;
+ double minElectronTime = *std::min_element(electronsTime.begin(), electronsTime.end());
+ if(startTime < minElectronTime) startTime = minElectronTime; // if smallest strip times are higher then the lower time window, just start the loop from the smallest electron time
+
+ double oldResponse = 0;
+ double oldDerivative = 0 , currentDerivative = 0;
+
+ for (double time = startTime; time < m_upperTimeWindow; time += m_timeStep) {
+
+ double response = vmmResponse(effectiveCharge, electronsTime, time);
+
+ oldDerivative = currentDerivative;
+ currentDerivative = (response-oldResponse) * m_inverseTimeStep;
+ // check if sign of derivative has not changed ==> no peak; or if response is below threshold
+ if (oldDerivative*currentDerivative >= 0 || oldResponse < electronicsThreshold) {oldResponse = response; continue;}
+
+ // from here one its assumed that a peak above threshold was found
+
+ // check if the derivative is monoton falling within the given window of 5 bins
+ bool checkDerivative = true;
+ double tmp_checkOldDerivative = 0, tmp_checkCurrentDerivative = 0;
+ double tmp_checkOldResponse = 0;
- // First the derivative gets calculated, then its zero crossing is searched for
- // if in 4 bins around the 0 crossing the derivative is falling, it is accepted as peak
+ int searchWindow = 5;
- for(int i_bin = 1; i_bin<m_nBins; i_bin++){
- derivative.SetBinContent(i_bin, (response.GetBinContent(i_bin+1) - response.GetBinContent(i_bin)) / response.GetXaxis()->GetBinWidth(i_bin));
+ for (int i_timeOfPeak = -searchWindow; i_timeOfPeak <= searchWindow; i_timeOfPeak ++) {
+ double response = vmmResponse(effectiveCharge, electronsTime, time + i_timeOfPeak * m_timeStep);
+
+ tmp_checkOldDerivative = tmp_checkCurrentDerivative;
+ tmp_checkCurrentDerivative = (response - tmp_checkOldResponse) * m_inverseTimeStep;
+ tmp_checkOldResponse = response;
+
+ if (i_timeOfPeak >= -searchWindow + 2 // needs two iterations to fill the variables
+ && tmp_checkOldDerivative < tmp_checkCurrentDerivative) { // derivative is not falling monothonic
+ checkDerivative = false;
+ break;
+ }
+ }
+ if (!checkDerivative) continue;
+ return time - m_timeStep;
}
+ return -9999; // no peak found
+}
+
+
+bool VMM_Shaper::hasChargeAboveThreshold(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold) const{
+ if (!aboveThresholdSimple(effectiveCharge, electronsTime, electronicsThreshold)) return false;
- for (int i_bin = 1; i_bin < m_nBins; i_bin++) {
- if (derivative.GetBinContent(i_bin) * derivative.GetBinContent(i_bin + 1) > 0.0) continue; // continue loop if no 0 crossing was there
- if (response.GetBinContent(i_bin + 1) < electronicsThreshold) continue; // continue if peak is below threshold
- bool derivativeCut = true;
- for (int j_bin = i_bin-2; j_bin <= i_bin + 2; j_bin++) {
- if (derivative.GetBinContent(j_bin) <= derivative.GetBinContent(j_bin+1)) { // check that the derivative is falling for 4 bins around the 0 crossing
- derivativeCut = false;
- }
- if (derivativeCut) {
- return i_bin + 1;
- }
+ if (effectiveCharge.size() == 0) return false; // protect min_element
+ double startTime = m_lowerTimeWindow;
+ double minElectronTime = *std::min_element(electronsTime.begin(), electronsTime.end());
+ // since we are only checking if signal is above threshold, we can start searching close to the peak
+ minElectronTime += m_peakTime * 0.8;
+ if (startTime < minElectronTime) startTime = minElectronTime; // if smallest strip times are higher then the lower time window, just start the loop from the smallest electron time
+
+ for (double time = startTime; time < m_upperTimeWindow; time += m_timeStep) {
+ if (vmmResponse(effectiveCharge, electronsTime, time) >= electronicsThreshold) {
+ return true;
}
}
- return -1;
+ return false;
+
}
-bool VMM_Shaper::hasChargeAboveThreshold(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double &electronicsThreshold) const {
- TH1F response("response", "response", m_nBins, 0, m_maxTime);
- vmmResponse(effectiveCharge, electronsTime, response);
- int i_aboveThreshold = response.FindFirstBinAbove(electronicsThreshold);
- return i_aboveThreshold > 0;
-}
\ No newline at end of file
+bool VMM_Shaper::aboveThresholdSimple(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold) const {
+ // check if total strip charge is above threshold, otherwise skip VMM
+ float chargeSum = 0;
+ for (unsigned int i_elec = 0; i_elec < effectiveCharge.size(); i_elec++) {
+ if (electronsTime.at(i_elec) >= m_lowerTimeWindow - m_peakTime && electronsTime.at(i_elec) <= m_upperTimeWindow) {
+ chargeSum += effectiveCharge.at(i_elec)*m_peakTimeChargeScaling;
+ }
+ }
+ if (chargeSum >= electronicsThreshold) return true;
+ return false;
+}
--
GitLab