From ebaf1d5f2db17618e23fce33e33d8d4cf91a7bac Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 18 Feb 2022 18:40:16 +0000
Subject: [PATCH] Update digitisation to have realistic normalisation and
baseline. Update reco with --isMC arg to run on digitised MC
---
.../CaloDigiAlgs/python/CaloDigiAlgsConfig.py | 4 +++
.../CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx | 7 ++++--
.../CaloDigiAlgs/src/CaloWaveformDigiAlg.h | 6 ++---
.../Reconstruction/scripts/faser_reco.py | 25 +++++++++++++------
.../python/ScintDigiAlgsConfig.py | 14 +++++++----
.../src/ScintWaveformDigiAlg.cxx | 8 +++---
.../ScintDigiAlgs/src/ScintWaveformDigiAlg.h | 3 +++
.../WaveDigiTools/IWaveformDigitisationTool.h | 10 +++++++-
.../IWaveformDigitisationTool.icc | 16 +++++++++---
.../src/WaveformDigitisationTool.cxx | 1 +
.../src/WaveformReconstructionTool.cxx | 1 +
11 files changed, 70 insertions(+), 25 deletions(-)
diff --git a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
index cd0c9a45..ce45e625 100644
--- a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
+++ b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
@@ -36,6 +36,10 @@ def CaloWaveformDigiCfg(flags, name="CaloWaveformDigiAlg", **kwargs):
digiAlg.CB_n = 10
digiAlg.CB_sigma = 4
digiAlg.CB_mean = 820
+ digiAlg.CB_norm = 2
+
+ digiAlg.base_mean = 15000
+ digiAlg.base_rms = 3
acc.addEventAlgo(digiAlg)
diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
index 8e60aefa..1abfcf9f 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
@@ -4,6 +4,7 @@
#include <vector>
#include <map>
+#include <utility>
CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name,
ISvcLocator* pSvcLocator)
@@ -30,12 +31,13 @@ CaloWaveformDigiAlg::initialize() {
// TODO: Change params compared to scint
// m_kernel = new TF1("PDF", " ROOT::Math::crystalball_pdf(x, -0.9, 10, 4, 900)", 0, 1200);
- m_kernel = new TF1("PDF", "ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
+ m_kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
//m_kernel->SetParameters(-0.25,10,4,900);
m_kernel->SetParameter(0, m_CB_alpha);
m_kernel->SetParameter(1, m_CB_n);
m_kernel->SetParameter(2, m_CB_sigma);
m_kernel->SetParameter(3, m_CB_mean);
+ m_kernel->SetParameter(4, m_CB_norm);
return StatusCode::SUCCESS;
@@ -76,7 +78,8 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
// Digitise the hits
CHECK( m_digiTool->digitise<CaloHitCollection>(caloHitHandle.ptr(),
- waveformContainerHandle.ptr(), m_kernel) );
+ waveformContainerHandle.ptr(), m_kernel,
+ std::pair<float, float>(m_base_mean, m_base_rms)) );
ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
index 9de68257..f99d13b5 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
@@ -52,10 +52,10 @@ class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
Gaudi::Property<double> m_CB_n {this, "CB_n", 0, "n of the crystal ball function"};
Gaudi::Property<double> m_CB_mean {this, "CB_mean", 0, "Mean of the crystal ball function"};
Gaudi::Property<double> m_CB_sigma {this, "CB_sigma", 0, "Sigma of the crystal ball function"};
+ Gaudi::Property<double> m_CB_norm {this, "CB_norm", 0, "Norm of the crystal ball function"};
-
-
-
+ Gaudi::Property<double> m_base_mean {this, "base_mean", 0, "Mean of the baseline"};
+ Gaudi::Property<double> m_base_rms {this, "base_rms", 0, "RMS of the baseline"};
/// Kernel PDF
TF1* m_kernel;
diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 9810f960..1579f4c1 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -30,6 +30,9 @@ parser.add_argument("-v", "--verbose", action='store_true',
help="Turn on DEBUG output")
parser.add_argument("--clusterFit", action='store_true',
help="Use ClusterFit (old) track finder - default: SegmentFit(new)")
+parser.add_argument("--isMC", action='store_true',
+ help="Running on digitised MC rather than data")
+
args = parser.parse_args()
@@ -86,7 +89,7 @@ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
Configurable.configurableRun3Behavior = True
# Flags for this job
-ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ConfigFlags.Input.isMC = args.isMC # Needed to bypass autoconfig
ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
ConfigFlags.Input.ProjectName = "data20"
@@ -139,8 +142,13 @@ acc.merge(PoolWriteCfg(ConfigFlags))
#
# Set up RAW data access
-from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg
-acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags))
+
+if args.isMC:
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ acc.merge(PoolReadCfg(ConfigFlags))
+else:
+ from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg
+ acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags))
#
# Needed, or move to MainServicesCfg?
@@ -204,11 +212,12 @@ print( "Writing out xAOD objects:" )
print( acc.getEventAlgo("OutputStreamxAOD").ItemList )
# Hack to avoid problem with our use of MC databases when isMC = False
-replicaSvc = acc.getService("DBReplicaSvc")
-replicaSvc.COOLSQLiteVetoPattern = ""
-replicaSvc.UseCOOLSQLite = True
-replicaSvc.UseCOOLFrontier = False
-replicaSvc.UseGeomSQLite = True
+if not args.isMC:
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
# Configure verbosity
# ConfigFlags.dump()
diff --git a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
index 243e4493..7f2ba853 100644
--- a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
+++ b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
@@ -10,10 +10,10 @@ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
# https://indico.cern.ch/event/1099652/contributions/4626975/attachments/2352595/4013927/Faser-Physics-run3933-plots.pdf (20/01/2022)
# Parameters are per scintillator source, but not per channel.
dict_CB_param = {}
-dict_CB_param["Trigger"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7)
-dict_CB_param["Timing"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3) # copy from Preshower; Timing was not in TestBeam
-dict_CB_param["Veto"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7) # copy from Trigger; Veto was not in TestBeam, but in sim "Veto" is the TestBeam Trigger component
-dict_CB_param["Preshower"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3)
+dict_CB_param["Trigger"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7, CB_norm = 500 )
+dict_CB_param["Timing"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3, CB_norm = 500) # copy from Preshower; Timing was not in TestBeam
+dict_CB_param["Veto"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7, CB_norm = 1000) # copy from Trigger; Veto was not in TestBeam, but in sim "Veto" is the TestBeam Trigger component
+dict_CB_param["Preshower"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3, CB_norm = 500)
# One stop shopping for normal FASER data
def ScintWaveformDigitizationCfg(flags):
@@ -46,7 +46,11 @@ def ScintWaveformDigiCfg(flags, name="ScintWaveformDigiAlg", source="", **kwargs
digiAlg.CB_n = dict_CB_param[source]["CB_n"]
digiAlg.CB_mean = dict_CB_param[source]["CB_mean"]
digiAlg.CB_sigma = dict_CB_param[source]["CB_sigma"]
-
+ digiAlg.CB_norm = dict_CB_param[source]["CB_norm"]
+
+ digiAlg.base_mean = 15000
+ digiAlg.base_rms = 3
+
kwargs.setdefault("WaveformDigitisationTool", tool)
acc.addEventAlgo(digiAlg)
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
index b2508041..91cffd93 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -5,6 +5,7 @@
#include <vector>
#include <map>
+#include <utility>
ScintWaveformDigiAlg::ScintWaveformDigiAlg(const std::string& name,
@@ -33,12 +34,12 @@ ScintWaveformDigiAlg::initialize() {
- m_kernel = new TF1("PDF", "ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
+ m_kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
m_kernel->SetParameter(0, m_CB_alpha);
m_kernel->SetParameter(1, m_CB_n);
m_kernel->SetParameter(2, m_CB_sigma);
m_kernel->SetParameter(3, m_CB_mean);
-
+ m_kernel->SetParameter(4, m_CB_norm);
return StatusCode::SUCCESS;
}
@@ -77,7 +78,8 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
// Digitise the hits
CHECK( m_digiTool->digitise<ScintHitCollection>(scintHitHandle.ptr(),
- waveformContainerHandle.ptr(), m_kernel) );
+ waveformContainerHandle.ptr(), m_kernel,
+ std::pair<float, float>(m_base_mean, m_base_rms)) );
ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
index af60e02d..9bf7843a 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
@@ -54,7 +54,10 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
Gaudi::Property<double> m_CB_n {this, "CB_n", 0, "n of the crystal ball function"};
Gaudi::Property<double> m_CB_mean {this, "CB_mean", 0, "Mean of the crystal ball function"};
Gaudi::Property<double> m_CB_sigma {this, "CB_sigma", 0, "Sigma of the crystal ball function"};
+ Gaudi::Property<double> m_CB_norm {this, "CB_norm", 0, "Norm of the crystal ball function"};
+ Gaudi::Property<double> m_base_mean {this, "base_mean", 0, "Mean of the baseline"};
+ Gaudi::Property<double> m_base_rms {this, "base_rms", 0, "RMS of the baseline"};
/// Kernel PDF
TF1* m_kernel;
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
index 3e85ff83..7241e816 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
@@ -24,6 +24,9 @@
#include "WaveRawEvent/RawWaveform.h"
#include "TF1.h"
+#include "TRandom3.h"
+
+#include <utility>
///Interface for waveform digitisation tools
class IWaveformDigitisationTool : virtual public IAlgTool
@@ -42,11 +45,16 @@ public:
// Digitise HITS to Raw waveform
template<class CONT>
StatusCode digitise(const CONT* hitCollection,
- RawWaveformContainer* waveContainer, TF1* kernel) const;
+ RawWaveformContainer* waveContainer,
+ TF1* kernel, std::pair<float, float> base
+ ) const;
private:
ServiceHandle<IMessageSvc> m_msgSvc;
+protected:
+ TRandom3* m_random;
+
};
#include "WaveDigiTools/IWaveformDigitisationTool.icc"
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
index 57d4839b..8d9a466d 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
@@ -3,7 +3,8 @@
template<class CONT>
StatusCode IWaveformDigitisationTool::digitise(const CONT* hitCollection,
- RawWaveformContainer* container, TF1* kernel) const {
+ RawWaveformContainer* container, TF1* kernel,
+ std::pair<float, float> base) const {
// Check the container
@@ -18,12 +19,16 @@ StatusCode IWaveformDigitisationTool::digitise(const CONT* hitCollection,
for (unsigned int i=0; i<size; i++) time[i] = 2.*i;
std::map<unsigned int, std::vector<uint16_t>> waveforms;
- unsigned int baseline = 8000; // TODO: vary this + add noise
+ //unsigned int baseline = 800;
+
+ // TODO: varying time to centre of edge of bin (odd = centre, even = edge)
// Loop over time samples
for (const auto& t : time) {
std::map<unsigned int, float> counts;
+ unsigned int baseline = m_random->Gaus(base.first, base.second);
+
// Convolve hit energy with kernel and sum for each ID (i.e. channel)
for (const auto& hit : *hitCollection) {
counts[hit.identify()] += kernel->Eval(t) * hit.energyLoss();
@@ -32,7 +37,12 @@ StatusCode IWaveformDigitisationTool::digitise(const CONT* hitCollection,
// Add count to correct waveform vec
for (const auto& c : counts) {
- waveforms[c.first].push_back(baseline - c.second);
+ int value = baseline - c.second;
+ if (value < 0) {
+ log << MSG::WARNING << "Found pulse " << c.second << " larger than baseline " << c.first << endmsg;
+ value = 0; // Protect against scaling signal above baseline
+ }
+ waveforms[c.first].push_back(value);
//std::cout << "ADC " << c.first << " @ " << t << ": " << baseline - c.second << std::endl;
}
}
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
index e4776da3..1c71fe19 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
@@ -20,6 +20,7 @@ WaveformDigitisationTool::WaveformDigitisationTool(const std::string& type, cons
StatusCode
WaveformDigitisationTool::initialize() {
ATH_MSG_INFO( name() << "::initalize()" );
+ m_random = new TRandom3();
return StatusCode::SUCCESS;
}
diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
index 7e564dae..ce9504ee 100644
--- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
+++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
@@ -646,6 +646,7 @@ WaveformReconstructionTool::findRawHitValues(const std::vector<float> time, cons
double sum = 0.;
double sum2 = 0.;
for (unsigned int i=0; i<time.size(); i++) {
+ //ATH_MSG_DEBUG("findRawHitValues Time: " << time[i] << " Wave: " << wave[i]);
tot += wave[i];
sum += time[i] * wave[i];
sum2 += time[i] * time[i] * wave[i];
--
GitLab