diff --git a/Calorimeter/CaloDigiAlgs/CMakeLists.txt b/Calorimeter/CaloDigiAlgs/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..42991aa03cc03af6e51ab161d36d1d61df1a7956
--- /dev/null
+++ b/Calorimeter/CaloDigiAlgs/CMakeLists.txt
@@ -0,0 +1,15 @@
+################################################################################
+# Package: CaloDigiAlgs
+################################################################################
+
+# Declare the package name:
+atlas_subdir( CaloDigiAlgs )
+
+# Component(s) in the package:
+atlas_add_component( CaloDigiAlgs
+ src/*.cxx src/*.h
+ src/components/*.cxx
+ LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent FaserCaloSimEvent WaveDigiToolsLib)
+
+atlas_install_python_modules( python/*.py )
+
diff --git a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..cd0c9a454139d5906dd2e01ee0d1cc09b1697f76
--- /dev/null
+++ b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
@@ -0,0 +1,55 @@
+# Copyright (C) 2020-2021 CERN for the benefit of the FASER collaboration
+
+from re import VERBOSE
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+
+# One stop shopping for normal FASER data
+def CaloWaveformDigitizationCfg(flags):
+ """ Return all algorithms and tools for Waveform digitization """
+ acc = ComponentAccumulator()
+
+ if not flags.Input.isMC:
+ return acc
+
+ acc.merge(CaloWaveformDigiCfg(flags, "CaloWaveformDigiAlg"))
+ acc.merge(CaloWaveformDigitizationOutputCfg(flags))
+ return acc
+
+# Return configured digitization algorithm from SIM hits
+def CaloWaveformDigiCfg(flags, name="CaloWaveformDigiAlg", **kwargs):
+
+ acc = ComponentAccumulator()
+
+ tool = CompFactory.WaveformDigitisationTool(name="CaloWaveformDigtisationTool", **kwargs)
+
+ kwargs.setdefault("CaloHitContainerKey", "EcalHits")
+ kwargs.setdefault("WaveformContainerKey", "CaloWaveforms")
+
+ digiAlg = CompFactory.CaloWaveformDigiAlg(name, **kwargs)
+ kwargs.setdefault("WaveformDigitisationTool", tool)
+
+ digiAlg.CB_alpha = -0.9
+ digiAlg.CB_n = 10
+ digiAlg.CB_sigma = 4
+ digiAlg.CB_mean = 820
+
+
+ acc.addEventAlgo(digiAlg)
+
+ return acc
+
+def CaloWaveformDigitizationOutputCfg(flags, **kwargs):
+ """ Return ComponentAccumulator with output for Waveform Digi"""
+ acc = ComponentAccumulator()
+ ItemList = [
+ "RawWaveformContainer#*"
+ ]
+ acc.merge(OutputStreamCfg(flags, "RDO"))
+ ostream = acc.getEventAlgo("OutputStreamRDO")
+ ostream.TakeItemsFromInput = True # Copies all data from input file to output
+ ostream.ItemList += ItemList
+ return acc
diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8e60aefaf9db417d6a583b2b2ee59dd279d1f05a
--- /dev/null
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
@@ -0,0 +1,84 @@
+#include "CaloWaveformDigiAlg.h"
+
+#include "Identifier/Identifier.h"
+
+#include <vector>
+#include <map>
+
+CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name,
+ ISvcLocator* pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator) {
+
+}
+
+StatusCode
+CaloWaveformDigiAlg::initialize() {
+ ATH_MSG_INFO(name() << "::initalize()" );
+
+ // Initalize tools
+ ATH_CHECK( m_digiTool.retrieve() );
+
+
+ // Set key to read waveform from
+ ATH_CHECK( m_caloHitContainerKey.initialize() );
+
+ // Set key to write container
+ ATH_CHECK( m_waveformContainerKey.initialize() );
+
+ // Will eventually depend on the type of detector
+ // TODO: Vary time at which centre it?
+ // 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->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);
+
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+CaloWaveformDigiAlg::finalize() {
+ ATH_MSG_INFO(name() << "::finalize()");
+
+ delete m_kernel;
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
+ ATH_MSG_DEBUG("Executing");
+
+ ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number()
+ << " Event: " << ctx.eventID().event_number());
+
+ // Find the input HIT collection
+ SG::ReadHandle<CaloHitCollection> caloHitHandle(m_caloHitContainerKey, ctx);
+
+ ATH_CHECK( caloHitHandle.isValid() );
+ ATH_MSG_DEBUG("Found ReadHandle for CaloHitCollection " << m_caloHitContainerKey);
+
+ // Find the output waveform container
+ SG::WriteHandle<RawWaveformContainer> waveformContainerHandle(m_waveformContainerKey, ctx);
+ ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) );
+
+ ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized");
+
+ if (caloHitHandle->size() == 0) {
+ ATH_MSG_DEBUG("CaloHitCollection found with zero length!");
+ return StatusCode::SUCCESS;
+ }
+
+ // Digitise the hits
+ CHECK( m_digiTool->digitise<CaloHitCollection>(caloHitHandle.ptr(),
+ waveformContainerHandle.ptr(), m_kernel) );
+
+ ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
+
+ return StatusCode::SUCCESS;
+}
diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..9de68257e540890d29143df6a054578334f12503
--- /dev/null
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
@@ -0,0 +1,92 @@
+#ifndef CALODIGIALGS_CALOWAVEFORMDIGIALG_H
+#define CALODIGIALGS_CALOWAVEFORMDIGIALG_H
+
+// Base class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Data classes
+#include "WaveRawEvent/RawWaveformContainer.h"
+#include "FaserCaloSimEvent/CaloHitCollection.h"
+
+// Tool classes
+#include "WaveDigiTools/IWaveformDigitisationTool.h"
+
+// Handles
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+
+// Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+// ROOT
+#include "TF1.h"
+
+// STL
+#include <string>
+
+class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
+
+ public:
+ // Constructor
+ CaloWaveformDigiAlg(const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~CaloWaveformDigiAlg() = default;
+
+ /** @name Usual algorithm methods */
+ //@{
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext& ctx) const override;
+ virtual StatusCode finalize() override;
+ //@}
+
+ private:
+
+ /** @name Disallow default instantiation, copy, assignment */
+ //@{
+ CaloWaveformDigiAlg() = delete;
+ CaloWaveformDigiAlg(const CaloWaveformDigiAlg&) = delete;
+ CaloWaveformDigiAlg &operator=(const CaloWaveformDigiAlg&) = delete;
+ //@}
+
+ Gaudi::Property<double> m_CB_alpha {this, "CB_alpha", 0, "Alpha of the crystal ball function"};
+ 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"};
+
+
+
+
+
+ /// Kernel PDF
+ TF1* m_kernel;
+
+
+ /**
+ * @name Digitisation tool
+ */
+ ToolHandle<IWaveformDigitisationTool> m_digiTool
+ {this, "WaveformDigitisationTool", "WaveformDigitisationTool"};
+
+
+ /**
+ * @name Input HITS using SG::ReadHandleKey
+ */
+ //@{
+
+ SG::ReadHandleKey<CaloHitCollection> m_caloHitContainerKey
+ {this, "CaloHitContainerKey", ""};
+
+ //@}
+
+
+ /**
+ * @name Output data using SG::WriteHandleKey
+ */
+ //@{
+ SG::WriteHandleKey<RawWaveformContainer> m_waveformContainerKey
+ {this, "WaveformContainerKey", ""};
+ //@}
+
+};
+
+#endif // CALODIGIALGS_CALODIGIALG_H
diff --git a/Calorimeter/CaloDigiAlgs/src/components/CaloDigiAlgs_entries.cxx b/Calorimeter/CaloDigiAlgs/src/components/CaloDigiAlgs_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..832cfbb096ffdb13430ba92dc2bfc7b1125bd53c
--- /dev/null
+++ b/Calorimeter/CaloDigiAlgs/src/components/CaloDigiAlgs_entries.cxx
@@ -0,0 +1,2 @@
+#include "../CaloWaveformDigiAlg.h"
+DECLARE_COMPONENT( CaloWaveformDigiAlg )
diff --git a/Control/CalypsoExample/BarcodeChecker/CMakeLists.txt b/Control/CalypsoExample/BarcodeChecker/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8f1e720b867398acbdbdcba483abcf61ed84fefc
--- /dev/null
+++ b/Control/CalypsoExample/BarcodeChecker/CMakeLists.txt
@@ -0,0 +1,11 @@
+################################################################################
+# Package: BarcodeChecker
+################################################################################
+
+# Declare the package name:
+atlas_subdir( BarcodeChecker )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
+
+atlas_install_joboptions( share/*.py )
\ No newline at end of file
diff --git a/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py b/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py
new file mode 100644
index 0000000000000000000000000000000000000000..cc92bd8ac574dbaccfe522366666462b6219d988
--- /dev/null
+++ b/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py
@@ -0,0 +1,38 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+import AthenaPython.PyAthena as PyAthena
+from AthenaPython.PyAthena import StatusCode, McEventCollection, HepMC, CLHEP
+import McParticleEvent.Pythonizations
+
+__author__ = "Dave Caser <dcasper@uci.edu>"
+
+class BarcodeCheckerAlg(PyAthena.Alg):
+ def __init__(self, name="BarcodeChecker", MCEventKey="TruthEvent"):
+ super(BarcodeCheckerAlg,self).__init__(name=name)
+ self.MCEventKey = MCEventKey
+ return
+
+ def initialize(self):
+ self.maxLow = 0
+ self.maxMid = 0
+ self.maxHi = 0
+ return StatusCode.Success
+
+
+ def execute(self):
+ evtCollection = self.evtStore[self.MCEventKey]
+ for mcEvt in evtCollection:
+ for mcParticle in mcEvt.particles:
+ barCode = mcParticle.barcode()
+ self.maxLow = max(self.maxLow, barCode%200000)
+ if barCode%1000000 > 200000:
+ self.maxMid = max(self.maxMid, barCode%1000000 - 200000)
+ self.maxHi = max(self.maxHi, barCode//1000000)
+ return StatusCode.Success
+
+ def finalize(self):
+ print("Low part: ", self.maxLow, " out of 200000 (",100*self.maxLow/200000,"% of overflow)")
+ print("Mid part: ", self.maxMid, " out of ", 1000000 - 200000, " (",100*self.maxMid/(1000000-200000),"% of overflow")
+ print("Hi part: ", self.maxHi, " out of ", (1<<31)//1000000, " (", 100*self.maxHi/((1<<31)//1000000),"% of overflow")
+
+ return StatusCode.Success
\ No newline at end of file
diff --git a/Control/CalypsoExample/BarcodeChecker/python/__init__.py b/Control/CalypsoExample/BarcodeChecker/python/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..d13ae164caa4e93bf2899cea4e67d0ec515784a2
--- /dev/null
+++ b/Control/CalypsoExample/BarcodeChecker/python/__init__.py
@@ -0,0 +1 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
diff --git a/Control/CalypsoExample/BarcodeChecker/share/BarcodeChecker_jobOptions.py b/Control/CalypsoExample/BarcodeChecker/share/BarcodeChecker_jobOptions.py
new file mode 100644
index 0000000000000000000000000000000000000000..1e739beb0c3865e924721b21fe5683b9f18fdaec
--- /dev/null
+++ b/Control/CalypsoExample/BarcodeChecker/share/BarcodeChecker_jobOptions.py
@@ -0,0 +1,29 @@
+#
+# Usage: athena.py -c'INPUT=["faser.1.gfaser.root","faser.2.gfaser.root"]' BarcodeChecker/BarcodeChecker_jobOptions.py >& BarcodeChecker.log
+#
+# INPUT is mandatory (INPUT can be a list, as shown above)
+#
+
+if not 'INPUT' in dir():
+ print("Missing INPUT parameter")
+ exit()
+
+if not isinstance(INPUT, (list,tuple)):
+ INPUT = [INPUT]
+ pass
+
+from AthenaCommon.GlobalFlags import globalflags
+
+globalflags.InputFormat.set_Value_and_Lock('pool')
+
+import AthenaPoolCnvSvc.ReadAthenaPool
+
+svcMgr.EventSelector.InputCollections = INPUT
+
+from BarcodeChecker.BarcodeCheckerAlg import BarcodeCheckerAlg
+
+from AthenaCommon.AlgSequence import AlgSequence
+job = AlgSequence()
+job += BarcodeCheckerAlg()
+
+theApp.EvtMax = -1
diff --git a/Control/CalypsoExample/Digitization/CMakeLists.txt b/Control/CalypsoExample/Digitization/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f8caebb598f6cea848edafea449a0bd19e657c6d
--- /dev/null
+++ b/Control/CalypsoExample/Digitization/CMakeLists.txt
@@ -0,0 +1,26 @@
+################################################################################
+# Package: Digitization
+################################################################################
+
+# Declare the package name:
+atlas_subdir( Digitization )
+
+# Component(s) in the package:
+#atlas_add_component( GeoModelTest
+# src/GeoModelTestAlg.cxx
+# src/components/GeoModelTest_entries.cxx
+# INCLUDE_DIRS ${GEOMODEL_INCLUDE_DIRS}
+# LINK_LIBRARIES ${GEOMODEL_LIBRARIES} AthenaBaseComps GeoModelFaserUtilities ScintReadoutGeometry TrackerReadoutGeometry MagFieldInterfaces MagFieldElements MagFieldConditions )
+
+# Install files from the package:
+#atlas_install_joboptions( share/*.py )
+#atlas_install_python_modules( python/*.py )
+atlas_install_scripts( scripts/*.sh scripts/*.py )
+
+#atlas_add_test( ProdRecoTI12
+# SCRIPT scripts/faser_reco.py ${CMAKE_CURRENT_SOURCE_DIR}/../rawdata/Faser-Physics-001920-filtered.raw TI12Data
+# PROPERTIES TIMEOUT 300 )
+#atlas_add_test( ProdRecoTestBeam
+# SCRIPT scripts/faser_reco.py ${CMAKE_CURRENT_SOURCE_DIR}/../RAWDATA/Faser-Physics-003613-filtered.raw TestBeamData
+# PROPERTIES TIMEOUT 300 )
+
diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi.py b/Control/CalypsoExample/Digitization/scripts/faser_digi.py
new file mode 100755
index 0000000000000000000000000000000000000000..5ef818999634d592b21142b88ad2cee5c1a324f1
--- /dev/null
+++ b/Control/CalypsoExample/Digitization/scripts/faser_digi.py
@@ -0,0 +1,207 @@
+#!/usr/bin/env python
+#
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+# Run with:
+# ./faser_digi.py filepath runtype
+#
+# filepath - fully qualified path, including url if needed, to the input HITS file
+# example: "root://eospublic.cern.ch//eos/experiment/faser/sim/GeniePilot/HITS/1/faser.150fbInv.1.001.HITS.pool.root"
+#
+# runtype - MANDATORY flag to specify the data type (TI12OldMC or TI12MC or TestBeamMC).
+# Not extracted (yet) from file path for MC data
+#
+import sys
+import argparse
+
+parser = argparse.ArgumentParser(description="Run FASER reconstruction")
+
+parser.add_argument("file_path",
+ help="Fully qualified path of the raw input file")
+parser.add_argument("run_type", nargs="?", default="",
+ help="Specify run type (if it can't be parsed from path)")
+parser.add_argument("-r", "--reco", default="",
+ help="Specify reco tag (to append to output filename)")
+parser.add_argument("-n", "--nevents", type=int, default=-1,
+ help="Specify number of events to process (default: all)")
+parser.add_argument("-v", "--verbose", action='store_true',
+ help="Turn on DEBUG output")
+
+args = parser.parse_args()
+
+from pathlib import Path
+
+filepath=Path(args.file_path)
+
+# runtype has been provided
+if len(args.run_type) > 0:
+ runtype=args.run_type
+
+# Extract runtype from path
+# Should be directory above run
+# i.e.: TestBeamData/Run-004150/Faser-Physics-004150-00000.raw"
+else:
+ if True or len(filepath.parts) < 3:
+ print("Can't determine run type from path - specify on command line ")
+ sys.exit(-1)
+
+# runtype = filepath.parts[-3]
+
+print(f"Starting digitization of {filepath.name} with type {runtype}")
+if args.nevents > 0:
+ print(f"Reconstructing {args.nevents} events by command-line option")
+
+# Start digitization
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from AthenaCommon.Constants import VERBOSE, INFO
+
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+
+Configurable.configurableRun3Behavior = True
+
+# Flags for this job
+ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
+ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
+
+ConfigFlags.Input.ProjectName = "mc20"
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+ConfigFlags.Digitization.TruthOutput = True
+
+# TI12 old geometry
+if runtype == "TI12OldMC":
+ ConfigFlags.GeoModel.FaserVersion = "FASER-01"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+
+# Testbeam setup
+elif runtype == "TestBeamMC" :
+ ConfigFlags.GeoModel.FaserVersion = "FASER-TB00"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-TB00"
+
+# New TI12 geometry (ugh)
+elif runtype == "TI12MC":
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-02"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02"
+
+else:
+ print("Invalid run type found:", runtype)
+ print("Specify correct type or update list")
+ sys.exit(-1)
+
+
+# Must use original input string here, as pathlib mangles double // in path names
+ConfigFlags.Input.Files = [ args.file_path ]
+
+filestem = filepath.stem
+if len(args.reco) > 0:
+ filestem += f"-{args.reco}"
+
+# ConfigFlags.addFlag("Output.xAODFileName", f"{filestem}-xAOD.root")
+ConfigFlags.Output.RDOFileName = f"{filestem}-RDO.root"
+
+#
+# Play around with this?
+# ConfigFlags.Concurrency.NumThreads = 2
+# ConfigFlags.Concurrency.NumConcurrentEvents = 2
+ConfigFlags.lock()
+
+#
+# Configure components
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+acc.merge(PoolWriteCfg(ConfigFlags))
+
+#
+# Needed, or move to MainServicesCfg?
+from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+acc.merge(FaserGeometryCfg(ConfigFlags))
+
+# Set up algorithms
+from FaserSCT_Digitization.FaserSCT_DigitizationConfigNew import FaserSCT_DigitizationCfg
+acc.merge(FaserSCT_DigitizationCfg(ConfigFlags))
+
+from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
+acc.merge(CaloWaveformDigitizationCfg(ConfigFlags))
+
+from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg
+acc.merge(ScintWaveformDigitizationCfg(ConfigFlags))
+
+# from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
+# acc.merge(WaveformReconstructionCfg(ConfigFlags))
+
+# # Not ready for primetime
+# # from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg
+# # acc.merge(CalorimeterReconstructionCfg(ConfigFlags))
+
+# # Tracker clusters
+# from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+# acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+
+# # SpacePoints
+# from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+# acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+
+# # Try Dave's fitter
+# from TrackerClusterFit.TrackerClusterFitConfig import ClusterFitAlgCfg
+# acc.merge(ClusterFitAlgCfg(ConfigFlags))
+
+#
+# Configure output
+# from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+# itemList = [ "xAOD::EventInfo#*"
+# , "xAOD::EventAuxInfo#*"
+# , "xAOD::FaserTriggerData#*"
+# , "xAOD::FaserTriggerDataAux#*"
+# , "FaserSCT_RDO_Container#*"
+# # , "Tracker::FaserSCT_ClusterContainer#*"
+# # , "FaserSCT_SpacePointContainer#*"
+# # , "FaserSCT_SpacePointOverlapCollection#*"
+# # , "TrackCollection#*"
+# ]
+# acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList))
+
+# Waveform reconstruction output
+# from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg
+# acc.merge(WaveformReconstructionOutputCfg(ConfigFlags))
+
+# Calorimeter reconstruction output
+# from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg
+# acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags))
+
+# Check what we have
+# 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
+
+# Configure verbosity
+# ConfigFlags.dump()
+if args.verbose:
+ acc.foreach_component("*").OutputLevel = VERBOSE
+
+ #acc.getService("FaserByteStreamInputSvc").DumpFlag = True
+ #acc.getService("FaserEventSelector").OutputLevel = VERBOSE
+ #acc.getService("FaserByteStreamInputSvc").OutputLevel = VERBOSE
+ #acc.getService("FaserByteStreamCnvSvc").OutputLevel = VERBOSE
+ #acc.getService("FaserByteStreamAddressProviderSvc").OutputLevel = VERBOSE
+
+else:
+ acc.foreach_component("*").OutputLevel = INFO
+
+acc.foreach_component("*ClassID*").OutputLevel = INFO
+
+acc.getService("MessageSvc").Format = "% F%40W%S%7W%R%T %0W%M"
+
+# Execute and finish
+sys.exit(int(acc.run(maxEvents=args.nevents).isFailure()))
diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 7c1f662f4cfdad38026a90a7c01c82e834a0ef01..f0fa35b4aaa071c3a0e699e74f8a949a8f7765c9 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -7,9 +7,11 @@
# filepath - fully qualified path, including url if needed, to the input raw data file
# example: "root://hepatl30//atlas/local/torrence/faser/commissioning/TestBeamData/Run-004150/Faser-Physics-004150-00000.raw"
#
-# runtype - optional flag to specify the data type (TI12Data or TestBeamData).
+# runtype - optionally specify the data type (TI12Data, TI12Data02, or TestBeamData).
# In a normal file system location, this will be extracted from the directory name,
# but runtype will override this assignment.
+# TI12Data02 is needed for the IFT geometry. Script will auto-detect this if read
+# from normal file system location.
#
import sys
import argparse
@@ -26,6 +28,8 @@ parser.add_argument("-n", "--nevents", type=int, default=-1,
help="Specify number of events to process (default: all)")
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)")
args = parser.parse_args()
@@ -47,6 +51,23 @@ else:
runtype = filepath.parts[-3]
+ # Fix TI12 geometry versions as well (needed in production)
+ # Probably better to do this from configuration in upstream production scripts,
+ # so lets call this a hack for now
+ if runtype == "TI12Data":
+
+ runname = filepath.parts[-2]
+ try:
+ runnumber = int(runname.split('-')[1])
+ except Exception as e:
+ print(f"Failed to find run number in {filepath}")
+ print(f"Couldn't parse {runname}")
+ print(f"Leave runtype as {runtype}!")
+ else:
+ if runnumber > 5302: # Last TI12 run on Nov. 23, 2021 without IFT
+ print(f"Found run number {runnumber}, using TI12 configuration with IFT")
+ runtype = "TI12Data02"
+
print(f"Starting reconstruction of {filepath.name} with type {runtype}")
if args.nevents > 0:
print(f"Reconstructing {args.nevents} events by command-line option")
@@ -75,10 +96,15 @@ if runtype == "TI12Data":
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
# Testbeam setup
-elif runtype == "TestBeamData":
+elif runtype == "TestBeamData" or runtype == "TestBeam2021":
ConfigFlags.GeoModel.FaserVersion = "FASER-TB00"
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-TB00"
+# New TI12 geometry (ugh)
+elif runtype == "TI12Data02":
+ ConfigFlags.GeoModel.FaserVersion = "FASER-02"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02"
+
else:
print("Invalid run type found:", runtype)
print("Specify correct type or update list")
@@ -135,9 +161,18 @@ acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
-# Try Dave's fitter
-from TrackerClusterFit.TrackerClusterFitConfig import ClusterFitAlgCfg
-acc.merge(ClusterFitAlgCfg(ConfigFlags))
+# Can't use both in the same job, as they write to the same output histograms
+if args.clusterFit:
+ print("Configuring TrackerClusterFit (old)")
+ # Try Dave's fitter
+ from TrackerClusterFit.TrackerClusterFitConfig import ClusterFitAlgCfg
+ acc.merge(ClusterFitAlgCfg(ConfigFlags))
+
+else:
+ print("Configuring TrackerSegmentFit (new)")
+ # Try Dave's new fitter
+ from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+ acc.merge(SegmentFitAlgCfg(ConfigFlags))
#
# Configure output
diff --git a/Control/CalypsoExample/TruthEventDumper/CMakeLists.txt b/Control/CalypsoExample/TruthEventDumper/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ceea9c78f5a8d088b1d17d9f3c6af2d325758f9c
--- /dev/null
+++ b/Control/CalypsoExample/TruthEventDumper/CMakeLists.txt
@@ -0,0 +1,11 @@
+################################################################################
+# Package: TruthEventDumper
+################################################################################
+
+# Declare the package name:
+atlas_subdir( TruthEventDumper )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
+
+atlas_install_joboptions( share/*.py )
\ No newline at end of file
diff --git a/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py b/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py
new file mode 100644
index 0000000000000000000000000000000000000000..934433e2f375d6f0883a5a66e66dba70f1515427
--- /dev/null
+++ b/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py
@@ -0,0 +1,39 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+import AthenaPython.PyAthena as PyAthena
+from AthenaPython.PyAthena import StatusCode, McEventCollection, HepMC, CLHEP
+import McParticleEvent.Pythonizations
+
+__author__ = "Dave Caser <dcasper@uci.edu>"
+
+class TruthEventDumperAlg(PyAthena.Alg):
+ def __init__(self, name="TruthEventDumper", MCEventKey="TruthEvent"):
+ super(TruthEventDumperAlg,self).__init__(name=name)
+ self.MCEventKey = MCEventKey
+ return
+
+ def initialize(self):
+ # self.maxLow = 0
+ # self.maxMid = 0
+ # self.maxHi = 0
+ return StatusCode.Success
+
+
+ def execute(self):
+ evtCollection = self.evtStore[self.MCEventKey]
+ for mcEvt in evtCollection:
+ mcEvt.print()
+ # for mcParticle in mcEvt.particles:
+ # barCode = mcParticle.barcode()
+ # self.maxLow = max(self.maxLow, barCode%200000)
+ # if barCode%1000000 > 200000:
+ # self.maxMid = max(self.maxMid, barCode%1000000 - 200000)
+ # self.maxHi = max(self.maxHi, barCode//1000000)
+ return StatusCode.Success
+
+ def finalize(self):
+ # print("Low part: ", self.maxLow, " out of 200000 (",100*self.maxLow/200000,"% of overflow)")
+ # print("Mid part: ", self.maxMid, " out of ", 1000000 - 200000, " (",100*self.maxMid/(1000000-200000),"% of overflow")
+ # print("Hi part: ", self.maxHi, " out of ", (1<<31)//1000000, " (", 100*self.maxHi/((1<<31)//1000000),"% of overflow")
+
+ return StatusCode.Success
\ No newline at end of file
diff --git a/Control/CalypsoExample/TruthEventDumper/python/__init__.py b/Control/CalypsoExample/TruthEventDumper/python/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..d13ae164caa4e93bf2899cea4e67d0ec515784a2
--- /dev/null
+++ b/Control/CalypsoExample/TruthEventDumper/python/__init__.py
@@ -0,0 +1 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
diff --git a/Control/CalypsoExample/TruthEventDumper/share/TruthEventDumper_jobOptions.py b/Control/CalypsoExample/TruthEventDumper/share/TruthEventDumper_jobOptions.py
new file mode 100644
index 0000000000000000000000000000000000000000..303f4b8d26db18b0cf96f30e500548ecc1a84208
--- /dev/null
+++ b/Control/CalypsoExample/TruthEventDumper/share/TruthEventDumper_jobOptions.py
@@ -0,0 +1,41 @@
+#
+# Usage: athena.py -c'INPUT=["faser.1.gfaser.root","faser.2.gfaser.root"]; COLLECTION="TruthEvent"; MAXEVT=-1; SKIPEVT=0' TruthEventDumper/TruthEventDumper_jobOptions.py >& TruthEventDumper.log
+#
+# INPUT is mandatory (INPUT can be a list, as shown above)
+# COLLECTION would normally be either "TruthEvent" (Geant4 particles) or "BeamTruthEvent" (generator particles)
+# MAXEVT and SKIPEVT are self-explanatory and optional
+#
+
+if not 'INPUT' in dir():
+ print("Missing INPUT parameter")
+ exit()
+
+if not 'MAXEVT' in dir():
+ MAXEVT = -1
+
+if not 'SKIPEVT' in dir():
+ SKIPEVT = 0
+
+if not 'COLLECTION' in dir():
+ COLLECTION = "TruthEvent"
+
+if not isinstance(INPUT, (list,tuple)):
+ INPUT = [INPUT]
+ pass
+
+from AthenaCommon.GlobalFlags import globalflags
+
+globalflags.InputFormat.set_Value_and_Lock('pool')
+
+import AthenaPoolCnvSvc.ReadAthenaPool
+
+svcMgr.EventSelector.InputCollections = INPUT
+
+from TruthEventDumper.TruthEventDumperAlg import TruthEventDumperAlg
+
+from AthenaCommon.AlgSequence import AlgSequence
+job = AlgSequence()
+job += TruthEventDumperAlg(MCEventKey=COLLECTION)
+
+theApp.EvtMax = MAXEVT
+theApp.SkipEvents = SKIPEVT
diff --git a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
index f5ba68bcba67791b3d7fa6792518ba74b90d1bb6..28bb6147a3276b510f3b4c3de83ad281a2997af4 100644
--- a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
+++ b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
@@ -12,6 +12,10 @@ from AthenaConfiguration.ComponentFactory import CompFactory
from AthenaCommon.SystemOfUnits import TeV
from AthenaCommon.PhysicalConstants import pi
+### add radial pos sampler ### with gaussian beam implemented
+from FaserParticleGun.RadialPosSampler import RadialPosSampler
+
+
def FaserParticleGunCommonCfg(ConfigFlags, **kwargs) :
cfg = ComponentAccumulator(AthSequencer("AthBeginSeq", Sequential = True))
@@ -34,6 +38,7 @@ def FaserParticleGunSingleParticleCfg(ConfigFlags, **kwargs) :
pg = cfg.getPrimary()
+ pg.sampler.n = kwargs.setdefault("n", 1)
pg.sampler.pid = kwargs.setdefault("pid", -13)
pg.sampler.mom = PG.EThetaMPhiSampler(energy = kwargs.setdefault("energy", 1*TeV),
theta = kwargs.setdefault("theta", [0, pi/20]),
@@ -46,6 +51,32 @@ def FaserParticleGunSingleParticleCfg(ConfigFlags, **kwargs) :
return cfg
+def FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs) :
+ cfg = FaserParticleGunCommonCfg(ConfigFlags, **kwargs)
+
+ pg = cfg.getPrimary()
+
+ pg.sampler.pid = kwargs.setdefault("pid", 13)
+ pg.sampler.mom = PG.EThetaMPhiSampler(energy = kwargs.setdefault("energy", 1*TeV),
+ theta = kwargs.setdefault("theta", 0.0),
+ phi = kwargs.setdefault("phi", 0.0),
+ mass = kwargs.setdefault("mass", 105.7) )
+
+ if "radius" in kwargs:
+ pg.sampler.pos = RadialPosSampler(x = kwargs.setdefault("x", 0.0),
+ y = kwargs.setdefault("y", 0.0),
+ z = kwargs.setdefault("z", 0.0),
+ r = kwargs.setdefault("radius", 1.0),
+ t = kwargs.setdefault("t", 0.0) )
+ else:
+ pg.sampler.pos = PG.PosSampler(x = kwargs.setdefault("x", 0.0),
+ y = kwargs.setdefault("y", 0.0),
+ z = kwargs.setdefault("z", 0.0),
+ t = kwargs.setdefault("t", 0.0) )
+
+ return cfg
+
+'''
def FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs) :
cfg = FaserParticleGunCommonCfg(ConfigFlags, **kwargs)
@@ -62,7 +93,7 @@ def FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs) :
t = kwargs.setdefault("t", 0.0) )
return cfg
-
+'''
def FaserParticleGunCosmicsCfg(ConfigFlags, **kwargs) :
# Supported keyword arguments:
#
diff --git a/Generators/GenieReader/python/GenieReaderAlg.py b/Generators/GenieReader/python/GenieReaderAlg.py
index fc1cb184f0016fe4f73599a9f1eb1edab209f5f1..c97ea494b089a25ff88e4c2a1f2abdd4bf472c6b 100644
--- a/Generators/GenieReader/python/GenieReaderAlg.py
+++ b/Generators/GenieReader/python/GenieReaderAlg.py
@@ -2,7 +2,8 @@
from AthenaCommon.AppMgr import ServiceMgr as svcMgr
from GeneratorModules.EvgenAlg import EvgenAlg
-from AthenaPython.PyAthena import StatusCode
+from AthenaPython.PyAthena import StatusCode, EventInfo, EventID, EventType
+from AthenaCommon.SystemOfUnits import GeV, m
import ROOT
__author__ = "Dave Caser <dcasper@uci.edu>"
@@ -20,7 +21,23 @@ class GenieReader(EvgenAlg):
from AthenaPython.PyAthena import HepMC as HepMC
evt.weights().push_back(1.0)
- pos = HepMC.FourVector(self.evtStore["vx"]*1000, self.evtStore["vy"]*1000, self.evtStore["vz"]*1000, 0)
+ treeEventInfo = self.evtStore["TTreeEventInfo"]
+ treeEventId = treeEventInfo.event_ID()
+ runNumber = treeEventId.run_number()
+ eventNumber = treeEventId.event_number()
+ # print("found run/event number ", runNumber, "/", eventNumber)
+
+ mcEventType = EventType()
+ mcEventType.add_type(EventType.IS_SIMULATION)
+
+ mcEventId = EventID(run_number = runNumber, event_number = eventNumber)
+ mcEventInfo = EventInfo(id = mcEventId, type = mcEventType)
+ self.evtStore.record(mcEventInfo, "McEventInfo", True, False)
+ ROOT.SetOwnership(mcEventType, False)
+ ROOT.SetOwnership(mcEventId, False)
+ ROOT.SetOwnership(mcEventInfo, False)
+
+ pos = HepMC.FourVector(self.evtStore["vx"]*m, self.evtStore["vy"]*m, self.evtStore["vz"]*m, 0)
gv = HepMC.GenVertex(pos)
ROOT.SetOwnership(gv, False)
evt.add_vertex(gv)
@@ -32,13 +49,13 @@ class GenieReader(EvgenAlg):
py = list(self.evtStore["py"])
pz = list(self.evtStore["pz"])
E = list(self.evtStore["E"])
- M = list(self.evtStore["M"])
+ M = list(self.evtStore["M"])
for i in range(nParticles):
gp = HepMC.GenParticle()
- mom = HepMC.FourVector(px[i], py[i], pz[i], E[i])
+ mom = HepMC.FourVector(px[i]*GeV, py[i]*GeV, pz[i]*GeV, E[i]*GeV)
gp.set_momentum(mom)
- gp.set_generated_mass(M[i])
+ gp.set_generated_mass(M[i]*GeV)
gp.set_pdg_id(pdgc[i])
genie_status = status[i]
if (genie_status == 0): # initial particle
diff --git a/Generators/GenieReader/share/GenieReader_jobOptions.py b/Generators/GenieReader/share/GenieReader_jobOptions.py
index 91ad7c1318b772fd8cb015eda03ea274b0a6815a..649e705b66eb05f7b5ac2b331890a914794bec60 100644
--- a/Generators/GenieReader/share/GenieReader_jobOptions.py
+++ b/Generators/GenieReader/share/GenieReader_jobOptions.py
@@ -30,9 +30,10 @@ job += GenieReader()
from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream
ostream = AthenaPoolOutputStream( "StreamEVGEN" , OUTPUT, noTag=True )
-ostream.ItemList += [ "xAOD::EventInfo#*", "xAOD::EventAuxInfo#*",
- "EventInfo#*",
+ostream.ItemList.remove("EventInfo#*")
+ostream.ItemList += [ "EventInfo#McEventInfo",
"McEventCollection#*" ]
+print(ostream.ItemList)
if not 'EVTMAX' in dir():
EVTMAX = -1
diff --git a/Scintillator/ScintDigiAlgs/CMakeLists.txt b/Scintillator/ScintDigiAlgs/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..43695c5f90d96053a8b44fa3ac41c0f4c159ad99
--- /dev/null
+++ b/Scintillator/ScintDigiAlgs/CMakeLists.txt
@@ -0,0 +1,15 @@
+################################################################################
+# Package: ScintDigiAlgs
+################################################################################
+
+# Declare the package name:
+atlas_subdir( ScintDigiAlgs )
+
+# Component(s) in the package:
+atlas_add_component( ScintDigiAlgs
+ src/*.cxx src/*.h
+ src/components/*.cxx
+ LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent ScintSimEvent WaveDigiToolsLib)
+
+atlas_install_python_modules( python/*.py )
+
diff --git a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..243e4493dd117e6758c894a7904115b3d3941ec3
--- /dev/null
+++ b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
@@ -0,0 +1,66 @@
+# Copyright (C) 2020-2021 CERN for the benefit of the FASER collaboration
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+# Crystallball parameter dictionary used in simulated digitized wave reconstruction.
+# Crystalball function Parameters estimated from Deion's slides uploaded at
+# 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)
+
+# One stop shopping for normal FASER data
+def ScintWaveformDigitizationCfg(flags):
+ """ Return all algorithms and tools for Waveform digitization """
+ acc = ComponentAccumulator()
+
+ if not flags.Input.isMC:
+ return acc
+
+ if "TB" not in flags.GeoModel.FaserVersion:
+ acc.merge(ScintWaveformDigiCfg(flags, "TimingWaveformDigiAlg", "Trigger"))
+ acc.merge(ScintWaveformDigiCfg(flags, "VetoWaveformDigiAlg", "Veto"))
+ acc.merge(ScintWaveformDigiCfg(flags, "PreshowerWaveformDigiAlg", "Preshower"))
+ acc.merge(ScintWaveformDigitizationOutputCfg(flags))
+ return acc
+
+# Return configured digitization algorithm from SIM hits
+# Specify data source (Veto, Trigger, Preshower)
+def ScintWaveformDigiCfg(flags, name="ScintWaveformDigiAlg", source="", **kwargs):
+
+ acc = ComponentAccumulator()
+
+ tool = CompFactory.WaveformDigitisationTool(name=source+"WaveformDigtisationTool", **kwargs)
+
+ kwargs.setdefault("ScintHitContainerKey", source+"Hits")
+ kwargs.setdefault("WaveformContainerKey", source+"Waveforms")
+
+ digiAlg = CompFactory.ScintWaveformDigiAlg(name, **kwargs)
+ digiAlg.CB_alpha = dict_CB_param[source]["CB_alpha"]
+ 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"]
+
+ kwargs.setdefault("WaveformDigitisationTool", tool)
+
+ acc.addEventAlgo(digiAlg)
+
+ return acc
+
+def ScintWaveformDigitizationOutputCfg(flags, **kwargs):
+ """ Return ComponentAccumulator with output for Waveform Digi"""
+ acc = ComponentAccumulator()
+ ItemList = [
+ "RawWaveformContainer#*"
+ ]
+ acc.merge(OutputStreamCfg(flags, "RDO"))
+ ostream = acc.getEventAlgo("OutputStreamRDO")
+ ostream.TakeItemsFromInput = True # Copies all data from input file to output
+ ostream.ItemList += ItemList
+ return acc
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b25080419ad0a7ea64f3933a0eb3acfdfe4c7eac
--- /dev/null
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -0,0 +1,85 @@
+
+#include "ScintWaveformDigiAlg.h"
+
+#include "Identifier/Identifier.h"
+
+#include <vector>
+#include <map>
+
+
+ScintWaveformDigiAlg::ScintWaveformDigiAlg(const std::string& name,
+ ISvcLocator* pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator) {
+
+}
+
+StatusCode
+ScintWaveformDigiAlg::initialize() {
+ ATH_MSG_INFO(name() << "::initalize()" );
+
+ // Initalize tools
+ ATH_CHECK( m_digiTool.retrieve() );
+
+
+ // Set key to read waveform from
+ ATH_CHECK( m_scintHitContainerKey.initialize() );
+
+ // Set key to write container
+ ATH_CHECK( m_waveformContainerKey.initialize() );
+
+ // Will eventually depend on the type of detector
+ // TODO: Vary time at which centre it?
+ // TODO: Better parameters
+
+
+
+ m_kernel = new TF1("PDF", "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);
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+ScintWaveformDigiAlg::finalize() {
+ ATH_MSG_INFO(name() << "::finalize()");
+
+ delete m_kernel;
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
+ ATH_MSG_DEBUG("Executing");
+
+ ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number()
+ << " Event: " << ctx.eventID().event_number());
+
+ // Find the input HIT collection
+ SG::ReadHandle<ScintHitCollection> scintHitHandle(m_scintHitContainerKey, ctx);
+
+ ATH_CHECK( scintHitHandle.isValid() );
+ ATH_MSG_DEBUG("Found ReadHandle for ScintHitCollection " << m_scintHitContainerKey);
+
+ // Find the output waveform container
+ SG::WriteHandle<RawWaveformContainer> waveformContainerHandle(m_waveformContainerKey, ctx);
+ ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) );
+
+ ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized");
+
+ if (scintHitHandle->size() == 0) {
+ ATH_MSG_DEBUG("ScintHitCollection found with zero length!");
+ return StatusCode::SUCCESS;
+ }
+
+ // Digitise the hits
+ CHECK( m_digiTool->digitise<ScintHitCollection>(scintHitHandle.ptr(),
+ waveformContainerHandle.ptr(), m_kernel) );
+
+ ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
+
+ return StatusCode::SUCCESS;
+}
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..af60e02d8169194ccc2459f0b7a69727b2544eba
--- /dev/null
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
@@ -0,0 +1,91 @@
+#ifndef SCINTDIGIALGS_SCINTWAVEFORMDIGIALG_H
+#define SCINTDIGIALGS_SCINTWAVEFORMDIGIALG_H
+
+// Base class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Data classes
+#include "WaveRawEvent/RawWaveformContainer.h"
+#include "ScintSimEvent/ScintHitCollection.h"
+
+// Tool classes
+#include "WaveDigiTools/IWaveformDigitisationTool.h"
+
+// Handles
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+
+// Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+// ROOT
+#include "TF1.h"
+
+
+// STL
+#include <string>
+
+class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
+
+ public:
+ // Constructor
+ ScintWaveformDigiAlg(const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~ScintWaveformDigiAlg() = default;
+
+ /** @name Usual algorithm methods */
+ //@{
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext& ctx) const override;
+ virtual StatusCode finalize() override;
+ //@}
+
+
+ private:
+
+ /** @name Disallow default instantiation, copy, assignment */
+ //@{
+ ScintWaveformDigiAlg() = delete;
+ ScintWaveformDigiAlg(const ScintWaveformDigiAlg&) = delete;
+ ScintWaveformDigiAlg &operator=(const ScintWaveformDigiAlg&) = delete;
+ //@}
+
+ Gaudi::Property<double> m_CB_alpha {this, "CB_alpha", 0, "Alpha of the crystal ball function"};
+ 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"};
+
+
+ /// Kernel PDF
+ TF1* m_kernel;
+
+
+ /**
+ * @name Digitisation tool
+ */
+ ToolHandle<IWaveformDigitisationTool> m_digiTool
+ {this, "WaveformDigitisationTool", "WaveformDigitisationTool"};
+
+
+ /**
+ * @name Input HITS using SG::ReadHandleKey
+ */
+ //@{
+
+ SG::ReadHandleKey<ScintHitCollection> m_scintHitContainerKey
+ {this, "ScintHitContainerKey", ""};
+
+ //@}
+
+
+ /**
+ * @name Output data using SG::WriteHandleKey
+ */
+ //@{
+ SG::WriteHandleKey<RawWaveformContainer> m_waveformContainerKey
+ {this, "WaveformContainerKey", ""};
+ //@}
+
+};
+
+#endif // SCINTDIGIALGS_SCINTDIGIALG_H
diff --git a/Scintillator/ScintDigiAlgs/src/components/ScintDigiAlgs_entries.cxx b/Scintillator/ScintDigiAlgs/src/components/ScintDigiAlgs_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..3a8deff1f06b692e1cf0928472fb12426274f059
--- /dev/null
+++ b/Scintillator/ScintDigiAlgs/src/components/ScintDigiAlgs_entries.cxx
@@ -0,0 +1,2 @@
+#include "../ScintWaveformDigiAlg.h"
+DECLARE_COMPONENT( ScintWaveformDigiAlg )
diff --git a/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.cxx b/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.cxx
index 88592d298b2eb974007a73c15ad2413beb193648..f9c937aa22c8a69153e7c32174d4c1009e9b5845 100644
--- a/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.cxx
+++ b/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.cxx
@@ -3,6 +3,7 @@
*/
#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1.h"
+#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h"
#include "ScintSimEventTPCnv/ScintHits/ScintHit_p1.h"
#include "ScintHitCollectionCnv.h"
@@ -18,11 +19,19 @@ ScintHitCollection_PERS* ScintHitCollectionCnv::createPersistent(ScintHitCollect
ScintHitCollection* ScintHitCollectionCnv::createTransient() {
MsgStream mlog(msgSvc(), "ScintHitCollectionConverter" );
ScintHitCollectionCnv_p1 converter_p1;
+ ScintHitCollectionCnv_p1a converter_p1a;
static const pool::Guid p1_guid("B2573A16-4B46-4E1E-98E3-F93421680779");
+ static const pool::Guid p1a_guid("0DFC461F-9FF5-4447-B4F0-7CC7157191D1");
ScintHitCollection *trans_cont(0);
- if( this->compareClassGuid(p1_guid)) {
+ if( this->compareClassGuid(p1a_guid))
+ {
+ std::unique_ptr< ScintHitCollection_p1a > col_vect( this->poolReadObject< ScintHitCollection_p1a >() );
+ trans_cont = converter_p1a.createTransient( col_vect.get(), mlog );
+ }
+ else if( this->compareClassGuid(p1_guid))
+ {
std::unique_ptr< ScintHitCollection_p1 > col_vect( this->poolReadObject< ScintHitCollection_p1 >() );
trans_cont = converter_p1.createTransient( col_vect.get(), mlog );
} else {
diff --git a/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.h b/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.h
index ca9c329ddad032c35766ad6a271fd7637734b1c0..4cdaf2046a595f763c424baae5f6942ab91a33fd 100644
--- a/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.h
+++ b/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.h
@@ -8,12 +8,14 @@
#include "ScintSimEvent/ScintHitCollection.h"
#include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1.h"
#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1.h"
+#include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h"
+#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h"
#include "AthenaPoolCnvSvc/T_AthenaPoolCustomCnv.h"
// Gaudi
#include "GaudiKernel/MsgStream.h"
// typedef to the latest persistent version
-typedef ScintHitCollection_p1 ScintHitCollection_PERS;
-typedef ScintHitCollectionCnv_p1 ScintHitCollectionCnv_PERS;
+typedef ScintHitCollection_p1a ScintHitCollection_PERS;
+typedef ScintHitCollectionCnv_p1a ScintHitCollectionCnv_PERS;
class ScintHitCollectionCnv : public T_AthenaPoolCustomCnv<ScintHitCollection, ScintHitCollection_PERS > {
friend class CnvFactory<ScintHitCollectionCnv>;
diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h
new file mode 100644
index 0000000000000000000000000000000000000000..4f8de2b4c704240402bdf5ea23a36cf3a12acefe
--- /dev/null
+++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h
@@ -0,0 +1,41 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef SCINTHITCOLLECTIONCNV_P1A_H
+#define SCINTHITCOLLECTIONCNV_P1A_H
+
+// ScintHitCollectionCnv_p1, T/P separation of Scint Hits
+// author D.Costanzo <davide.costanzo@cern.ch>
+// author O.Arnaez <olivier.arnaez@cern.ch>
+
+#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h"
+#include "ScintSimEvent/ScintHitCollection.h"
+#include "ScintHitCollection_p1a.h"
+
+
+class ScintHitCollectionCnv_p1a : public T_AthenaPoolTPCnvBase<ScintHitCollection, ScintHitCollection_p1a>
+{
+ public:
+
+ ScintHitCollectionCnv_p1a() {};
+
+ virtual ScintHitCollection* createTransient(const ScintHitCollection_p1a* persObj, MsgStream &log);
+
+ virtual void persToTrans(const ScintHitCollection_p1a* persCont,
+ ScintHitCollection* transCont,
+ MsgStream &log) ;
+ virtual void transToPers(const ScintHitCollection* transCont,
+ ScintHitCollection_p1a* persCont,
+ MsgStream &log) ;
+
+ private:
+
+ static const double m_persEneUnit;
+ static const double m_persLenUnit;
+ static const double m_persAngUnit;
+ static const double m_2bHalfMaximum;
+ static const int m_2bMaximum;
+};
+
+#endif
diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h
new file mode 100644
index 0000000000000000000000000000000000000000..8e714d6e571c746ae6376536ed62c5681476232e
--- /dev/null
+++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h
@@ -0,0 +1,57 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef SCINTHITCOLLECTION_P1A_H
+#define SCINTHITCOLLECTION_P1A_H
+
+/*
+
+Authors: Davide Costanzo Rob Duxfield
+
+*/
+
+#include <vector>
+#include <string>
+
+class ScintHitCollection_p1a
+{
+ public:
+/// Default constructor
+ ScintHitCollection_p1a ();
+ //private:
+
+ std::vector<float> m_hit1_meanTime; // 1 element per string
+ std::vector<float> m_hit1_x0; //
+ std::vector<float> m_hit1_y0; //
+ std::vector<float> m_hit1_z0; //
+ std::vector<float> m_hit1_theta; //
+ std::vector<float> m_hit1_phi; //
+ std::vector<unsigned long> m_nHits; //
+
+ std::vector<unsigned short> m_hitEne_2b; // 1 element per hit
+ std::vector<unsigned short> m_hitLength_2b; //
+
+ std::vector<unsigned short> m_dTheta; // 1 element per hit except for first hit in string
+ std::vector<unsigned short> m_dPhi; //
+
+ std::vector<float> m_hitEne_4b; // 1 element per hit with m_hitEne_2b[i] == 2**16
+
+ std::vector<float> m_hitLength_4b; // 1 element per hit with m_hitLength_2b[i] == 2**16
+
+ std::vector<unsigned long> m_barcode;
+ std::vector<unsigned long> m_mcEvtIndex;
+ std::vector<char> m_evtColl;
+ std::vector<unsigned long> m_nBC;
+
+ std::vector<unsigned long> m_id;
+ std::vector<unsigned long> m_nId;
+};
+
+
+// inlines
+
+inline
+ScintHitCollection_p1a::ScintHitCollection_p1a () {}
+
+#endif
diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintSimEventTPCnvDict.h b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintSimEventTPCnvDict.h
index f0042a97ec671d78b8fb783af5f516116570952a..1606077dd0efcd14371c1c7445d7a553b0b0579d 100644
--- a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintSimEventTPCnvDict.h
+++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintSimEventTPCnvDict.h
@@ -15,6 +15,8 @@
#include "ScintSimEventTPCnv/ScintHits/ScintHitCnv_p1.h"
#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1.h"
#include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1.h"
+#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h"
+#include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h"
#include "ScintSimEventTPCnv/ScintHits/ScintHit_p1.h"
#endif // SCINTEVENTTPCNV_SCINTSIMEVENTTPCNVDICT_H
diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/selection.xml b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/selection.xml
index 6b0d915eac9a1e2fedb5f8254bbbf419dec16dac..345d08165aeb11723e94f19e7195f42fcb5a88af 100644
--- a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/selection.xml
+++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/selection.xml
@@ -4,4 +4,5 @@
<class name="ScintHit_p1" />
<class name="std::vector<ScintHit_p1>" />
<class name="ScintHitCollection_p1" id="B2573A16-4B46-4E1E-98E3-F93421680779" />
+ <class name="ScintHitCollection_p1a" id="0DFC461F-9FF5-4447-B4F0-7CC7157191D1" />
</lcgdict>
diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/src/ScintHits/ScintHitCollectionCnv_p1a.cxx b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/src/ScintHits/ScintHitCollectionCnv_p1a.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..cd162d1ad3c8064966a59c3517b19b233bc29c80
--- /dev/null
+++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/src/ScintHits/ScintHitCollectionCnv_p1a.cxx
@@ -0,0 +1,334 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "ScintSimEvent/ScintHit.h"
+#include "ScintSimEvent/ScintHitCollection.h"
+#include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h"
+#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h"
+#include "GeneratorObjects/HepMcParticleLink.h"
+
+#include <cmath>
+
+//CLHEP
+#include "CLHEP/Geometry/Point3D.h"
+// Gaudi
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/ThreadLocalContext.h"
+
+// Athena
+#include "AthenaKernel/ExtendedEventContext.h"
+#include "StoreGate/StoreGateSvc.h"
+
+// * * * stolen from eflowRec * * * //
+inline double phicorr(double a)
+{
+ if (a <= -M_PI)
+ {
+ return a+(2*M_PI*floor(-(a-M_PI)/(2*M_PI)));
+ }
+ else if (a > M_PI)
+ {
+ return a-(2*M_PI*floor((a+M_PI)/(2*M_PI)));
+ }
+ else
+ {
+ return a;
+ }
+}
+
+// * * * stolen from eflowRec * * * //
+inline double cycle(double a, double b)
+{
+ double del = b-a;
+ if (del > M_PI)
+ {
+ return a+2.0*M_PI;
+ }
+ else if (del < -M_PI)
+ {
+ return a-2.0*M_PI;
+ }
+ else
+ {
+ return a;
+ }
+}
+
+
+const double ScintHitCollectionCnv_p1a::m_persEneUnit = 1.0e-5;
+const double ScintHitCollectionCnv_p1a::m_persLenUnit = 1.0e-5;
+const double ScintHitCollectionCnv_p1a::m_persAngUnit = 1.0e-5;
+const double ScintHitCollectionCnv_p1a::m_2bHalfMaximum = pow(2.0, 15.0);
+const int ScintHitCollectionCnv_p1a::m_2bMaximum = (unsigned short)(-1);
+
+
+void ScintHitCollectionCnv_p1a::transToPers(const ScintHitCollection* transCont, ScintHitCollection_p1a* persCont, MsgStream &log)
+{
+ // Finds hits belonging to a "string" (in which the end point of one hit is the same as the start point of the next) and
+ // persistifies the end point of each hit plus the start point of the first hit in each string.
+ //
+ // Further compression is achieved by optimising the storage of the position vectors:- start (x,y,z) and (theta,phi) of
+ // first hit are stored as floats, (delta_theta,delta_phi) relative to the fisrst hit are stored as 2 byte numbers and
+ // used to specify the hit direction. All hit lengths are stored as 2 byte numbers.
+ //
+ // Additional savings are achieved by storing the energy loss for each hit as a 2 byte number and only storing the mean
+ // time of the first hit per string.
+ //
+ // See http://indico.cern.ch/getFile.py/access?contribId=11&resId=2&materialId=slides&confId=30893 for more info.
+
+ static const double dRcut = 1.0e-7;
+ static const double dTcut = 1.0;
+
+ const EventContext& ctx = Gaudi::Hive::currentContext();
+ const IProxyDict* proxy = Atlas::getExtendedEventContext(ctx).proxy();
+ const HepMcParticleLink * lastLink=nullptr;
+ int lastId = -1;
+ double stringFirstTheta = 0.0;
+ double stringFirstPhi = 0.0;
+ double lastT = 0.0;
+ double persSumE = 0.0;
+ double transSumE = 0.0;
+ unsigned int idx = 0;
+ unsigned int endBC = 0;
+ unsigned int endId = 0;
+ unsigned int endHit = 0;
+ HepGeom::Point3D<double> lastTransEnd(0.0, 0.0, 0.0);
+ HepGeom::Point3D<double> lastPersEnd(0.0, 0.0, 0.0);
+
+ for (ScintHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) {
+
+ ScintHitCollection::const_iterator siHit = it;
+
+
+ if ( !lastLink || (siHit->particleLink() != *lastLink) ) {
+
+ // store barcode, eventIndex and McEventCollection once for set of consecutive hits with same barcode
+
+ lastLink = &(siHit->particleLink());
+ persCont->m_barcode.push_back(lastLink->barcode());
+ unsigned short index{0};
+ const HepMcParticleLink::index_type position =
+ HepMcParticleLink::getEventPositionInCollection(lastLink->eventIndex(),
+ lastLink->getEventCollection(),
+ proxy).at(0);
+ if (position!=0) {
+ index = lastLink->eventIndex();
+ if(lastLink->eventIndex()!=static_cast<HepMcParticleLink::index_type>(index)) {
+ log << MSG::WARNING << "Attempting to persistify an eventIndex larger than max unsigned short!" << endmsg;
+ }
+ }
+ persCont->m_mcEvtIndex.push_back(index);
+ persCont->m_evtColl.push_back(lastLink->getEventCollectionAsChar());
+
+ if (idx > 0) {
+ persCont->m_nBC.push_back(idx - endBC);
+ endBC = idx;
+ }
+ }
+
+ if ( (int)siHit->identify() != lastId ) {
+
+ // store id once for set of consecutive hits with same barcode
+
+ lastId = siHit->identify();
+ persCont->m_id.push_back(lastId);
+
+ if (idx > 0) {
+ persCont->m_nId.push_back(idx - endId);
+ endId = idx;
+ }
+ }
+
+ HepGeom::Point3D<double> st = siHit->localStartPosition();
+ HepGeom::Point3D<double> en = siHit->localEndPosition();
+
+ const double dx = st.x() - lastTransEnd.x();
+ const double dy = st.y() - lastTransEnd.y();
+ const double dz = st.z() - lastTransEnd.z();
+ const double t = siHit->meanTime();
+
+ const double dRLast = sqrt(dx * dx + dy * dy + dz * dz); // dR between end of previous hit and start of current one
+ const double dTLast = fabs(t - lastT);
+
+ CLHEP::Hep3Vector direction(0.0, 0.0, 0.0);
+ double theta = 0.0;
+ double phi = 0.0;
+ bool startNewString = false;
+
+ if (dRLast < dRcut && dTLast < dTcut) {
+
+ // hit is part of existing string
+
+ direction = CLHEP::Hep3Vector( en.x() - lastPersEnd.x(), en.y() - lastPersEnd.y(), en.z() - lastPersEnd.z() );
+
+ theta = direction.theta();
+ phi = phicorr( direction.phi() );
+
+ const int dTheta_2b = (int)( (theta - stringFirstTheta) / m_persAngUnit + m_2bHalfMaximum + 0.5 );
+ const int dPhi_2b = (int)( (cycle(phi, stringFirstPhi) - stringFirstPhi) / m_persAngUnit + m_2bHalfMaximum + 0.5 );
+
+ if ( dTheta_2b < m_2bMaximum && dTheta_2b >= 0 && dPhi_2b < m_2bMaximum && dPhi_2b >= 0) {
+ persCont->m_dTheta.push_back(dTheta_2b);
+ persCont->m_dPhi.push_back(dPhi_2b);
+ theta = stringFirstTheta + ( (double)dTheta_2b - m_2bHalfMaximum ) * m_persAngUnit;
+ phi = stringFirstPhi + ( (double)dPhi_2b - m_2bHalfMaximum ) * m_persAngUnit;
+ phi = phicorr(phi);
+ }
+ else {
+ startNewString = true;
+ }
+ }
+
+ if (startNewString || dRLast >= dRcut || dTLast >= dTcut) {
+
+ // begin new hit string
+
+ direction = CLHEP::Hep3Vector( en.x() - st.x(), en.y() - st.y(), en.z() - st.z() );
+
+ theta = direction.theta();
+ phi = phicorr( direction.phi() );
+
+ persCont->m_hit1_meanTime.push_back(t);
+ persCont->m_hit1_x0.push_back(st.x());
+ persCont->m_hit1_y0.push_back(st.y());
+ persCont->m_hit1_z0.push_back(st.z());
+ persCont->m_hit1_theta.push_back(theta);
+ persCont->m_hit1_phi.push_back(phi);
+
+ lastPersEnd = HepGeom::Point3D<double>(st.x(), st.y(), st.z());
+
+ stringFirstTheta = theta;
+ stringFirstPhi = phi;
+
+ if (idx > 0) {
+ persCont->m_nHits.push_back(idx - endHit);
+ endHit = idx;
+ }
+ }
+
+ lastTransEnd = HepGeom::Point3D<double>(en.x(), en.y(), en.z());
+ transSumE += siHit->energyLoss();
+
+ const int eneLoss_2b = (int)((transSumE - persSumE) / m_persEneUnit + 0.5); // calculated to allow recovery sum over
+ // whole hit string to chosen precision
+
+ const int hitLength_2b = (int)(direction.mag() / m_persLenUnit + 0.5); // calculated to give the correct position to
+ // the chosen precision, NOT the length of the
+ // hit (small difference in practice).
+ double eneLoss = 0.0;
+
+ if (eneLoss_2b >= m_2bMaximum) {
+ eneLoss = siHit->energyLoss();
+ persCont->m_hitEne_2b.push_back(m_2bMaximum);
+ persCont->m_hitEne_4b.push_back(eneLoss);
+ }
+ else {
+ eneLoss = eneLoss_2b * m_persEneUnit;
+ persCont->m_hitEne_2b.push_back(eneLoss_2b);
+ }
+
+ double length = 0.0;
+
+ if (hitLength_2b >= m_2bMaximum) {
+ length = direction.mag();
+ persCont->m_hitLength_2b.push_back(m_2bMaximum);
+ persCont->m_hitLength_4b.push_back(direction.mag());
+ }
+ else {
+ length = hitLength_2b * m_persLenUnit;
+ persCont->m_hitLength_2b.push_back(hitLength_2b);
+ }
+
+ CLHEP::Hep3Vector persDir(length, 0.0, 0.0);
+ persDir.setTheta(theta);
+ persDir.setPhi(phi);
+
+ lastPersEnd = (CLHEP::Hep3Vector)lastPersEnd + persDir;
+ persSumE += eneLoss;
+ lastT = t;
+
+ ++idx;
+ }
+
+ persCont->m_nBC.push_back(idx - endBC);
+ persCont->m_nId.push_back(idx - endId);
+ persCont->m_nHits.push_back(idx - endHit);
+}
+
+
+ScintHitCollection* ScintHitCollectionCnv_p1a::createTransient(const ScintHitCollection_p1a* persObj, MsgStream &log) {
+ std::unique_ptr<ScintHitCollection> trans(std::make_unique<ScintHitCollection>("DefaultCollectionName",persObj->m_nHits.size()));
+ persToTrans(persObj, trans.get(), log);
+ return(trans.release());
+}
+
+
+void ScintHitCollectionCnv_p1a::persToTrans(const ScintHitCollection_p1a* persCont, ScintHitCollection* transCont, MsgStream &/*log*/)
+{
+ const EventContext& ctx = Gaudi::Hive::currentContext();
+
+ unsigned int hitCount = 0;
+ unsigned int angleCount = 0;
+ unsigned int idxBC = 0;
+ unsigned int idxId = 0;
+ unsigned int idxEne4b = 0;
+ unsigned int idxLen4b = 0;
+ unsigned int endHit = 0;
+ unsigned int endBC = 0;
+ unsigned int endId = 0;
+
+ for (unsigned int i = 0; i < persCont->m_nHits.size(); i++) {
+
+ if (persCont->m_nHits[i]) {
+
+ const unsigned int start = endHit;
+ endHit += persCont->m_nHits[i];
+
+ const double t0 = persCont->m_hit1_meanTime[i];
+ const double theta0 = persCont->m_hit1_theta[i];
+ const double phi0 = persCont->m_hit1_phi[i];
+ HepGeom::Point3D<double> endLast(persCont->m_hit1_x0[i], persCont->m_hit1_y0[i], persCont->m_hit1_z0[i]);
+
+ for (unsigned int j = start; j < endHit; j++) {
+
+ if (j >= endBC + persCont->m_nBC[idxBC])
+ endBC += persCont->m_nBC[idxBC++];
+
+ if (j >= endId + persCont->m_nId[idxId])
+ endId += persCont->m_nId[idxId++];
+
+ const double eneLoss_2b = persCont->m_hitEne_2b[hitCount];
+ const double hitLength_2b = persCont->m_hitLength_2b[hitCount];
+
+ const double eneLoss = (eneLoss_2b < m_2bMaximum) ? eneLoss_2b * m_persEneUnit : persCont->m_hitEne_4b[idxEne4b++];
+ const double length = (hitLength_2b < m_2bMaximum) ? hitLength_2b * m_persLenUnit : persCont->m_hitLength_4b[idxLen4b++];
+
+ const double dTheta = (j > start) ? ((double)persCont->m_dTheta[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0;
+ const double dPhi = (j > start) ? ((double)persCont->m_dPhi[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0;
+
+ const double meanTime = t0;
+ const double theta = theta0 + dTheta;
+ const double phi = phicorr(phi0 + dPhi);
+
+ CLHEP::Hep3Vector r(length, 0.0, 0.0);
+ r.setTheta(theta);
+ r.setPhi(phi);
+
+ HepGeom::Point3D<double> endThis( endLast + r );
+
+ HepMcParticleLink::PositionFlag flag = HepMcParticleLink::IS_INDEX;
+ if (persCont->m_mcEvtIndex[idxBC] == 0) {
+ flag = HepMcParticleLink::IS_POSITION;
+ }
+ HepMcParticleLink partLink( persCont->m_barcode[idxBC], persCont->m_mcEvtIndex[idxBC], HepMcParticleLink::ExtendedBarCode::eventCollectionFromChar(persCont->m_evtColl[idxBC]), flag, ctx );
+ transCont->Emplace( endLast, endThis, eneLoss, meanTime, partLink, persCont->m_id[idxId]);
+
+ endLast = endThis;
+
+ ++hitCount;
+ if (j > start) ++angleCount;
+ }
+ }
+ }
+}
diff --git a/Simulation/ISF/ISF_Core/FaserISF_Services/python/FaserISF_ServicesConfigNew.py b/Simulation/ISF/ISF_Core/FaserISF_Services/python/FaserISF_ServicesConfigNew.py
index 7260d59fc1e6df1b327196817f8a44c865d72a16..e507545cd536d1a85889e6768330a77c5755c5c5 100644
--- a/Simulation/ISF/ISF_Core/FaserISF_Services/python/FaserISF_ServicesConfigNew.py
+++ b/Simulation/ISF/ISF_Core/FaserISF_Services/python/FaserISF_ServicesConfigNew.py
@@ -10,7 +10,8 @@ from AthenaConfiguration.ComponentFactory import CompFactory
from BarcodeServices.BarcodeServicesConfigNew import MC15aPlusBarcodeSvcCfg
from ISF_HepMC_Tools.ISF_HepMC_ToolsConfigNew import ParticleFinalStateFilterCfg, GenParticleInteractingFilterCfg
-from FaserISF_HepMC_Tools.FaserISF_HepMC_ToolsConfigNew import FaserTruthStrategyCfg, FaserDipoleTruthStrategyCfg
+# from FaserISF_HepMC_Tools.FaserISF_HepMC_ToolsConfigNew import FaserTruthStrategyCfg, FaserDipoleTruthStrategyCfg
+from FaserISF_HepMC_Tools.FaserISF_HepMC_ToolsConfigNew import TruthStrategyGroupCfg
ISF__FaserTruthSvc, ISF__FaserGeoIDSvc, ISF__FaserInputConverter = CompFactory.getComps("ISF::FaserTruthSvc","ISF::FaserGeoIDSvc","ISF::FaserInputConverter")
@@ -60,10 +61,12 @@ def FaserTruthServiceCfg(ConfigFlags, name="FaserISF_TruthService", **kwargs):
result = MC15aPlusBarcodeSvcCfg(ConfigFlags)
kwargs.setdefault('BarcodeSvc', result.getService("Barcode_MC15aPlusBarcodeSvc") )
- acc = FaserTruthStrategyCfg(ConfigFlags)
- acc2= FaserDipoleTruthStrategyCfg(ConfigFlags)
- kwargs.setdefault('TruthStrategies',[result.popToolsAndMerge(acc), result.popToolsAndMerge(acc2)])
-
+ # acc = FaserTruthStrategyCfg(ConfigFlags)
+ # acc2= FaserDipoleTruthStrategyCfg(ConfigFlags)
+ # kwargs.setdefault('TruthStrategies',[result.popToolsAndMerge(acc), result.popToolsAndMerge(acc2)])
+ acc = TruthStrategyGroupCfg(ConfigFlags)
+ kwargs.setdefault('TruthStrategies', [result.popToolsAndMerge(acc)])
+
kwargs.setdefault('SkipIfNoChildren', True)
kwargs.setdefault('SkipIfNoParentBarcode', True)
diff --git a/Simulation/ISF/ISF_Core/FaserISF_Services/src/FaserTruthSvc.cxx b/Simulation/ISF/ISF_Core/FaserISF_Services/src/FaserTruthSvc.cxx
index b1a64408271187de69318c448fb727056c66f704..c43929542e093404fe9f53c884b59196f510e8d5 100644
--- a/Simulation/ISF/ISF_Core/FaserISF_Services/src/FaserTruthSvc.cxx
+++ b/Simulation/ISF/ISF_Core/FaserISF_Services/src/FaserTruthSvc.cxx
@@ -128,6 +128,8 @@ StatusCode ISF::FaserTruthSvc::finalize()
StatusCode ISF::FaserTruthSvc::initializeTruthCollection()
{
ATH_CHECK( m_barcodeSvc->initializeBarcodes() );
+ m_incrementProcs.clear();
+ m_secondaryProcs.clear();
return StatusCode::SUCCESS;
}
@@ -264,13 +266,20 @@ void ISF::FaserTruthSvc::recordIncidentToMCTruth( ISF::IFaserTruthIncident& ti)
newPrimBC = this->maxGeneratedParticleBarcode(ti.parentParticle()->parent_event())+1;
}
else {
+ if (m_incrementProcs.count(processCode) == 0)
+ m_incrementProcs[processCode] = 0;
+ m_incrementProcs[processCode]++;
newPrimBC = m_barcodeSvc->incrementBarcode( parentBC, processCode);
}
if ( newPrimBC == Barcode::fUndefinedBarcode) {
+ for (auto pair : m_incrementProcs)
+ {
+ ATH_MSG_ALWAYS("Process code: " << pair.first << " : " << pair.second);
+ }
if (m_ignoreUndefinedBarcodes) {
ATH_MSG_WARNING("Unable to generate new Particle Barcode. Continuing due to 'IgnoreUndefinedBarcodes'==True");
} else {
- ATH_MSG_FATAL("Unable to generate new Particle Barcode in region " << ti.geoID() << ", at vertex (" << vtx->position().x() << ", " << vtx->position().y() << ", " << vtx->position().z() << "). Aborting");
+ ATH_MSG_FATAL("Unable to generate new Particle Barcode from " << parentBC << "( id= " << ti.parentPdgCode() << ", Ekin = " << ti.parentEkin()/Gaudi::Units::GeV << " GeV) with process code " << processCode << " in region " << ti.geoID() << ", at vertex (" << vtx->position().x() << ", " << vtx->position().y() << ", " << vtx->position().z() << "). Aborting");
abort();
}
}
@@ -370,12 +379,20 @@ void ISF::FaserTruthSvc::recordIncidentToMCTruth( ISF::IFaserTruthIncident& ti)
}
else {
// generate a new barcode for the child particle
+ if (m_secondaryProcs.count(processCode) == 0)
+ m_secondaryProcs[processCode] = 0;
+ m_secondaryProcs[processCode]++;
Barcode::ParticleBarcode secBC = m_barcodeSvc->newSecondary( parentBC, processCode);
if ( secBC == Barcode::fUndefinedBarcode) {
+ for (auto pair : m_secondaryProcs)
+ {
+ ATH_MSG_ALWAYS("Process code: " << pair.first << " : " << pair.second);
+ }
+
if (m_ignoreUndefinedBarcodes)
ATH_MSG_WARNING("Unable to generate new Secondary Particle Barcode. Continuing due to 'IgnoreUndefinedBarcodes'==True");
else {
- ATH_MSG_ERROR("Unable to generate new Secondary Particle Barcode. Aborting");
+ ATH_MSG_FATAL("Unable to generate new Secondary Particle Barcode from " << parentBC << " with process code " << processCode <<" in region " << ti.geoID() << ", at vertex (" << vtx->position().x() << ", " << vtx->position().y() << ", " << vtx->position().z() << "). Aborting");
abort();
}
}
diff --git a/Simulation/ISF/ISF_Core/FaserISF_Services/src/FaserTruthSvc.h b/Simulation/ISF/ISF_Core/FaserISF_Services/src/FaserTruthSvc.h
index d674d80782cd3cbc5dfbef9e27939a1f4923dea8..21b4d2717eab17f96cf6e5a40cd25e45499c2b8e 100644
--- a/Simulation/ISF/ISF_Core/FaserISF_Services/src/FaserTruthSvc.h
+++ b/Simulation/ISF/ISF_Core/FaserISF_Services/src/FaserTruthSvc.h
@@ -126,6 +126,10 @@ namespace ISF {
bool m_quasiStableParticlesIncluded; //!< does this job simulate quasi-stable particles.
+ mutable std::map<int, int> m_incrementProcs;
+ mutable std::map<int, int> m_secondaryProcs;
+
+
};
}
diff --git a/Simulation/ISF/ISF_HepMC/FaserISF_HepMC_Tools/python/FaserISF_HepMC_ToolsConfigNew.py b/Simulation/ISF/ISF_HepMC/FaserISF_HepMC_Tools/python/FaserISF_HepMC_ToolsConfigNew.py
index 0f251682c98dfba7a08fd7e9a7427fb6a330e91b..4f739f6c0ca89afab6020eab6f42e096e1f53545 100644
--- a/Simulation/ISF/ISF_HepMC/FaserISF_HepMC_Tools/python/FaserISF_HepMC_ToolsConfigNew.py
+++ b/Simulation/ISF/ISF_HepMC/FaserISF_HepMC_Tools/python/FaserISF_HepMC_ToolsConfigNew.py
@@ -153,44 +153,65 @@ def FaserParticleGenericFilterCfg(ConfigFlags, name="ISF_FaserGenericFilter", **
# http://www-geant4.kek.jp/lxr/source//processes/hadronic/management/include/G4HadronicProcessType.hh#L46
-def FaserDipoleTruthStrategyCfg(ConfigFlags, name="ISF_FaserDipoleTruthStrategy", **kwargs):
- result = ComponentAccumulator()
+def TruthStrategyGroupCfg(ConfigFlags, name="ISF_MCTruthStrategyGroupID", **kwargs):
import ROOT, cppyy
cppyy.load_library('FaserDetDescrDict')
FaserRegion = ROOT.FaserDetDescr.FaserRegion
- #
- # Save truth in Dipole region
- #
- kwargs.setdefault('Regions', [FaserRegion.fFaserDipole,
- FaserRegion.fFaserNeutrino,
+
+ result = ComponentAccumulator()
+ kwargs.setdefault("ParentMinEkin", 100.*MeV)
+ kwargs.setdefault("ChildMinEkin" , 100.*MeV)
+ kwargs.setdefault("VertexTypes", [3, 14, 15, 4, 5, 6, 7, 2, 12, 13])
+ kwargs.setdefault("VertexTypeRangeLow" , 201) # All kinds of decay processes
+ kwargs.setdefault("VertexTypeRangeHigh" , 298) # ...
+ kwargs.setdefault("Regions", [FaserRegion.fFaserNeutrino,
+ FaserRegion.fFaserScintillator,
+ FaserRegion.fFaserTracker,
+ FaserRegion.fFaserCalorimeter,
FaserRegion.fFaserCavern])
- kwargs.setdefault('ParentMinEkin', 1000.0*MeV)
- kwargs.setdefault('ChildMinEkin', 1000.0*MeV)
- result.setPrivateTools(ISF__FaserTruthStrategy(name, **kwargs))
+ result.setPrivateTools(CompFactory.ISF.FaserTruthStrategy(name, **kwargs))
return result
-def FaserTruthStrategyCfg(ConfigFlags, name="ISF_FaserTruthStrategy", **kwargs):
- result = ComponentAccumulator()
+# def FaserDipoleTruthStrategyCfg(ConfigFlags, name="ISF_FaserDipoleTruthStrategy", **kwargs):
+# result = ComponentAccumulator()
- import ROOT, cppyy
- cppyy.load_library('FaserDetDescrDict')
- FaserRegion = ROOT.FaserDetDescr.FaserRegion
- #
- # Save truth in all regions except Dipole
- #
- kwargs.setdefault('Regions', [
- # FaserRegion.fFaserNeutrino,
- FaserRegion.fFaserScintillator,
- FaserRegion.fFaserTracker])
- # FaserRegion.fFaserDipole,
- # FaserRegion.fFaserCalorimeter,
- # FaserRegion.fFaserCavern])
- # kwargs.setdefault('ParentMinEkin', 0.1*MeV)
- # kwargs.setdefault('ChildMinEkin', 0.1*MeV)
- result.setPrivateTools(ISF__FaserTruthStrategy(name, **kwargs))
- return result
+# import ROOT, cppyy
+# cppyy.load_library('FaserDetDescrDict')
+# FaserRegion = ROOT.FaserDetDescr.FaserRegion
+# #
+# # Save truth in Dipole region
+# #
+# kwargs.setdefault('Regions', [FaserRegion.fFaserDipole,
+# FaserRegion.fFaserNeutrino,
+# FaserRegion.fFaserCavern])
+# kwargs.setdefault('ParentMinEkin', 1000.0*MeV)
+# kwargs.setdefault('ChildMinEkin', 1000.0*MeV)
+# result.setPrivateTools(ISF__FaserTruthStrategy(name, **kwargs))
+# return result
+
+
+# def FaserTruthStrategyCfg(ConfigFlags, name="ISF_FaserTruthStrategy", **kwargs):
+# result = ComponentAccumulator()
+
+# import ROOT, cppyy
+# cppyy.load_library('FaserDetDescrDict')
+# FaserRegion = ROOT.FaserDetDescr.FaserRegion
+# #
+# # Save truth in all regions except Dipole
+# #
+# kwargs.setdefault('Regions', [
+# # FaserRegion.fFaserNeutrino,
+# FaserRegion.fFaserScintillator,
+# FaserRegion.fFaserTracker])
+# # FaserRegion.fFaserDipole,
+# # FaserRegion.fFaserCalorimeter,
+# # FaserRegion.fFaserCavern])
+# # kwargs.setdefault('ParentMinEkin', 0.1*MeV)
+# # kwargs.setdefault('ChildMinEkin', 0.1*MeV)
+# result.setPrivateTools(ISF__FaserTruthStrategy(name, **kwargs))
+# return result
# def TruthStrategyGroupID_MC15Cfg(ConfigFlags, name="ISF_MCTruthStrategyGroupID_MC15", **kwargs):
# result = ComponentAccumulator()
diff --git a/Tracker/TrackerDetDescr/TrackerReadoutGeometry/TrackerReadoutGeometry/SiDetectorElement.h b/Tracker/TrackerDetDescr/TrackerReadoutGeometry/TrackerReadoutGeometry/SiDetectorElement.h
index 652add49be92b54bcf1b273d8d1da11f1ad4fd52..3286b43b3f8b66d91ee8fe864cd3218a8301802c 100644
--- a/Tracker/TrackerDetDescr/TrackerReadoutGeometry/TrackerReadoutGeometry/SiDetectorElement.h
+++ b/Tracker/TrackerDetDescr/TrackerReadoutGeometry/TrackerReadoutGeometry/SiDetectorElement.h
@@ -221,6 +221,8 @@ namespace TrackerDD {
// Position
/// Local (simulation/hit frame) to global transform
virtual const GeoTrf::Transform3D& transformHit() const;
+ /// Local (simulation/hit frame) to module transform
+ virtual const GeoTrf::Transform3D& transformModule() const;
/// Local (reconstruction frame) to global transform
virtual const Amg::Transform3D& transform() const override final;
/**
@@ -340,40 +342,37 @@ namespace TrackerDD {
@name Module Frame
Methods to help work with the module frame.
- This is mainly of of use in the SCT where the module frame can
- in general be different from the element frame. It is actully
- defined as the frame of one of the sides (currently the axial
- side), but using these methods one does not need to make any
- assumptions about what the actual frame is. In the following
- methods the local element frame is the local reconstruction
- frame of this element.
+ This is mainly of of use in the SCT where the module frame is
+ different from the element frame. We have removed the ATLAS convention
+ in which the module frame was one of the wafers.
*/
//@{
+
/// Module to global frame transform.
/// Includes misalignment. The module frame is defined to be the
- /// local reconstruction frame of the axial layer in the SCT. For
- /// the pixel it is the same as the element.
+ /// frame in which the wafers are centered, rotated by +/- stereo angle, and equally offset in depth.
+ /// For the pixel it is the same as the element.
//const HepGeom::Transform3D & moduleTransform() const;
+ const HepGeom::Transform3D& moduleTransformCLHEP() const;
const Amg::Transform3D& moduleTransform() const;
/// Default module to global frame transform, ie with no misalignment.
/// The module frame is defined to be the
- /// local reconstruction frame of the axial layer in the SCT. For
- /// the pixel it is the same as the element.
+ /// frame in which each wafer is centered, rotated by +/- stereo angle, and offset in depth.
+ /// For the pixel it is the same as the element.
Amg::Transform3D defModuleTransform() const;
+ /// Default to global transform
+ /// ie with no misalignment.
+ const HepGeom::Transform3D defModuleTransformCLHEP() const;
/// Take a transform of the local element frame and return its equivalent in the module frame.
//HepGeom::Transform3D localToModuleFrame(const HepGeom::Transform3D & localTransform) const;
Amg::Transform3D localToModuleFrame(const Amg::Transform3D& localTransform) const;
- /// Transformation from local element to module frame. This can be
- /// used to take a local position in the element frame and transform
- /// it to a position in the module frame. If one is already in the
- /// module frame it will return the Identity transform.
- //HepGeom::Transform3D localToModuleTransform() const;
+ /// Transformation from local element to module frame.
Amg::Transform3D localToModuleTransform() const;
/// Check if the element and module frame are the same.
@@ -695,6 +694,8 @@ namespace TrackerDD {
mutable Amg::Transform3D m_transform ATLAS_THREAD_SAFE; // Guarded by m_mutex
mutable HepGeom::Transform3D m_transformCLHEP ATLAS_THREAD_SAFE; // Guarded by m_mutex
+ mutable Amg::Transform3D m_moduleTransform ATLAS_THREAD_SAFE; // Guarded by m_mutex
+ mutable HepGeom::Transform3D m_moduleTransformCLHEP ATLAS_THREAD_SAFE; // Guarded by m_mutex
mutable Amg::Vector3D m_normal ATLAS_THREAD_SAFE; // Guarded by m_mutex
mutable Amg::Vector3D m_etaAxis ATLAS_THREAD_SAFE; // Guarded by m_mutex
diff --git a/Tracker/TrackerDetDescr/TrackerReadoutGeometry/src/SiDetectorElement.cxx b/Tracker/TrackerDetDescr/TrackerReadoutGeometry/src/SiDetectorElement.cxx
index c9a29684d559305ac76ad06b2b632701df8cdc3c..07a1fadae9b0560bd722ed281f0a01b83d514477 100644
--- a/Tracker/TrackerDetDescr/TrackerReadoutGeometry/src/SiDetectorElement.cxx
+++ b/Tracker/TrackerDetDescr/TrackerReadoutGeometry/src/SiDetectorElement.cxx
@@ -122,6 +122,9 @@ namespace TrackerDD {
SiDetectorElement::updateCache() const
{
const GeoTrf::Transform3D& geoTransform = transformHit();
+ const GeoTrf::Transform3D& geoModuleTransform = transformModule();
+ // ATH_MSG_ALWAYS("Wafer transform: " << geoTransform.translation() << "//" << geoTransform.rotation() );
+ // ATH_MSG_ALWAYS("Module transform: " << geoModuleTransform.translation() << "//" << geoModuleTransform.rotation() );
double radialShift = 0.;
@@ -215,6 +218,9 @@ namespace TrackerDD {
m_transformCLHEP = Amg::EigenTransformToCLHEP(geoTransform) * recoToHitTransformImpl();
m_transform = Amg::CLHEPTransformToEigen(m_transformCLHEP);
+ m_moduleTransformCLHEP = Amg::EigenTransformToCLHEP(geoModuleTransform) * recoToHitTransformImpl();
+ m_moduleTransform = Amg::CLHEPTransformToEigen(m_moduleTransformCLHEP);
+
// Check that local frame is right-handed. (ie transform has no reflection)
// This can be done by checking that the determinant is >0.
if (m_firstTime) { // Only need to check this once.
@@ -286,6 +292,28 @@ namespace TrackerDD {
return getMaterialGeom()->getAbsoluteTransform();
}
+ const GeoTrf::Transform3D&
+ SiDetectorElement::transformModule() const
+ {
+ PVConstLink parent = getMaterialGeom()->getParent()->getParent();
+ const GeoVFullPhysVol* fullParent = dynamic_cast<const GeoVFullPhysVol*>(&(*parent));
+ if (fullParent == nullptr)
+ {
+ ATH_MSG_FATAL("Unable to reach parent module volume");
+ if (m_geoAlignStore) {
+ const GeoTrf::Transform3D* ptrXf = m_geoAlignStore->getAbsPosition(getMaterialGeom());
+ if (ptrXf) return *ptrXf;
+ }
+ return getMaterialGeom()->getAbsoluteTransform();
+ }
+ // ATH_MSG_ALWAYS("Found full parent named: " << fullParent->getLogVol()->getName() << " with id " << fullParent->getId());
+ if (m_geoAlignStore) {
+ const GeoTrf::Transform3D* ptrXf = m_geoAlignStore->getAbsPosition(fullParent);
+ if (ptrXf) return *ptrXf;
+ }
+ return fullParent->getAbsoluteTransform();
+ }
+
const Amg::Transform3D&
SiDetectorElement::transform() const
{
@@ -377,40 +405,81 @@ namespace TrackerDD {
const Amg::Transform3D&
SiDetectorElement::moduleTransform() const
{
- return (isModuleFrame()) ? transform() : m_otherSide->transform();
+ if (!m_cacheValid) {
+ std::lock_guard<std::mutex> lock(m_mutex);
+ if (!m_cacheValid) updateCache();
+ }
+
+ return m_moduleTransform;
+ }
+
+ const HepGeom::Transform3D&
+ SiDetectorElement::moduleTransformCLHEP() const
+ {
+ //stuff to get the CLHEP version of the local to global transform.
+ if (!m_cacheValid) {
+ std::lock_guard<std::mutex> lock(m_mutex);
+ if (!m_cacheValid) updateCache();
+ }
+
+ return m_moduleTransformCLHEP;
}
Amg::Transform3D
SiDetectorElement::defModuleTransform() const
{
- return (isModuleFrame()) ? defTransform() : m_otherSide->defTransform();
+ HepGeom::Transform3D tmpTransform = defModuleTransformCLHEP();
+ return Amg::CLHEPTransformToEigen(tmpTransform);
}
+ const HepGeom::Transform3D
+ SiDetectorElement::defModuleTransformCLHEP() const
+ {
+ PVConstLink parent = getMaterialGeom()->getParent()->getParent();
+ const GeoVFullPhysVol* fullParent = dynamic_cast<const GeoVFullPhysVol*>(&(*parent));
+ if (fullParent == nullptr)
+ {
+ ATH_MSG_FATAL("Unable to reach parent module volume");
+ if (m_geoAlignStore) {
+ const GeoTrf::Transform3D* ptrXf = m_geoAlignStore->getDefAbsPosition(getMaterialGeom());
+ if (ptrXf) return Amg::EigenTransformToCLHEP(*ptrXf) * recoToHitTransform();
+ }
+ return Amg::EigenTransformToCLHEP(getMaterialGeom()->getDefAbsoluteTransform()) * recoToHitTransform();
+ }
+ // ATH_MSG_ALWAYS("Found full parent named: " << fullParent->getLogVol()->getName() << " with id " << fullParent->getId());
+ if (m_geoAlignStore) {
+ const GeoTrf::Transform3D* ptrXf = m_geoAlignStore->getDefAbsPosition(fullParent);
+ if (ptrXf) return Amg::EigenTransformToCLHEP(*ptrXf) * recoToHitTransform();
+ }
+ return Amg::EigenTransformToCLHEP(fullParent->getDefAbsoluteTransform()) * recoToHitTransform();
+ }
+
+
// Take a transform in the local reconstruction and return it in the module frame
// For a given transform l in frame A. The equivalent transform in frame B is
// B.inverse() * A * l * A.inverse() * B
// Here A is the local to global transform of the element and B is the local to global
// transform of the module.
- // If we are already in the module frame then there is nothing to do, we just return the
- // transform that is input. Otherwise we use the above formula.
Amg::Transform3D
SiDetectorElement::localToModuleFrame(const Amg::Transform3D& localTransform) const
{
- if (isModuleFrame()) {
- return localTransform;
- } else {
- return m_otherSide->transform().inverse() * transform() * localTransform * transform().inverse() * m_otherSide->transform();
- }
+ return moduleTransform().inverse() * transform() * localTransform * transform().inverse() * moduleTransform();
+ // if (isModuleFrame()) {
+ // return localTransform;
+ // } else {
+ // return m_otherSide->transform().inverse() * transform() * localTransform * transform().inverse() * m_otherSide->transform();
+ // }
}
Amg::Transform3D
SiDetectorElement::localToModuleTransform() const
{
- if (isModuleFrame()) {
- return Amg::Transform3D(); // Identity
- } else {
- return m_otherSide->transform().inverse() * transform();
- }
+ return transform().inverse() * moduleTransform();
+ // if (isModuleFrame()) {
+ // return Amg::Transform3D(); // Identity
+ // } else {
+ // return m_otherSide->transform().inverse() * transform();
+ // }
}
bool
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/python/FaserSCT_DigitizationConfigNew.py b/Tracker/TrackerDigitization/FaserSCT_Digitization/python/FaserSCT_DigitizationConfigNew.py
index a4bd9fb2c54f97e9983106eb6743c00bb772e585..c9f267da7d9d2ba0c6bc973c3fa949b65c217ecc 100644
--- a/Tracker/TrackerDigitization/FaserSCT_Digitization/python/FaserSCT_DigitizationConfigNew.py
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/python/FaserSCT_DigitizationConfigNew.py
@@ -240,9 +240,10 @@ def FaserSCT_OutputCfg(flags):
if flags.Digitization.TruthOutput:
ItemList += ["TrackerSimDataCollection#*"]
acc.merge(TruthDigitizationOutputCfg(flags))
- acc.merge(OutputStreamCfg(flags, "RDO", ItemList))
+ acc.merge(OutputStreamCfg(flags, "RDO"))
ostream = acc.getEventAlgo("OutputStreamRDO")
ostream.TakeItemsFromInput = True
+ ostream.ItemList += ItemList
return acc
diff --git a/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.cxx b/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.cxx
index 51514fb5b599630e9c32d04c406439c5ea891991..2149bdb33ed2ef55ea57a0c1bf55b025dd69dc16 100644
--- a/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.cxx
+++ b/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.cxx
@@ -3,7 +3,6 @@
*/
#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3.h"
-//#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHit_p2.h"
#include "FaserSiHitCollectionCnv.h"
@@ -18,14 +17,23 @@ FaserSiHitCollection_PERS* FaserSiHitCollectionCnv::createPersistent(FaserSiHitC
FaserSiHitCollection* FaserSiHitCollectionCnv::createTransient() {
MsgStream mlog(msgSvc(), "FaserSiHitCollectionConverter" );
FaserSiHitCollectionCnv_p3 converter_p3;
+ FaserSiHitCollectionCnv_p3a converter_p3a;
static const pool::Guid p3_guid("FF9508DE-3E25-425D-9556-16D319DCE0E1");
+ static const pool::Guid p3a_guid("72FD9F51-AB1B-4DF7-B430-6CCAE0A994DB");
FaserSiHitCollection *trans_cont(0);
- if( this->compareClassGuid(p3_guid)) {
+ if( this->compareClassGuid(p3a_guid))
+ {
+ std::unique_ptr< FaserSiHitCollection_p3a > col_vect( this->poolReadObject< FaserSiHitCollection_p3a >() );
+ trans_cont = converter_p3a.createTransient( col_vect.get(), mlog );
+ }
+ else if( this->compareClassGuid(p3_guid))
+ {
std::unique_ptr< FaserSiHitCollection_p3 > col_vect( this->poolReadObject< FaserSiHitCollection_p3 >() );
trans_cont = converter_p3.createTransient( col_vect.get(), mlog );
- } else {
+ }
+ else {
throw std::runtime_error("Unsupported persistent version of Data container");
}
return trans_cont;
diff --git a/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.h b/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.h
index 3a81e19fa0671a028948f7217aa256c06e94c9a1..47892fe8a5cd1738428aa90858c95e83d46fe232 100644
--- a/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.h
+++ b/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.h
@@ -8,12 +8,14 @@
#include "TrackerSimEvent/FaserSiHitCollection.h"
#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3.h"
#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3.h"
+#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h"
+#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h"
#include "AthenaPoolCnvSvc/T_AthenaPoolCustomCnv.h"
// Gaudi
#include "GaudiKernel/MsgStream.h"
// typedef to the latest persistent version
-typedef FaserSiHitCollection_p3 FaserSiHitCollection_PERS;
-typedef FaserSiHitCollectionCnv_p3 FaserSiHitCollectionCnv_PERS;
+typedef FaserSiHitCollection_p3a FaserSiHitCollection_PERS;
+typedef FaserSiHitCollectionCnv_p3a FaserSiHitCollectionCnv_PERS;
class FaserSiHitCollectionCnv : public T_AthenaPoolCustomCnv<FaserSiHitCollection, FaserSiHitCollection_PERS > {
friend class CnvFactory<FaserSiHitCollectionCnv>;
diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h
new file mode 100644
index 0000000000000000000000000000000000000000..8f453cf57ab52efb41fbe739d6543d75ba77c72a
--- /dev/null
+++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h
@@ -0,0 +1,41 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef FASERSIHITCOLLECTIONCNV_P3A_H
+#define FASERSIHITCOLLECTIONCNV_P3A_H
+
+// FaserSiHitCollectionCnv_p3, T/P separation of FaserSi Hits
+// author D.Costanzo <davide.costanzo@cern.ch>
+// author O.Arnaez <olivier.arnaez@cern.ch>
+
+#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h"
+#include "TrackerSimEvent/FaserSiHitCollection.h"
+#include "FaserSiHitCollection_p3a.h"
+
+
+class FaserSiHitCollectionCnv_p3a : public T_AthenaPoolTPCnvBase<FaserSiHitCollection, FaserSiHitCollection_p3a>
+{
+ public:
+
+ FaserSiHitCollectionCnv_p3a() {};
+
+ virtual FaserSiHitCollection* createTransient(const FaserSiHitCollection_p3a* persObj, MsgStream &log);
+
+ virtual void persToTrans(const FaserSiHitCollection_p3a* persCont,
+ FaserSiHitCollection* transCont,
+ MsgStream &log) ;
+ virtual void transToPers(const FaserSiHitCollection* transCont,
+ FaserSiHitCollection_p3a* persCont,
+ MsgStream &log) ;
+
+ private:
+
+ static const double m_persEneUnit;
+ static const double m_persLenUnit;
+ static const double m_persAngUnit;
+ static const double m_2bHalfMaximum;
+ static const int m_2bMaximum;
+};
+
+#endif
diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h
new file mode 100644
index 0000000000000000000000000000000000000000..137ebd68b84c0ea6e6643132023d90c49b23df0f
--- /dev/null
+++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h
@@ -0,0 +1,57 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef FASERSIHITCOLLECTION_P3A_H
+#define FASERSIHITCOLLECTION_P3A_H
+
+/*
+
+Authors: Davide Costanzo Rob Duxfield
+
+*/
+
+#include <vector>
+#include <string>
+
+class FaserSiHitCollection_p3a
+{
+ public:
+/// Default constructor
+ FaserSiHitCollection_p3a ();
+ //private:
+
+ std::vector<float> m_hit1_meanTime; // 1 element per string
+ std::vector<float> m_hit1_x0; //
+ std::vector<float> m_hit1_y0; //
+ std::vector<float> m_hit1_z0; //
+ std::vector<float> m_hit1_theta; //
+ std::vector<float> m_hit1_phi; //
+ std::vector<unsigned long> m_nHits; //
+
+ std::vector<unsigned short> m_hitEne_2b; // 1 element per hit
+ std::vector<unsigned short> m_hitLength_2b; //
+
+ std::vector<unsigned short> m_dTheta; // 1 element per hit except for first hit in string
+ std::vector<unsigned short> m_dPhi; //
+
+ std::vector<float> m_hitEne_4b; // 1 element per hit with m_hitEne_2b[i] == 2**16
+
+ std::vector<float> m_hitLength_4b; // 1 element per hit with m_hitLength_2b[i] == 2**16
+
+ std::vector<unsigned long> m_barcode;
+ std::vector<unsigned long> m_mcEvtIndex;
+ std::vector<char> m_evtColl;
+ std::vector<unsigned long> m_nBC;
+
+ std::vector<unsigned long> m_id;
+ std::vector<unsigned long> m_nId;
+};
+
+
+// inlines
+
+inline
+FaserSiHitCollection_p3a::FaserSiHitCollection_p3a () {}
+
+#endif
diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnvDict.h b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnvDict.h
index a730fca90074a9da71420e41651384173aa298cf..a898e3095bee7e4b25d2090b0e7365cbdcc4d00e 100644
--- a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnvDict.h
+++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnvDict.h
@@ -15,6 +15,8 @@
#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCnv_p2.h"
#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3.h"
#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3.h"
+#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h"
+#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h"
#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHit_p2.h"
#endif // TRACKERSIMEVENTTPCNV_TRACKERSIMEVENTTPCNVDICT_H
diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/selection.xml b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/selection.xml
index 7684404421c37267df94656b57ff9067dcce9503..6b34f88fcb38eb134efb2a283ac6635860a76928 100644
--- a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/selection.xml
+++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/selection.xml
@@ -3,4 +3,5 @@
<class name="FaserSiHit_p2" />
<class name="std::vector<FaserSiHit_p2>" />
<class name="FaserSiHitCollection_p3" id="FF9508DE-3E25-425D-9556-16D319DCE0E1" />
+ <class name="FaserSiHitCollection_p3a" id="72FD9F51-AB1B-4DF7-B430-6CCAE0A994DB" />
</lcgdict>
diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/src/TrackerHits/FaserSiHitCollectionCnv_p3a.cxx b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/src/TrackerHits/FaserSiHitCollectionCnv_p3a.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0aebca47dc6ae6eb0d6662d4a2b4abe1a52dc1ae
--- /dev/null
+++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/src/TrackerHits/FaserSiHitCollectionCnv_p3a.cxx
@@ -0,0 +1,351 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "TrackerSimEvent/FaserSiHit.h"
+#include "TrackerSimEvent/FaserSiHitCollection.h"
+#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h"
+#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h"
+#include "GeneratorObjects/HepMcParticleLink.h"
+
+#include <cmath>
+
+//CLHEP
+#include "CLHEP/Geometry/Point3D.h"
+// Gaudi
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/ThreadLocalContext.h"
+
+// Athena
+#include "StoreGate/StoreGateSvc.h"
+#include "AthenaKernel/ExtendedEventContext.h"
+
+
+// * * * stolen from eflowRec * * * //
+inline double phicorr(double a)
+{
+ if (a <= -M_PI)
+ {
+ return a+(2*M_PI*floor(-(a-M_PI)/(2*M_PI)));
+ }
+ else if (a > M_PI)
+ {
+ return a-(2*M_PI*floor((a+M_PI)/(2*M_PI)));
+ }
+ else
+ {
+ return a;
+ }
+}
+
+// * * * stolen from eflowRec * * * //
+inline double cycle(double a, double b)
+{
+ double del = b-a;
+ if (del > M_PI)
+ {
+ return a+2.0*M_PI;
+ }
+ else if (del < -M_PI)
+ {
+ return a-2.0*M_PI;
+ }
+ else
+ {
+ return a;
+ }
+}
+
+
+const double FaserSiHitCollectionCnv_p3a::m_persEneUnit = 1.0e-5;
+const double FaserSiHitCollectionCnv_p3a::m_persLenUnit = 1.0e-5;
+const double FaserSiHitCollectionCnv_p3a::m_persAngUnit = 1.0e-5;
+const double FaserSiHitCollectionCnv_p3a::m_2bHalfMaximum = pow(2.0, 15.0);
+const int FaserSiHitCollectionCnv_p3a::m_2bMaximum = (unsigned short)(-1);
+
+
+void FaserSiHitCollectionCnv_p3a::transToPers(const FaserSiHitCollection* transCont, FaserSiHitCollection_p3a* persCont, MsgStream &log)
+{
+ // Finds hits belonging to a "string" (in which the end point of one hit is the same as the start point of the next) and
+ // persistifies the end point of each hit plus the start point of the first hit in each string.
+ //
+ // Further compression is achieved by optimising the storage of the position vectors:- start (x,y,z) and (theta,phi) of
+ // first hit are stored as floats, (delta_theta,delta_phi) relative to the fisrst hit are stored as 2 byte numbers and
+ // used to specify the hit direction. All hit lengths are stored as 2 byte numbers.
+ //
+ // Additional savings are achieved by storing the energy loss for each hit as a 2 byte number and only storing the mean
+ // time of the first hit per string.
+ //
+ // See http://indico.cern.ch/getFile.py/access?contribId=11&resId=2&materialId=slides&confId=30893 for more info.
+
+ static const double dRcut = 1.0e-7;
+ static const double dTcut = 1.0;
+
+ const EventContext& ctx = Gaudi::Hive::currentContext();
+ const IProxyDict* proxy = Atlas::getExtendedEventContext(ctx).proxy();
+ const HepMcParticleLink * lastLink=nullptr;
+ int lastId = -1;
+ double stringFirstTheta = 0.0;
+ double stringFirstPhi = 0.0;
+ double lastT = 0.0;
+ double persSumE = 0.0;
+ double transSumE = 0.0;
+ unsigned int idx = 0;
+ unsigned int endBC = 0;
+ unsigned int endId = 0;
+ unsigned int endHit = 0;
+ HepGeom::Point3D<double> lastTransEnd(0.0, 0.0, 0.0);
+ HepGeom::Point3D<double> lastPersEnd(0.0, 0.0, 0.0);
+
+ for (FaserSiHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) {
+
+ FaserSiHitCollection::const_iterator siHit = it;
+
+
+ if ( !lastLink || (siHit->particleLink() != *lastLink) ) {
+
+ // store barcode, eventIndex and McEventCollection once for set of consecutive hits with same barcode
+
+ lastLink = &(siHit->particleLink());
+ persCont->m_barcode.push_back(lastLink->barcode());
+ // log << MSG::ALWAYS << "TP: Event collection from link: " << lastLink->getEventCollection() << " ; char value: " << lastLink->getEventCollectionAsChar() << " barcode: " << lastLink->barcode() << endmsg;
+ unsigned short index{0};
+ const HepMcParticleLink::index_type position =
+ HepMcParticleLink::getEventPositionInCollection(lastLink->eventIndex(),
+ lastLink->getEventCollection(),
+ proxy).at(0);
+ if (position!=0) {
+ index = lastLink->eventIndex();
+ if(lastLink->eventIndex()!=static_cast<HepMcParticleLink::index_type>(index)) {
+ log << MSG::WARNING << "Attempting to persistify an eventIndex larger than max unsigned short!" << endmsg;
+ }
+ }
+ persCont->m_mcEvtIndex.push_back(index);
+ persCont->m_evtColl.push_back(lastLink->getEventCollectionAsChar());
+
+ if (idx > 0) {
+ persCont->m_nBC.push_back(idx - endBC);
+ endBC = idx;
+ }
+ }
+
+ if ( (int)siHit->identify() != lastId ) {
+
+ // store id once for set of consecutive hits with same barcode
+
+ lastId = siHit->identify();
+ persCont->m_id.push_back(lastId);
+
+ if (idx > 0) {
+ persCont->m_nId.push_back(idx - endId);
+ endId = idx;
+ }
+ }
+
+ HepGeom::Point3D<double> st = siHit->localStartPosition();
+ HepGeom::Point3D<double> en = siHit->localEndPosition();
+
+ const double dx = st.x() - lastTransEnd.x();
+ const double dy = st.y() - lastTransEnd.y();
+ const double dz = st.z() - lastTransEnd.z();
+ const double t = siHit->meanTime();
+
+ const double dRLast = sqrt(dx * dx + dy * dy + dz * dz); // dR between end of previous hit and start of current one
+ const double dTLast = fabs(t - lastT);
+
+ CLHEP::Hep3Vector direction(0.0, 0.0, 0.0);
+ double theta = 0.0;
+ double phi = 0.0;
+ bool startNewString = false;
+
+ if (dRLast < dRcut && dTLast < dTcut) {
+
+ // hit is part of existing string
+
+ direction = CLHEP::Hep3Vector( en.x() - lastPersEnd.x(), en.y() - lastPersEnd.y(), en.z() - lastPersEnd.z() );
+
+ theta = direction.theta();
+ phi = phicorr( direction.phi() );
+
+ const int dTheta_2b = (int)( (theta - stringFirstTheta) / m_persAngUnit + m_2bHalfMaximum + 0.5 );
+ const int dPhi_2b = (int)( (cycle(phi, stringFirstPhi) - stringFirstPhi) / m_persAngUnit + m_2bHalfMaximum + 0.5 );
+
+ if ( dTheta_2b < m_2bMaximum && dTheta_2b >= 0 && dPhi_2b < m_2bMaximum && dPhi_2b >= 0) {
+ persCont->m_dTheta.push_back(dTheta_2b);
+ persCont->m_dPhi.push_back(dPhi_2b);
+ theta = stringFirstTheta + ( (double)dTheta_2b - m_2bHalfMaximum ) * m_persAngUnit;
+ phi = stringFirstPhi + ( (double)dPhi_2b - m_2bHalfMaximum ) * m_persAngUnit;
+ phi = phicorr(phi);
+ }
+ else {
+ startNewString = true;
+ }
+ }
+
+ if (startNewString || dRLast >= dRcut || dTLast >= dTcut) {
+
+ // begin new hit string
+
+ direction = CLHEP::Hep3Vector( en.x() - st.x(), en.y() - st.y(), en.z() - st.z() );
+
+ theta = direction.theta();
+ phi = phicorr( direction.phi() );
+
+ persCont->m_hit1_meanTime.push_back(t);
+ persCont->m_hit1_x0.push_back(st.x());
+ persCont->m_hit1_y0.push_back(st.y());
+ persCont->m_hit1_z0.push_back(st.z());
+ persCont->m_hit1_theta.push_back(theta);
+ persCont->m_hit1_phi.push_back(phi);
+
+ lastPersEnd = HepGeom::Point3D<double>(st.x(), st.y(), st.z());
+
+ stringFirstTheta = theta;
+ stringFirstPhi = phi;
+
+ if (idx > 0) {
+ persCont->m_nHits.push_back(idx - endHit);
+ endHit = idx;
+ }
+ }
+
+ lastTransEnd = HepGeom::Point3D<double>(en.x(), en.y(), en.z());
+ transSumE += siHit->energyLoss();
+
+ const int eneLoss_2b = (int)((transSumE - persSumE) / m_persEneUnit + 0.5); // calculated to allow recovery sum over
+ // whole hit string to chosen precision
+
+ const int hitLength_2b = (int)(direction.mag() / m_persLenUnit + 0.5); // calculated to give the correct position to
+ // the chosen precision, NOT the length of the
+ // hit (small difference in practice).
+ double eneLoss = 0.0;
+
+ if (eneLoss_2b >= m_2bMaximum) {
+ eneLoss = siHit->energyLoss();
+ persCont->m_hitEne_2b.push_back(m_2bMaximum);
+ persCont->m_hitEne_4b.push_back(eneLoss);
+ }
+ else {
+ eneLoss = eneLoss_2b * m_persEneUnit;
+ persCont->m_hitEne_2b.push_back(eneLoss_2b);
+ }
+
+ double length = 0.0;
+
+ if (hitLength_2b >= m_2bMaximum) {
+ length = direction.mag();
+ persCont->m_hitLength_2b.push_back(m_2bMaximum);
+ persCont->m_hitLength_4b.push_back(direction.mag());
+ }
+ else {
+ length = hitLength_2b * m_persLenUnit;
+ persCont->m_hitLength_2b.push_back(hitLength_2b);
+ }
+
+ CLHEP::Hep3Vector persDir(length, 0.0, 0.0);
+ persDir.setTheta(theta);
+ persDir.setPhi(phi);
+
+ lastPersEnd = (CLHEP::Hep3Vector)lastPersEnd + persDir;
+ persSumE += eneLoss;
+ lastT = t;
+
+ ++idx;
+ }
+
+ persCont->m_nBC.push_back(idx - endBC);
+ persCont->m_nId.push_back(idx - endId);
+ persCont->m_nHits.push_back(idx - endHit);
+
+ // log << MSG::ALWAYS << "nBC: " << persCont->m_nBC << endmsg;
+ // log << MSG::ALWAYS << "nId: " << persCont->m_nId << endmsg;
+ // log << MSG::ALWAYS << "nHits:" << persCont->m_nHits << endmsg;
+
+}
+
+
+FaserSiHitCollection* FaserSiHitCollectionCnv_p3a::createTransient(const FaserSiHitCollection_p3a* persObj, MsgStream &log) {
+ // for (size_t i = 0; i < persObj->m_evtColl.size(); i++)
+ // {
+ // if (persObj->m_evtColl[i] != 'a')
+ // log << MSG::ALWAYS << "Corrupted in createTransient: " << persObj->m_evtColl[i] << endmsg;
+ // }
+ std::unique_ptr<FaserSiHitCollection> trans(std::make_unique<FaserSiHitCollection>("DefaultCollectionName",persObj->m_nHits.size()));
+ persToTrans(persObj, trans.get(), log);
+ return(trans.release());
+}
+
+
+void FaserSiHitCollectionCnv_p3a::persToTrans(const FaserSiHitCollection_p3a* persCont, FaserSiHitCollection* transCont, MsgStream &/*log*/)
+{
+ const EventContext& ctx = Gaudi::Hive::currentContext();
+
+ unsigned int hitCount = 0;
+ unsigned int angleCount = 0;
+ unsigned int idxBC = 0;
+ unsigned int idxId = 0;
+ unsigned int idxEne4b = 0;
+ unsigned int idxLen4b = 0;
+ unsigned int endHit = 0;
+ unsigned int endBC = 0;
+ unsigned int endId = 0;
+
+ // log << MSG::ALWAYS << "nHits.size(): " << persCont->m_nHits.size() << " nBC.size(): " << persCont->m_nBC.size() << endmsg;
+ // log << MSG::ALWAYS << " nBC: " << persCont->m_nBC << endmsg;
+
+ for (unsigned int i = 0; i < persCont->m_nHits.size(); i++) {
+ if (persCont->m_nHits[i]) {
+ const unsigned int start = endHit;
+ endHit += persCont->m_nHits[i];
+ // log << MSG::ALWAYS << "i: " << i << " persCont->m_nHits[i]: " << persCont->m_nHits[i] << " start: " << start << " endHit: " << endHit << endmsg;
+
+ const double t0 = persCont->m_hit1_meanTime[i];
+ const double theta0 = persCont->m_hit1_theta[i];
+ const double phi0 = persCont->m_hit1_phi[i];
+ HepGeom::Point3D<double> endLast(persCont->m_hit1_x0[i], persCont->m_hit1_y0[i], persCont->m_hit1_z0[i]);
+
+ for (unsigned int j = start; j < endHit; j++) {
+
+ if (j >= endBC + persCont->m_nBC[idxBC])
+ endBC += persCont->m_nBC[idxBC++];
+
+ if (j >= endId + persCont->m_nId[idxId])
+ endId += persCont->m_nId[idxId++];
+
+ const double eneLoss_2b = persCont->m_hitEne_2b[hitCount];
+ const double hitLength_2b = persCont->m_hitLength_2b[hitCount];
+
+ const double eneLoss = (eneLoss_2b < m_2bMaximum) ? eneLoss_2b * m_persEneUnit : persCont->m_hitEne_4b[idxEne4b++];
+ const double length = (hitLength_2b < m_2bMaximum) ? hitLength_2b * m_persLenUnit : persCont->m_hitLength_4b[idxLen4b++];
+
+ const double dTheta = (j > start) ? ((double)persCont->m_dTheta[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0;
+ const double dPhi = (j > start) ? ((double)persCont->m_dPhi[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0;
+
+ const double meanTime = t0;
+ const double theta = theta0 + dTheta;
+ const double phi = phicorr(phi0 + dPhi);
+
+ CLHEP::Hep3Vector r(length, 0.0, 0.0);
+ r.setTheta(theta);
+ r.setPhi(phi);
+
+ HepGeom::Point3D<double> endThis( endLast + r );
+
+ HepMcParticleLink::PositionFlag flag = HepMcParticleLink::IS_INDEX;
+ if (persCont->m_mcEvtIndex[idxBC] == 0) {
+ flag = HepMcParticleLink::IS_POSITION;
+ }
+ // log << MSG::ALWAYS << "PT: Event collection from persCont: " << persCont->m_evtColl[idxBC] << " ; idxBC: " << idxBC << " evtColl.size(): " << persCont->m_evtColl.size() << " evtIndex.size(): " << persCont->m_mcEvtIndex.size() << endmsg;
+
+ HepMcParticleLink partLink( persCont->m_barcode[idxBC], persCont->m_mcEvtIndex[idxBC], HepMcParticleLink::ExtendedBarCode::eventCollectionFromChar(persCont->m_evtColl[idxBC]), flag, ctx );
+ transCont->Emplace( endLast, endThis, eneLoss, meanTime, partLink, persCont->m_id[idxId]);
+
+ endLast = endThis;
+
+ ++hitCount;
+ if (j > start) ++angleCount;
+ }
+ }
+ }
+ // log << MSG::ALWAYS << "hitCount: " << hitCount << endmsg;
+}
diff --git a/Waveform/WaveDigiTools/CMakeLists.txt b/Waveform/WaveDigiTools/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..692fdb69bc14451ba5a6a660d011bad5b14b66e5
--- /dev/null
+++ b/Waveform/WaveDigiTools/CMakeLists.txt
@@ -0,0 +1,25 @@
+################################################################################
+# Package: WaveDigiTools
+################################################################################
+
+# Declare the package name:
+atlas_subdir( WaveDigiTools )
+
+# External dependencies:
+find_package( ROOT )
+
+# Component(s) in the package:
+atlas_add_library( WaveDigiToolsLib
+ WaveDigiTools/*.h src/*.cxx src/*.h
+ PUBLIC_HEADERS WaveDigiTools
+ PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives WaveRawEvent
+ PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES}
+ )
+
+atlas_add_component( WaveDigiTools
+ src/components/*.cxx
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel WaveDigiToolsLib )
+
+
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..3e85ff839a04f6d02c794fa74567841450c9cfeb
--- /dev/null
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
@@ -0,0 +1,55 @@
+/*
+ Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file IWaveformDigitisationTool.h
+ * Header file for the IWaveformDigitisationTool class
+ * @author Carl Gwilliam, 2021
+ */
+
+
+#ifndef WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H
+#define WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H
+
+// Base class
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/IMessageSvc.h"
+#include "GaudiKernel/MsgStream.h"
+
+#include "WaveRawEvent/RawWaveformContainer.h"
+#include "WaveRawEvent/RawWaveform.h"
+
+#include "TF1.h"
+
+///Interface for waveform digitisation tools
+class IWaveformDigitisationTool : virtual public IAlgTool
+{
+public:
+
+ // InterfaceID
+ DeclareInterfaceID(IWaveformDigitisationTool, 1, 0);
+
+ IWaveformDigitisationTool():
+ m_msgSvc ( "MessageSvc", "ITrkEventCnvTool" )
+ {}
+
+ virtual ~IWaveformDigitisationTool() = default;
+
+ // Digitise HITS to Raw waveform
+ template<class CONT>
+ StatusCode digitise(const CONT* hitCollection,
+ RawWaveformContainer* waveContainer, TF1* kernel) const;
+
+private:
+ ServiceHandle<IMessageSvc> m_msgSvc;
+
+};
+
+#include "WaveDigiTools/IWaveformDigitisationTool.icc"
+
+
+#endif //WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
new file mode 100644
index 0000000000000000000000000000000000000000..57d4839bda5f2d30286ec412e3f01c92ce353b11
--- /dev/null
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
@@ -0,0 +1,51 @@
+#include <vector>
+#include <map>
+
+template<class CONT>
+StatusCode IWaveformDigitisationTool::digitise(const CONT* hitCollection,
+ RawWaveformContainer* container, TF1* kernel) const {
+
+
+ // Check the container
+ if (!container) {
+ MsgStream log(&(*m_msgSvc), name());
+ log << MSG::ERROR << "HitCollection passed to digitise() is null!" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ unsigned int size = 600; // TODO: how know the correct number of time samples?
+ std::vector<float> time(size);
+ 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
+
+ // Loop over time samples
+ for (const auto& t : time) {
+ std::map<unsigned int, float> counts;
+
+ // 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();
+ //std::cout << "HIT " << hit.identify() << " @ " << t << ": " << kernel->Eval(t) << " " << hit.energyLoss() << " -> " << counts[hit.identify()] << std::endl;
+ }
+
+ // Add count to correct waveform vec
+ for (const auto& c : counts) {
+ waveforms[c.first].push_back(baseline - c.second);
+ //std::cout << "ADC " << c.first << " @ " << t << ": " << baseline - c.second << std::endl;
+ }
+ }
+
+ // Loop over wavefrom vecs to make and store waveform
+ for (const auto& w : waveforms) {
+ RawWaveform* wfm = new RawWaveform();
+ wfm->setWaveform(0, w.second);
+ wfm->setIdentifier(Identifier(w.first));
+ wfm->setSamples(size);
+ container->push_back(wfm);
+ }
+
+
+ return StatusCode::SUCCESS;
+}
diff --git a/Waveform/WaveDigiTools/share/WaveformDigiAndRecoExample_jobOptions.py b/Waveform/WaveDigiTools/share/WaveformDigiAndRecoExample_jobOptions.py
new file mode 100644
index 0000000000000000000000000000000000000000..779031828c9c5736f9403e403514e19285029499
--- /dev/null
+++ b/Waveform/WaveDigiTools/share/WaveformDigiAndRecoExample_jobOptions.py
@@ -0,0 +1,87 @@
+#!/usr/bin/env python
+
+import sys
+
+if __name__ == "__main__":
+
+ fileroot = "output"
+ if len(sys.argv) > 1:
+ filename = sys.argv[1]
+
+ doRDO = False
+ if len(sys.argv) > 2:
+ filename = sys.argv[2]
+
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaCommon.Configurable import Configurable
+
+ Configurable.configurableRun3Behavior = True
+
+ log.setLevel(VERBOSE)
+
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion
+ ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
+ ConfigFlags.Input.ProjectName = "mc21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASER-TB00" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+
+ ConfigFlags.Input.Files = [
+ "/eos/project-f/faser-commissioning/Simulation/Test/TB.Elec.200GeV.SIM.root"
+ ]
+
+ if doRDO:
+ ConfigFlags.Output.RDOFileName = f"{fileroot}.RDO.root"
+ else:
+ ConfigFlags.addFlag("Output.xAODFileName", f"{fileroot}.xAOD.root")
+ ConfigFlags.Output.ESDFileName = f"{fileroot}.ESD.root"
+
+ ConfigFlags.lock()
+
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ acc = MainServicesCfg(ConfigFlags)
+
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+
+ acc.merge(PoolReadCfg(ConfigFlags))
+ acc.merge(PoolWriteCfg(ConfigFlags))
+
+ if doRDO:
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ itemList = [
+ "RawWaveformContainer#*"
+ ]
+ acc.merge(OutputStreamCfg(ConfigFlags, "RDO", itemList,disableEventTag=True))
+
+
+ else:
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ itemList = [
+ "xAOD::EventInfo#*",
+ "xAOD::WaveformHitContainer#*",
+ "xAOD::WaveformHitAuxContainer#*",
+ ]
+
+ acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList, disableEventTag=True))
+
+ from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg
+ acc.merge(ScintWaveformDigitizationCfg(ConfigFlags))
+
+ from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
+ acc.merge(CaloWaveformDigitizationCfg(ConfigFlags))
+
+ if not doRDO:
+ from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
+ acc.merge(WaveformReconstructionCfg(ConfigFlags))
+
+ #acc.foreach_component("*").OutputLevel = VERBOSE
+
+ # Execute and finish
+ sc = acc.run(maxEvents=100)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e4776da3bf1b140fa1c614642d939f46557b0ba6
--- /dev/null
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
@@ -0,0 +1,26 @@
+/*
+ Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file WaveformDigitisationTool.cxx
+ * Implementation file for the WaveformDigitisationTool class
+ * @ author C. Gwilliam, 2021
+ **/
+
+#include "WaveformDigitisationTool.h"
+
+// Constructor
+WaveformDigitisationTool::WaveformDigitisationTool(const std::string& type, const std::string& name, const IInterface* parent) :
+ base_class(type, name, parent)
+{
+}
+
+// Initialization
+StatusCode
+WaveformDigitisationTool::initialize() {
+ ATH_MSG_INFO( name() << "::initalize()" );
+ return StatusCode::SUCCESS;
+}
+
+
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a5ba71f3dd124fcdd2c6b4b8124ee96591512da
--- /dev/null
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
@@ -0,0 +1,36 @@
+/*
+ Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/** @file WaveformDigitisationTool.h
+ * Header file for WaveformDigitisationTool.h
+ *
+ */
+#ifndef WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
+#define WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
+
+//Athena
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "WaveDigiTools/IWaveformDigitisationTool.h"
+
+//Gaudi
+#include "GaudiKernel/ToolHandle.h"
+
+//STL
+
+class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisationTool> {
+ public:
+
+ /// Normal constructor for an AlgTool; 'properties' are also declared here
+ WaveformDigitisationTool(const std::string& type,
+ const std::string& name, const IInterface* parent);
+
+ /// Retrieve the necessary services in initialize
+ StatusCode initialize();
+
+ private:
+ // None
+
+};
+
+#endif // WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
diff --git a/Waveform/WaveDigiTools/src/components/WaveDigiTools_entries.cxx b/Waveform/WaveDigiTools/src/components/WaveDigiTools_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..169476d905bb3e3a56de0d15aac7b7d13c3bdeb1
--- /dev/null
+++ b/Waveform/WaveDigiTools/src/components/WaveDigiTools_entries.cxx
@@ -0,0 +1,3 @@
+#include "../WaveformDigitisationTool.h"
+
+DECLARE_COMPONENT( WaveformDigitisationTool )
diff --git a/Waveform/WaveEventCnv/WaveByteStream/python/WaveByteStreamConfig.py b/Waveform/WaveEventCnv/WaveByteStream/python/WaveByteStreamConfig.py
index cc13137d8e59eef7801d0ffa8fa60230721b372f..840a4d0f091a7e195f4e40f818c2b0b99c0b346b 100644
--- a/Waveform/WaveEventCnv/WaveByteStream/python/WaveByteStreamConfig.py
+++ b/Waveform/WaveEventCnv/WaveByteStream/python/WaveByteStreamConfig.py
@@ -37,6 +37,17 @@ def WaveByteStreamCfg(configFlags, **kwargs):
waveform_tool.PreshowerChannels = [12, 13]
acc.addPublicTool(waveform_tool)
+ elif configFlags.GeoModel.FaserVersion == "FASER-02":
+ print(" - setting up TI12 with IFT detector")
+
+ # Need to fix this!
+ waveform_tool = CompFactory.RawWaveformDecoderTool("RawWaveformDecoderTool")
+ waveform_tool.CaloChannels = [1, 0, 3, 2]
+ waveform_tool.VetoChannels = [4, 5, 6, 7]
+ waveform_tool.TriggerChannels = [9, 8, 11, 10]
+ waveform_tool.PreshowerChannels = [12, 13]
+ acc.addPublicTool(waveform_tool)
+
else:
print(" - unknown version: user must set up Waveform channel mapping by hand!")
diff --git a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py
index 39a99ed0c39d01e08c5f555e4e3efd98e0bc70a1..64764cd1c30ca7f039a28799dd16cc22bb50720c 100644
--- a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py
+++ b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py
@@ -11,24 +11,27 @@ WaveformReconstructionTool = CompFactory.WaveformReconstructionTool
ClockReconstructionTool = CompFactory.ClockReconstructionTool
# One stop shopping for normal FASER data
-def WaveformReconstructionCfg(flags, naive = True):
+def WaveformReconstructionCfg(flags, naive = False):
""" Return all algorithms and tools for Waveform reconstruction """
acc = ComponentAccumulator()
+ if not flags.Input.isMC:
+ acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg"))
+
if flags.Input.isMC and naive:
- if not "TB" in flags.GeoModel.FaserVersion:
+ if "TB" not in flags.GeoModel.FaserVersion:
acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoTimingHitWaveformRecAlg", "Trigger"))
acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoVetoHitToWaveformRecAlg", "Veto"))
acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoPresehowerHitWaveformRecAlg", "Preshower"))
acc.merge(PseudoCaloHitToWaveformRecCfg(flags, "PseudoCaloHitWaveformRecAlg"))
return acc
- acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg"))
-
+ if "TB" not in flags.GeoModel.FaserVersion:
+ acc.merge(WaveformHitRecCfg(flags, "TimingWaveformRecAlg", "Trigger"))
acc.merge(WaveformHitRecCfg(flags, "VetoWaveformRecAlg", "Veto"))
- acc.merge(WaveformHitRecCfg(flags, "TimingWaveformRecAlg", "Trigger"))
acc.merge(WaveformHitRecCfg(flags, "PreshowerWaveformRecAlg", "Preshower"))
acc.merge(WaveformHitRecCfg(flags, "CaloWaveformRecAlg", "Calo"))
+
return acc
# Return configured WaveformClock reconstruction algorithm
@@ -52,7 +55,13 @@ def WaveformHitRecCfg(flags, name="WaveformRecAlg", source="", **kwargs):
acc = ComponentAccumulator()
+ if flags.Input.isMC:
+ kwargs.setdefault("PeakThreshold", 5)
+
tool = WaveformReconstructionTool(name=source+"WaveformRecTool", **kwargs)
+
+ if "PeakThreshold" in kwargs: kwargs.pop("PeakThreshold")
+
kwargs.setdefault("WaveformContainerKey", source+"Waveforms")
kwargs.setdefault("WaveformHitContainerKey", source+"WaveformHits")
kwargs.setdefault("WaveformReconstructionTool", tool)
diff --git a/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.cxx b/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.cxx
index 5ead65452727c36307d8b4c038812ca9c3683f47..2074b566e044e8bc647237fa7ff5c127f9a94ffe 100644
--- a/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.cxx
+++ b/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.cxx
@@ -42,12 +42,6 @@ PseudoCaloHitToWaveformRecAlg::execute(const EventContext& ctx) const {
ATH_CHECK( caloHitHandle.isValid() );
ATH_MSG_DEBUG("Found ReadHandle for CaloHitCollection " << m_caloHitContainerKey);
- if (caloHitHandle->size() == 0) {
- ATH_MSG_DEBUG("CaloHitCollection found with zero length!");
- return StatusCode::SUCCESS;
- }
-
-
// Find the output waveform container
SG::WriteHandle<xAOD::WaveformHitContainer> waveformHitContainerHandle(m_waveformHitContainerKey, ctx);
ATH_CHECK( waveformHitContainerHandle.record( std::make_unique<xAOD::WaveformHitContainer>(),
@@ -55,6 +49,10 @@ PseudoCaloHitToWaveformRecAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("WaveformsHitContainer '" << waveformHitContainerHandle.name() << "' initialized");
+ if (caloHitHandle->size() == 0) {
+ ATH_MSG_DEBUG("CaloHitCollection found with zero length!");
+ return StatusCode::SUCCESS;
+ }
// "Reconstruct" the hits
CHECK( m_recoTool->reconstruct<CaloHitCollection>(caloHitHandle.ptr(),
diff --git a/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.cxx b/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.cxx
index 2b896e63c8a67b845072b74489ffe5a39b162751..4b8b719781f95a8881391738f076acf8e5c9204b 100644
--- a/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.cxx
+++ b/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.cxx
@@ -42,12 +42,6 @@ PseudoScintHitToWaveformRecAlg::execute(const EventContext& ctx) const {
ATH_CHECK( scintHitHandle.isValid() );
ATH_MSG_DEBUG("Found ReadHandle for ScintHitCollection " << m_scintHitContainerKey);
- if (scintHitHandle->size() == 0) {
- ATH_MSG_DEBUG("ScintHitCollection found with zero length!");
- return StatusCode::SUCCESS;
- }
-
-
// Find the output waveform container
SG::WriteHandle<xAOD::WaveformHitContainer> waveformHitContainerHandle(m_waveformHitContainerKey, ctx);
ATH_CHECK( waveformHitContainerHandle.record( std::make_unique<xAOD::WaveformHitContainer>(),
@@ -55,6 +49,10 @@ PseudoScintHitToWaveformRecAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("WaveformsHitContainer '" << waveformHitContainerHandle.name() << "' initialized");
+ if (scintHitHandle->size() == 0) {
+ ATH_MSG_DEBUG("ScintHitCollection found with zero length!");
+ return StatusCode::SUCCESS;
+ }
// "Reconstruct" the hits
CHECK( m_recoTool->reconstruct<ScintHitCollection>(scintHitHandle.ptr(),
diff --git a/Waveform/WaveRecTools/src/ClockReconstructionTool.cxx b/Waveform/WaveRecTools/src/ClockReconstructionTool.cxx
index aaae0861b6fab59b7144afedb72797ad9da91f86..c45ea51007362eecd0ab7827f50123c05ddebab7 100644
--- a/Waveform/WaveRecTools/src/ClockReconstructionTool.cxx
+++ b/Waveform/WaveRecTools/src/ClockReconstructionTool.cxx
@@ -79,8 +79,8 @@ ClockReconstructionTool::reconstruct(const RawWaveform& raw_wave,
// Get the coefficients
std::vector<double> re_full(N);
std::vector<double> im_full(N);
- std::vector<double> magnitude(N/2);
- fftr2c->GetPointsComplex(&re_full[0],&im_full[0]);
+ std::vector<double> magnitude(N/2+1); // From fftw manual, output array is N/2+1 long
+ fftr2c->GetPointsComplex(&re_full[0], &im_full[0]);
// Normalize the magnitude (just using the positive complex frequencies)
unsigned int i=0;
@@ -106,7 +106,7 @@ ClockReconstructionTool::reconstruct(const RawWaveform& raw_wave,
clockdata->set_dc_offset(magnitude[0]);
clockdata->set_frequency(imax * freqmult);
clockdata->set_amplitude(magnitude[imax]);
- clockdata->set_phase(atan2(im_full[imax], re_full[imax]));
+ clockdata->set_phase(atan2(im_full[imax], re_full[imax])); // Not a bug, atan2(y,x)!
ATH_MSG_DEBUG("Before correcting for finite resolution:");
ATH_MSG_DEBUG(*clockdata);
@@ -125,10 +125,11 @@ ClockReconstructionTool::reconstruct(const RawWaveform& raw_wave,
// Improved values
+ // Not a bug, atan2(y,x)!
double phase = atan2(im_full[imax], re_full[imax]) - dm * M_PI;
- // Fix any overflows
- if (phase < M_PI) phase += (2*M_PI);
- if (phase > M_PI) phase -= (2*M_PI);
+ // Fix any overflows caused by adding dm
+ if (phase < -M_PI) phase += (2*M_PI);
+ if (phase > +M_PI) phase -= (2*M_PI);
clockdata->set_frequency( (imax+dm) * freqmult );
clockdata->set_phase (phase);