diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_muon_fasernu_logE-101307.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_muon_fasernu_logE-101307.json
index db2954f994dd7c74c9fb30b8e9bda28ac3eea966..3e57b2d300eac4ae2caa74cdccd746b2f24c997a 100644
--- a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_muon_fasernu_logE-101307.json
+++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_muon_fasernu_logE-101307.json
@@ -4,7 +4,7 @@
"maxE": 5000.0,
"minE": 10.0,
"pid": [-13, 13],
- "radius": -25.,
+ "radius": -25.0,
"angle": 0.0006,
"run": 101307,
"sampler": "log",
diff --git a/Derviation/DerivationAlgs/CMakeLists.txt b/Derviation/DerivationAlgs/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d5ef7b4b3e64410b5db7c39a08bbd16acec22562
--- /dev/null
+++ b/Derviation/DerivationAlgs/CMakeLists.txt
@@ -0,0 +1,15 @@
+################################################################################
+# Package: DerivationAlgs
+################################################################################
+
+# Declare the package name:
+atlas_subdir( DerivationAlgs )
+
+atlas_add_component( DerivationAlgs
+ src/*.cxx src/*.h
+ src/components/*.cxx
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel StoreGateLib DerivationToolsLib)
+
+atlas_install_python_modules( python/*.py )
+
+
diff --git a/Derviation/DerivationAlgs/python/DerivationAlgsConfig.py b/Derviation/DerivationAlgs/python/DerivationAlgsConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..a38fb966c9bcfb155b7a7951386988a4fff3235b
--- /dev/null
+++ b/Derviation/DerivationAlgs/python/DerivationAlgsConfig.py
@@ -0,0 +1,79 @@
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+def DerivationAlgCfg(flags, stream, tools, **kwargs):
+
+ acc = ComponentAccumulator()
+
+ kwargs.setdefault("Tools", tools)
+ ##print ("SETTING", kwargs["Tools"], type(kwargs["Tools"][0]))
+ acc.addEventAlgo(CompFactory.Derive(stream + "_DerivationAlg", **kwargs))
+
+ return acc
+
+def DerivationOutputCfg(flags, stream, accept, items = [], extra_items = [], exclude_items = [], **kwargs):
+
+ acc = ComponentAccumulator()
+
+ if not items:
+ items = [ "xAOD::EventInfo#*"
+ , "xAOD::EventAuxInfo#*"
+ , "xAOD::FaserTriggerData#*"
+ , "xAOD::FaserTriggerDataAux#*"
+ , "FaserSiHitCollection#*" # Strip hits, do we want this?
+ , "FaserSCT_RDO_Container#*"
+ , "xAOD::WaveformHitContainer#*"
+ , "xAOD::WaveformHitAuxContainer#*"
+ , "xAOD::WaveformClock#*"
+ , "xAOD::WaveformClockAuxInfo#*"
+ ]
+
+ if not items and extra_items:
+ items.append(extra_items)
+
+ # from PrimaryDPDMaker/python/Primary_DPD_OutputDefinitions.py in athena
+ # Once can use TakeItemsFromInput = True probably need to use
+ # acc.getEventAlgo(f"OutputStream{stream}").RemoveItem(exclude_list)
+ exclude_list = []
+ if exclude_items:
+ for ex in exclude_items:
+ if ex in items:
+ exclude_list.append(ex)
+
+ if ex.endswith("*"):
+ for it in items:
+ if it.startswith(ex.rstrip("*")):
+ exclude_list.append(it)
+
+ # items = list(set(items) - set(exclude_list))
+
+
+ #flags.unlock()
+ #flags.addFlag(f"Output.AOD{stream}FileName", f"my.{stream}.xAOD.root")
+ #flags.lock()
+
+ acc.merge(OutputStreamCfg(flags, stream, items))
+ acc.getEventAlgo(f"OutputStream{stream}").AcceptAlgs = [accept]
+ ## if not items:
+ ##acc.getEventAlgo(f"OutputStream{stream}").TakeItemsFromInput = True # crashes
+ ## if extra_items:
+ ##cc.getEventAlgo(f"OutputStream{stream}").RemoveItem(extra_items)
+ ## if exclude_list:
+ ##acc.getEventAlgo(f"OutputStream{stream}").AddItem(exclude_list)
+
+ return acc
+
+def FullyConfiguredStream(flags, stream, tools, items = [], extra_items = [], **kwargs):
+ # TODO:
+ # - get items from input + why crash
+
+ acc = ComponentAccumulator()
+
+ acc.merge(DerivationAlgCfg(flags, stream, tools, **kwargs))
+ acc.merge(DerivationOutputCfg(flags, stream, stream + "_DerivationAlg", items, extra_items))
+
+ return acc
+
+
+
diff --git a/Derviation/DerivationAlgs/share/runDerive.py b/Derviation/DerivationAlgs/share/runDerive.py
new file mode 100644
index 0000000000000000000000000000000000000000..b8b22c187c15a39993a68aa722e5d6bae3362f7b
--- /dev/null
+++ b/Derviation/DerivationAlgs/share/runDerive.py
@@ -0,0 +1,147 @@
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+
+
+def DerivationAlgCfg(flags, name, frac, **kwargs):
+
+ acc = ComponentAccumulator()
+
+ tool = CompFactory.ExampleDerivationTool(name + "_TestTool", SaveFraction = frac)
+ print ("ZEBRA", tool.SaveFraction)
+
+ kwargs.setdefault("Tools", [tool])
+ acc.addEventAlgo(CompFactory.Derive(name, **kwargs))
+
+ return acc
+
+def DerivationAlgCfg2(flags, name, **kwargs):
+
+ acc = ComponentAccumulator()
+
+ tool = CompFactory.TriggerStreamTool(name + "_TriggerSteamTool")
+ kwargs.setdefault("Tools", [tool])
+ acc.addEventAlgo(CompFactory.Derive(name, **kwargs))
+
+ return acc
+
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ ConfigFlags.Input.Files = [
+ "/eos/experiment/faser/rec/2022/p0008/007984/Faser-Physics-007984-00000-p0008-xAOD.root"
+ #"/bundle/data/FASER/Ti12data/filter/r0008/007983/Faser-Physics-007983-TrigMask08-r0008-xAOD.root"
+ ]
+
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.addFlag("Output.AODSTREAM1FileName", "my.STREAM1.xAOD.root")
+ ConfigFlags.addFlag("Output.AODSTREAM2FileName", "my.STREAM2.xAOD.root")
+ ConfigFlags.addFlag("Output.AODSTREAM3FileName", "my.STREAM3.xAOD.root")
+ #ConfigFlags.Output.STREAM1FileName = fileName
+
+ ConfigFlags.lock()
+
+ # Core components
+ cfg = MainServicesCfg(ConfigFlags)
+ cfg.merge(PoolReadCfg(ConfigFlags))
+ cfg.merge(PoolWriteCfg(ConfigFlags))
+
+ # Derivation alg
+ cfg.merge(DerivationAlgCfg(ConfigFlags, "DerivationAlg1", 10))
+ cfg.merge(DerivationAlgCfg(ConfigFlags, "DerivationAlg2", 90))
+ cfg.merge(DerivationAlgCfg2(ConfigFlags, "DerivationAlg3"))
+
+ # Writing
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ streamName1 = "AODSTREAM1" # Needs to have AOD in name
+ itemList1 = [ "xAOD::EventInfo#*"
+ , "xAOD::EventAuxInfo#*"
+ , "xAOD::FaserTriggerData#*"
+ , "xAOD::FaserTriggerDataAux#*"
+ , "FaserSiHitCollection#*" # Strip hits, do we want this?
+ , "FaserSCT_RDO_Container#*"
+ , "xAOD::WaveformHitContainer#*"
+ , "xAOD::WaveformHitAuxContainer#*"
+ , "xAOD::WaveformClock#*"
+ , "xAOD::WaveformClockAuxInfo#*"
+# , "FaserSCT_SpacePointContainer#*" # Crashes
+# , "Tracker::FaserSCT_ClusterContainer#*"
+# , "TrackCollection#*"
+ ]
+
+ cfg.merge(OutputStreamCfg(ConfigFlags, streamName1, itemList1)) #, disableEventTag = True))
+ cfg.getEventAlgo("OutputStreamAODSTREAM1").AcceptAlgs = ["DerivationAlg1"]
+ #cfg.getEventAlgo("OutputStreamAODSTREAM1").TakeItemsFromInput = True
+
+
+ # Writing
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ streamName2 = "AODSTREAM2" # Needs to have AOD in name
+ itemList2 = [ "xAOD::EventInfo#*"
+ , "xAOD::EventAuxInfo#*"
+ , "xAOD::FaserTriggerData#*"
+ , "xAOD::FaserTriggerDataAux#*"
+ ]
+ cfg.merge(OutputStreamCfg(ConfigFlags, streamName2, itemList2)) #, disableEventTag = True))
+ cfg.getEventAlgo("OutputStreamAODSTREAM2").AcceptAlgs = ["DerivationAlg2"]
+
+ # Writing
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ streamName3 = "AODSTREAM3" # Needs to have AOD in name
+ itemList3 = [ "xAOD::EventInfo#*"
+ , "xAOD::EventAuxInfo#*"
+ , "xAOD::FaserTriggerData#*"
+ , "xAOD::FaserTriggerDataAux#*"
+ ]
+ cfg.merge(OutputStreamCfg(ConfigFlags, streamName3, itemList3)) #, disableEventTag = True))
+ cfg.getEventAlgo("OutputStreamAODSTREAM3").AcceptAlgs = ["DerivationAlg3"]
+
+
+
+# from OutputStreamAthenaPool.MultipleStreamManager import MSMgr
+# streamName = "STREAM1"
+# fileName = "streaming.STREAM1.root"
+# testStream = MSMgr.NewPoolRootStream(streamName, fileName)
+# testStream.AcceptAlgs(["DerivationAlg1"])
+# cfg.addEventAlgo(testStream)
+
+# # 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
+
+ # Execute and finish
+ cfg.printConfig()
+
+
+ sc = cfg.run(maxEvents=1000)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Derviation/DerivationAlgs/share/run_streaming.py b/Derviation/DerivationAlgs/share/run_streaming.py
new file mode 100644
index 0000000000000000000000000000000000000000..7ec4d11aa0f18643dfc8ee12e8ebf11e12b29166
--- /dev/null
+++ b/Derviation/DerivationAlgs/share/run_streaming.py
@@ -0,0 +1,92 @@
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ ConfigFlags.Input.Files = [
+ "/eos/experiment/faser/rec/2022/p0008/007984/Faser-Physics-007984-00000-p0008-xAOD.root"
+ ]
+
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ for stream in ["STREAM1", "STREAM2", "STREAM3"]:
+ ConfigFlags.addFlag(f"Output.AOD{stream}FileName", f"my.{stream}.xAOD.root")
+
+ ConfigFlags.lock()
+
+ # Core components
+ cfg = MainServicesCfg(ConfigFlags)
+ cfg.merge(PoolReadCfg(ConfigFlags))
+ cfg.merge(PoolWriteCfg(ConfigFlags))
+
+ # Derivations
+
+ from DerivationAlgs.DerivationAlgsConfig import FullyConfiguredStream
+
+ name = "STREAM1"
+
+
+ cfg.merge(FullyConfiguredStream(ConfigFlags, stream = name,
+ tools = [CompFactory.ExampleDerivationTool(name + "_TestTool", SaveFraction = 10.)])
+ )
+
+ name = "STREAM2"
+ cfg.merge(FullyConfiguredStream(ConfigFlags, stream = name,
+ tools = [ CompFactory.ExampleDerivationTool(name + "_TestTool", SaveFraction = 90.)],
+ items = [ "xAOD::EventInfo#*"
+ , "xAOD::EventAuxInfo#*"
+ , "xAOD::FaserTriggerData#*"
+ , "xAOD::FaserTriggerDataAux#*"
+ ])
+ )
+
+ name = "STREAM3"
+ cfg.merge(FullyConfiguredStream(ConfigFlags, stream = name,
+ tools = [CompFactory.TriggerStreamTool(name + "_TriggerTool")],
+ items = [ "xAOD::EventInfo#*"
+ , "xAOD::EventAuxInfo#*"
+ , "xAOD::FaserTriggerData#*"
+ , "xAOD::FaserTriggerDataAux#*"
+ ])
+ )
+
+
+# # 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
+
+ # Execute and finish
+ cfg.printConfig(withDetails = True, summariseProps = True, printDefaults = True)
+
+
+ sc = cfg.run(maxEvents=1000)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Derviation/DerivationAlgs/src/Derive.cxx b/Derviation/DerivationAlgs/src/Derive.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b34492787d313c955bff97b12c82fb402cdbc98f
--- /dev/null
+++ b/Derviation/DerivationAlgs/src/Derive.cxx
@@ -0,0 +1,53 @@
+#include "Derive.h"
+
+
+Derive::Derive(const std::string& name,
+ ISvcLocator* pSvcLocator)
+ : AthFilterAlgorithm(name, pSvcLocator) {
+
+ //declareProperty("Tools", m_tools);
+
+}
+
+StatusCode
+Derive::initialize() {
+ ATH_MSG_INFO(name() << "::initalize()" );
+
+ ATH_CHECK( m_tools.retrieve() );
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+Derive::finalize() {
+ ATH_MSG_INFO(name() << "::finalize()");
+ ATH_MSG_INFO("Derivation" << name() << " accepted " << m_passed << " out of " << m_events << " events");
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+Derive::execute() {
+ ATH_MSG_DEBUG("Executing ... ");
+
+ m_events++;
+
+ bool acceptEvent(true);
+
+ for (auto& tool : m_tools) {
+
+ // Skimming - remove events
+ if (!tool->passed()) acceptEvent = false;
+
+ // Thinning - remove info from event
+ ATH_CHECK(tool->removeBranch());
+
+ // Augmenting - add info to an event
+ ATH_CHECK(tool->addBranch());
+
+ }
+
+ setFilterPassed(acceptEvent);
+ if (acceptEvent) m_passed++;
+
+ return StatusCode::SUCCESS;
+}
diff --git a/Derviation/DerivationAlgs/src/Derive.h b/Derviation/DerivationAlgs/src/Derive.h
new file mode 100644
index 0000000000000000000000000000000000000000..03b971dd0852e54ea7a90e5ffb1135aedcba07a7
--- /dev/null
+++ b/Derviation/DerivationAlgs/src/Derive.h
@@ -0,0 +1,64 @@
+#ifndef DERIVATIONALGS_DERIVE_H
+#define DERIVATIONALGS_DERIVE_H
+
+// Base class
+#include "AthenaBaseComps/AthFilterAlgorithm.h"
+
+// FASER
+#include "DerivationTools/IDerivationTool.h"
+
+// Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/IAlgTool.h"
+
+// Athena
+#include "xAODEventInfo/EventInfo.h"
+#include "StoreGate/ReadHandleKey.h"
+
+// ROOT
+#include <TRandom3.h>
+
+// STL
+
+class Derive : public AthFilterAlgorithm {
+
+public:
+ // Constructor
+ Derive(const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~Derive() = default;
+
+ /** @name Usual algorithm methods */
+ //@{
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute() override;
+ virtual StatusCode finalize() override;
+ //@}
+
+
+ private:
+
+ /** @name Disallow default instantiation, copy, assignment */
+ //@{
+ Derive() = delete;
+ Derive(const Derive&) = delete;
+ Derive &operator=(const Derive&) = delete;
+ //@}
+
+ ///
+
+ /** @name Steerable tools */
+ //@
+ ToolHandleArray<IDerivationTool> m_tools {this, "Tools", {}, "List of tools"};
+ //@}
+
+ /** Number of events processed */
+ int m_events {0};
+
+ /** Number of events selected */
+ int m_passed {0};
+
+};
+
+
+#endif // DERIVATIONALGS_DERIVE_H
diff --git a/Derviation/DerivationAlgs/src/components/DerivationAlgs_entries.cxx b/Derviation/DerivationAlgs/src/components/DerivationAlgs_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c46a518fefa207cc65c90781eb1660a9438c8e32
--- /dev/null
+++ b/Derviation/DerivationAlgs/src/components/DerivationAlgs_entries.cxx
@@ -0,0 +1,3 @@
+#include "../Derive.h"
+
+DECLARE_COMPONENT( Derive )
diff --git a/Derviation/DerivationTools/CMakeLists.txt b/Derviation/DerivationTools/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c167cb1598233f25fb51a5bcbea28c6f17a22c40
--- /dev/null
+++ b/Derviation/DerivationTools/CMakeLists.txt
@@ -0,0 +1,24 @@
+################################################################################
+# Package: DerivationTools
+################################################################################
+
+# Declare the package name:
+atlas_subdir( DerivationTools )
+
+# External dependencies:
+find_package( ROOT )
+
+# Component(s) in the package:
+atlas_add_library( DerivationToolsLib
+ DerivationTools/*.h src/*.cxx src/*.h
+ PUBLIC_HEADERS DerivationTools
+ PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES AthenaBaseComps AthenaKernel xAODFaserTrigger
+ PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES}
+ )
+
+atlas_add_component( DerivationTools
+ src/components/*.cxx
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel xAODEventInfo xAODFaserTrigger DerivationToolsLib)
+
diff --git a/Derviation/DerivationTools/DerivationTools/IDerivationTool.h b/Derviation/DerivationTools/DerivationTools/IDerivationTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..29970fe09230f02c11aea8b5c8663540bc8312fc
--- /dev/null
+++ b/Derviation/DerivationTools/DerivationTools/IDerivationTool.h
@@ -0,0 +1,48 @@
+/*
+ Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file IDerivationTool.h
+ * Header file for the IDerivationTool class
+ * @author Carl Gwilliam, 2021
+ */
+
+
+#ifndef DERIVATIONTOOLS_IDERIVATIONTOOL_H
+#define DERIVATIONTOOLS_IDERIVATIONTOOL_H
+
+// Gaudi
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+
+///Interface for derivation tools
+class IDerivationTool : virtual public IAlgTool
+{
+public:
+
+ // InterfaceID
+ DeclareInterfaceID(IDerivationTool, 1, 0);
+
+ virtual ~IDerivationTool() = default;
+
+ // Apply skimming
+ virtual bool passed() = 0;
+
+ /// Apply thinning
+ virtual StatusCode removeBranch() = 0;
+
+ /// Apply augmentation
+ virtual StatusCode addBranch() = 0;
+
+
+private:
+ // None
+
+};
+
+
+
+
+#endif //DERIVATIONTOOLS_IDERIVATIONTOOL_H
diff --git a/Derviation/DerivationTools/src/ExampleDerivationTool.cxx b/Derviation/DerivationTools/src/ExampleDerivationTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6f27c1861eba2d5ac40c617dc8ed17be2085b80b
--- /dev/null
+++ b/Derviation/DerivationTools/src/ExampleDerivationTool.cxx
@@ -0,0 +1,48 @@
+/*
+ Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file ExamplederivationTool.cxx
+ * Implementation file for the ExamplederivationTool class
+ * @ author C. Gwilliam, 2021
+ **/
+
+#include "ExampleDerivationTool.h"
+
+// Constructor
+ExampleDerivationTool::ExampleDerivationTool(const std::string& type, const std::string& name, const IInterface* parent) :
+ base_class(type, name, parent)
+{
+ //std::cout << "CTOR " << name << " with SaveFraction = " << m_fraction << std::endl;
+}
+
+// Initialization
+StatusCode
+ExampleDerivationTool::initialize() {
+ ATH_MSG_INFO( name() << "::initalize()" );
+ //std::cout << "INIT " << name() << " with SaveFraction = " << m_fraction << std::endl;
+
+ return StatusCode::SUCCESS;
+}
+
+bool
+ExampleDerivationTool::passed(){
+
+ bool accept(false);
+
+ m_events++;
+
+ float frac = ((float)(m_passed+1))/(float)m_events * 100.0;
+
+ if (frac > m_fraction) {
+ ATH_MSG_DEBUG("Filter failed " << m_passed << " " << m_events << " " << frac << " " << m_fraction);
+ accept = false;
+ } else {
+ ATH_MSG_DEBUG("Filter passed " << m_passed << " " << m_events << " " << frac << " " << m_fraction);
+ accept = true;
+ m_passed++;
+ }
+
+ return accept;
+}
diff --git a/Derviation/DerivationTools/src/ExampleDerivationTool.h b/Derviation/DerivationTools/src/ExampleDerivationTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..d89657fff0d590ec21320d318d3123a5553fd307
--- /dev/null
+++ b/Derviation/DerivationTools/src/ExampleDerivationTool.h
@@ -0,0 +1,53 @@
+/*
+ Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/** @file TriggerStreamTool.h
+ * Header file for TriggerStreamTool.h
+ *
+ */
+#ifndef DERIVATIONTOOLS_EXAMPLEDERIVATIONTOOL_H
+#define DERIVATIONTOOLS_EXAMPLEDERIVATIONTOOL_H
+
+//Athena
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "DerivationTools/IDerivationTool.h"
+
+//Gaudi
+#include "GaudiKernel/ToolHandle.h"
+
+//STL
+
+class ExampleDerivationTool: public extends<AthAlgTool, IDerivationTool> {
+ public:
+
+ /// Normal constructor for an AlgTool; 'properties' are also declared here
+ ExampleDerivationTool(const std::string& type,
+ const std::string& name, const IInterface* parent);
+
+ /// Retrieve the necessary services in initialize
+ StatusCode initialize();
+
+ // Apply skimming
+ bool passed();
+
+ /// Apply thinning
+ StatusCode removeBranch() {return StatusCode::SUCCESS;}
+
+ /// Apply augmentation
+ StatusCode addBranch() {return StatusCode::SUCCESS;}
+
+ private:
+
+ /** Fraction of events to save */
+ Gaudi::Property<float> m_fraction {this, "SaveFraction", 100, "Fraction of events to save"};
+
+ /** Number of events processed */
+ int m_events {0};
+
+ /** Number of events selected */
+ int m_passed {0};
+
+};
+
+#endif // WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
diff --git a/Derviation/DerivationTools/src/TriggerStreamTool.cxx b/Derviation/DerivationTools/src/TriggerStreamTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8345766320effb1ea0cac4fa1e70d155e407c069
--- /dev/null
+++ b/Derviation/DerivationTools/src/TriggerStreamTool.cxx
@@ -0,0 +1,35 @@
+/*
+ Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file TriggerStreamTool.cxx
+ * Implementation file for the TriggerStreamTool class
+ * @ author C. Gwilliam, 2021
+ **/
+
+#include "TriggerStreamTool.h"
+
+// Constructor
+TriggerStreamTool::TriggerStreamTool(const std::string& type, const std::string& name, const IInterface* parent) :
+ base_class(type, name, parent)
+{
+}
+
+// Initialization
+StatusCode
+TriggerStreamTool::initialize() {
+ ATH_MSG_INFO( name() << "::initalize() with mask = " << m_mask );
+
+ ATH_CHECK( m_triggerDataKey.initialize() );
+
+ return StatusCode::SUCCESS;
+}
+
+bool
+TriggerStreamTool::passed() {
+
+ SG::ReadHandle<xAOD::FaserTriggerData> triggerData(m_triggerDataKey);
+
+ return triggerData->tap() & m_mask;
+}
diff --git a/Derviation/DerivationTools/src/TriggerStreamTool.h b/Derviation/DerivationTools/src/TriggerStreamTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..bb229f60418707c2ad56eae717b7e0daca6dc26e
--- /dev/null
+++ b/Derviation/DerivationTools/src/TriggerStreamTool.h
@@ -0,0 +1,56 @@
+/*
+ Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/** @file TriggerStreamTool.h
+ * Header file for TriggerStreamTool.h
+ *
+ */
+#ifndef DERIVATIONTOOLS_TRIGGERSTREAMTOOL_H
+#define DERIVATIONTOOLS_TRIGGERSTREAMTOOL_H
+
+
+// FASER
+#include "StoreGate/ReadHandleKey.h"
+#include "xAODFaserTrigger/FaserTriggerData.h"
+
+//Athena
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "DerivationTools/IDerivationTool.h"
+
+//Gaudi
+#include "GaudiKernel/ToolHandle.h"
+
+//STL
+
+class TriggerStreamTool: public extends<AthAlgTool, IDerivationTool> {
+ public:
+
+ /// Normal constructor for an AlgTool; 'properties' are also declared here
+ TriggerStreamTool(const std::string& type,
+ const std::string& name, const IInterface* parent);
+
+ /// Retrieve the necessary services in initialize
+ StatusCode initialize();
+
+ // Apply skimming
+ bool passed();
+
+ /// Apply thinning
+ StatusCode removeBranch() {return StatusCode::SUCCESS;}
+
+ /// Apply augmentation
+ StatusCode addBranch() {return StatusCode::SUCCESS;}
+
+ private:
+
+ /** Trigger condition to apply */
+ Gaudi::Property<int> m_mask {this, "TriggerMask", 0x8, "Trigger mask to apply"};
+
+ /// StoreGate key
+ SG::ReadHandleKey<xAOD::FaserTriggerData> m_triggerDataKey
+ { this, "FaserTriggerDataKey", "FaserTriggerData", "ReadHandleKey for xAOD::FaserTriggerData"};
+
+};
+
+#endif // WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
diff --git a/Derviation/DerivationTools/src/components/DerivationTools_entries.cxx b/Derviation/DerivationTools/src/components/DerivationTools_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..91e58e594e7864bd61ddf7d196f1625d361a4c8b
--- /dev/null
+++ b/Derviation/DerivationTools/src/components/DerivationTools_entries.cxx
@@ -0,0 +1,5 @@
+#include "../TriggerStreamTool.h"
+#include "../ExampleDerivationTool.h"
+
+DECLARE_COMPONENT( TriggerStreamTool )
+DECLARE_COMPONENT( ExampleDerivationTool )
diff --git a/DetectorDescription/GeoModel/FaserGeoModel/python/GeoModelInit.py b/DetectorDescription/GeoModel/FaserGeoModel/python/GeoModelInit.py
index d17c8e39d932df7be5a7265d7b4c2403a0546430..fbee2144d985412d2624d31e745eec881f3fb54d 100644
--- a/DetectorDescription/GeoModel/FaserGeoModel/python/GeoModelInit.py
+++ b/DetectorDescription/GeoModel/FaserGeoModel/python/GeoModelInit.py
@@ -97,7 +97,7 @@ def _setupGeoModel():
#conddb.addFolderSplitOnline("SCT","/Tracker/Onl/Align","/Tracker/Align",className="AlignableTransformContainer")
print("Override Alignment dbname to OFLP200, fix this when alignment available in CONDBR3")
- conddb.addFolder("/Tracker/Align", "SCT_OFL",className="AlignableTransformContainer",db="OFLP200")
+ conddb.addFolder("/Tracker/Align", "SCT_OFL",className="AlignableTransformContainer")
from AthenaCommon.AlgSequence import AthSequencer
condSeq = AthSequencer("AthCondSeq")
if not hasattr(condSeq, "FaserSCT_AlignCondAlg"):
diff --git a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..944ed9689053915dd526251fabf1c5b1af5dfa4e
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
@@ -0,0 +1,12 @@
+atlas_subdir(NtupleDumper)
+
+atlas_add_component(
+ NtupleDumper
+ src/NtupleDumperAlg.h
+ src/NtupleDumperAlg.cxx
+ src/component/NtupleDumper_entries.cxx
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserTrigger ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint
+)
+
+atlas_install_python_modules(python/*.py)
+atlas_install_scripts(scripts/*.py)
diff --git a/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py b/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..5bc52c41df03c934465f97658a39aeb342350631
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py
@@ -0,0 +1,104 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+
+def NtupleDumperAlgCfg(flags, **kwargs):
+ # Initialize GeoModel
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+
+ acc.merge(MagneticFieldSvcCfg(flags))
+ # acc.merge(FaserActsTrackingGeometrySvcCfg(flags))
+ # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+
+ actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
+ actsExtrapolationTool.MaxSteps = 1000
+ actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+
+ NtupleDumperAlg = CompFactory.NtupleDumperAlg("NtupleDumperAlg",**kwargs)
+ NtupleDumperAlg.ExtrapolationTool = actsExtrapolationTool
+ acc.addEventAlgo(NtupleDumperAlg)
+
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output += ["HIST2 DATAFILE='Data-tuple.root' OPT='RECREATE'"]
+ acc.addService(thistSvc)
+
+ return acc
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+ # Set up logging and new style config
+ log.setLevel(INFO)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ ConfigFlags.Input.Files = [
+ '/eos/experiment/faser/rec/2022/p0008//008119/Faser-Physics-008119-00168-p0008-xAOD.root',
+
+
+ ]
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, UseFlukaWeights=True))
+
+ # silencio
+ AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
+ AthenaEventLoopMgr.EventPrintoutInterval=500
+ acc.addService(AthenaEventLoopMgr)
+
+ # # 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
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/analyzeNtuple.py b/PhysicsAnalysis/NtupleDumper/scripts/analyzeNtuple.py
new file mode 100755
index 0000000000000000000000000000000000000000..7d735f1599e455001162fff0e5a33acad16d2e35
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/analyzeNtuple.py
@@ -0,0 +1,235 @@
+#!/usr/bin/env python
+
+# Set up (Py)ROOT.
+import ROOT
+import glob
+import sys
+import pandas as pd
+
+
+# Define a Landau convoluted with a gaussian for MIP fitting
+landguas_conv = ROOT.TF1Convolution("landau","gaus",-10,100,True) # the functions have to be close to zero at min and max bin of convolution or else circular Fourier transform will move convolve values at max and min
+landguas_conv.SetNofPointsFFT(10000)
+landgaus = ROOT.TF1("landgaus",landguas_conv, -10, 100, landguas_conv.GetNpar())
+landgaus.SetParNames("Landau constant","Landau MPV","Landau width","Gaussian mean","Gaussian width")
+
+user_input = str(sys.argv[1]) # set to either 'all_high', 'all_low', or a run number
+
+t = ROOT.TChain("nt")
+nfiles = 0
+all_run_paths = glob.glob("/eos/project/f/faser-commissioning/DeionsNtuples/*")
+
+if user_input=="all_high":
+ runconfig = "High_gain"
+ print("processing high-gain runs")
+ gain = 30.0
+ for run_path in all_run_paths:
+ nfiles += t.Add(run_path+"/Data-tuple-High_gain*.root") # chain all ntuples from all runs that are high gain
+ rootFile = ROOT.TFile("/eos/project/f/faser-commissioning/DeionsNtuples/7930/Data-tuple-High_gain-007930-00000-100.root"); # load file from largest high gain run to get noise histograms
+elif user_input=="all_low":
+ runconfig = "Low_gain"
+ print("processing low-gain runs")
+ gain = 1.0
+ for run_path in all_run_paths:
+ nfiles += t.Add(run_path+"/Data-tuple-Low_gain*.root") # chain all ntuples from all runs that are high gain
+ rootFile = ROOT.TFile("/eos/project/f/faser-commissioning/DeionsNtuples/8137/Data-tuple-Low_gain-008137-00000-100.root"); # load file from largest low gain run to get noise histograms
+else: # assume user_input is a run number
+ # get run configuration from table oon Brian's website
+ table_runs = pd.read_html('http://aagaard.web.cern.ch/aagaard/FASERruns.html') # load in run tables from website
+ df = table_runs[0] # table_runs is a list of all tables on the website, we only want the first one
+ runconfig=str(df.at[df.loc[df['Run'] == int(user_input)].index[0],'Configuration'].replace(' ','_')) # get config from website run log telling if run is High_gain or Low_gain calo
+ print("processing run "+runconfig+" ("+runconfig+")")
+ if runconfig=="High_gain":
+ gain = 30.0
+ elif runconfig=="Low_gain":
+ gain = 1.0
+ else:
+ print("run config is neither 'High_gain' nor 'Low_gain', calo histogram ranges may be messed up")
+ gain = 1.0 # assume low gain
+
+ nfiles += t.Add("/eos/project/f/faser-commissioning/DeionsNtuples/"+user_input+"/*.root") # chain all ntuples from all runs that are high gain
+ rootFile = ROOT.TFile("/eos/project/f/faser-commissioning/DeionsNtuples/"+user_input+"/Data-tuple-"+runconfig+"-00"+user_input+"-00000-100.root"); # load file from largest low gain run to get noise histograms
+
+
+
+
+print("number of files chained together = ",nfiles)
+
+#ROOT.gROOT.SetStyle("ATLAS")
+#ROOT.gStyle.SetOptStat(111110) #take away option box in histograms
+#ROOT.gStyle.SetOptTitle(1)
+#ROOT.gStyle.SetOptFit(1)
+
+# Define histograms here
+hCaloCharge = []
+hCaloPeak = []
+hXYvsEcalo = []
+for chan in range(4):
+ hCaloCharge.append(ROOT.TH1F("hCalo"+str(chan)+"charge", "Charge in calo ch"+str(chan)+";Q (pC);# of events",100,0.2*gain,2.0*gain))
+ hCaloPeak.append(ROOT.TH1F("hCalo"+str(chan)+"peak", "Peak in calo ch"+str(chan)+";peak (mV);# of events",100,1.0*gain,5.0*gain))
+ hXYvsEcalo.append(ROOT.TProfile2D("hXYvsEcalo"+str(chan)+"" , "Calo ch"+str(chan)+" Charge vs Pos;X pos (mm);Y pos (mm)",26, -130.0, 130.0, 26, -130.0, 130.0))
+
+hCaloChargeTotal = ROOT.TH1F("hCaloChargeTotal", "Charge in Calo;Charge (pC);# of events",100,0.2*gain,2.0*gain)
+hCaloEdep = ROOT.TH1F("hCaloEdep", "Edep in Calo;Edep (GeV);# of events",100,0.0,1.8)
+
+hCaloThetaX = ROOT.TH1F("hCaloThetaX", "Track #theta_{x} at Calo face;#theta_{x} (radians);# of tracks",100,-0.1,0.1)
+hCaloThetaY = ROOT.TH1F("hCaloThetaY", "Track #theta_{y} at Calo face;#theta_{y} (radians);# of tracks",100,-0.1,0.1)
+
+hTrackPvsPYdiff = ROOT.TProfile("hTrackPvsPYdiff" , "Track #Deltap_{Y}/p vs p;Track p (MeV);(pY_{upstream} - pY_{downstream}) / p_{total}",100, 1000.0, 200000.0)
+hTrackPvsPXdiff = ROOT.TProfile("hTrackPvsPXdiff" , "Track #Deltap_{X}/p vs p;Track p (MeV);(pX_{upstream} - pX_{downstream}) / p_{total}",100, 1000.0, 200000.0)
+
+#t.Print() # will show you all variables in ttree
+
+i = 0
+for event in t:
+ i += 1
+
+ if i%1000 == 0:
+ print( "Processing event #%i of %i" % (i, t.GetEntries() ) )
+
+ if event.longTracks > 0: # only process events with at least one track that has hits in last 3 tracking stations
+ for j in range(event.longTracks): # loop over all long tracks in the event (long = has hits in last 3 tracking stations)
+ if event.Track_p0[j] != 0.0:
+ hTrackPvsPYdiff.Fill(event.Track_p0[j],(event.Track_py0[j] - event.Track_py1[j])/event.Track_p0[j])
+ hTrackPvsPXdiff.Fill(event.Track_p0[j],(event.Track_px0[j] - event.Track_px1[j])/event.Track_p0[j])
+
+ #print("track charge = %i and nLayers = %i" % (event.Track_charge[j],event.Track_nLayers[j]))
+ #print("track upstream (x,y,z) (px,py,pz) = (%f,%f,%f) (%f,%f,%f)" % (event.Track_x0[j],event.Track_y0[j],event.Track_z0[j],event.Track_px0[j],event.Track_py0[j],event.Track_pz0[j]))
+ #print("track downstream (x,y,z) (px,py,pz) = (%f,%f,%f) (%f,%f,%f)" % (event.Track_x1[j],event.Track_y1[j],event.Track_z1[j],event.Track_px1[j],event.Track_py1[j],event.Track_pz1[j]))
+
+ #print("track at vetoNu (x,y) (thetaX,thetaY) = (%f,%f) (%f,%f)" % (event.Track_X_atVetoNu[j],event.Track_Y_atVetoNu[j],event.Track_ThetaX_atVetoNu[j],event.Track_ThetaY_atVetoNu[j]))
+ #print("track at Calo (x,y) (thetaX,thetaY) = (%f,%f) (%f,%f)" % (event.Track_X_atCalo[j],event.Track_Y_atCalo[j],event.Track_ThetaX_atCalo[j],event.Track_ThetaY_atCalo[j]))
+
+ #print("number of track segments = ",event.TrackSegments)
+ #for j in range(event.TrackSegments):
+ #print("trackseg (x,y,z) (px,py,pz) = (%f,%f,%f) (%f,%f,%f)" % (event.TrackSegment_x[j],event.TrackSegment_y[j],event.TrackSegment_z[j],event.TrackSegment_px[j],event.TrackSegment_py[j],event.TrackSegment_pz[j]))
+ #print("trackseg chi2 = %i and ndof = %i" % (event.TrackSegment_Chi2[j],event.TrackSegment_nDoF[j]))
+
+ #print("number of SpacePoints = ",event.SpacePoints)
+ #for j in range(event.SpacePoints):
+ #print("Spacepoint #",j)
+ #print("SpacePoint (x,y,z) = (%f,%f,%f)" % (event.SpacePoint_x[j],event.SpacePoint_y[j],event.SpacePoint_z[j]))
+
+ hCaloEdep.Fill(event.Calo_total_Edep)
+ hCaloChargeTotal.Fill(event.Calo_total_charge)
+
+ x_calo = event.Track_X_atCalo[0]
+ y_calo = event.Track_Y_atCalo[0]
+
+ hCaloThetaX.Fill(event.Track_ThetaX_atCalo[0])
+ hCaloThetaY.Fill(event.Track_ThetaY_atCalo[0])
+
+ if abs(event.Track_ThetaX_atCalo[0]) > 0.1 or abs(event.Track_ThetaX_atCalo[0]) > 0.1: continue
+
+ for chan,charge in enumerate([event.Calo0_raw_charge,event.Calo1_raw_charge,event.Calo2_raw_charge,event.Calo3_raw_charge]):
+ if charge > 0.2*gain and charge < 2.0*gain:
+ hXYvsEcalo[chan].Fill(x_calo,y_calo,charge)
+
+ if x_calo > -60.0 and x_calo < -20.0 and y_calo > -80.0 and y_calo < -10.0:
+ hCaloCharge[0].Fill(event.Calo0_raw_charge)
+ hCaloPeak[0].Fill(event.Calo0_raw_peak)
+ elif x_calo > 70.0 and x_calo < 100.0 and y_calo > -90.0 and y_calo < -10.0:
+ hCaloCharge[1].Fill(event.Calo1_raw_charge)
+ hCaloPeak[1].Fill(event.Calo1_raw_peak)
+ elif x_calo > -60.0 and x_calo < -20.0 and y_calo > 20.0 and y_calo < 110.0:
+ hCaloCharge[2].Fill(event.Calo2_raw_charge)
+ hCaloPeak[2].Fill(event.Calo2_raw_peak)
+ elif x_calo > 70.0 and x_calo < 100.0 and y_calo > 20.0 and y_calo < 110.0:
+ hCaloCharge[3].Fill(event.Calo3_raw_charge)
+ hCaloPeak[3].Fill(event.Calo3_raw_peak)
+
+# if i > 10000:
+# break
+
+# create a list of histograms of random event integrals
+hRandomCharge = []
+for chan in range(15):
+ hRandomCharge.append(rootFile.Get("hRandomCharge"+str(chan)))
+
+# Now make some plots
+filename = "analyze-"+runconfig+"-Ntuples.pdf"
+
+c = ROOT.TCanvas()
+c.Print(filename+'[')
+hCaloEdep.Draw()
+ROOT.gPad.SetLogy()
+c.Print(filename)
+
+c = ROOT.TCanvas()
+hCaloChargeTotal.Draw()
+c.Print(filename)
+
+c = ROOT.TCanvas()
+c.Divide(2,2)
+for chan in range(4):
+ c.cd(1+chan)
+ hXYvsEcalo[chan].GetZaxis().SetRangeUser(hCaloCharge[chan].GetMean() - 0.3*hCaloCharge[chan].GetStdDev(),hCaloCharge[chan].GetMean() + 0.4*hCaloCharge[chan].GetStdDev())
+ hXYvsEcalo[chan].Draw('COLZ')
+c.Print(filename)
+
+leg = []
+c = ROOT.TCanvas()
+c.Divide(2,2)
+for chan in range(4):
+ c.cd(1+chan)
+ hCaloCharge[chan].Fit("landau")
+ landgaus.SetParameters(hCaloCharge[chan].GetFunction("landau").GetParameter(0),hCaloCharge[chan].GetFunction("landau").GetParameter(1),hCaloCharge[chan].GetFunction("landau").GetParameter(2),0.0,hRandomCharge[chan].GetStdDev())
+ landgaus.SetParLimits(0,0.1*hCaloCharge[chan].GetFunction("landau").GetParameter(0),20.0*hCaloCharge[chan].GetFunction("landau").GetParameter(0))
+ landgaus.SetParLimits(1,0.5*hCaloCharge[chan].GetFunction("landau").GetParameter(1),1.2*hCaloCharge[chan].GetFunction("landau").GetParameter(1))
+ landgaus.SetParLimits(2,0.1*hCaloCharge[chan].GetFunction("landau").GetParameter(2),1.2*hCaloCharge[chan].GetFunction("landau").GetParameter(2))
+ landgaus.FixParameter(3,0.0)
+ landgaus.FixParameter(4,hRandomCharge[chan].GetStdDev()) # fix gaussian smearing to the noise seen in randomly triggered events
+ hCaloCharge[chan].Fit("landgaus","+")
+ hCaloCharge[chan].GetFunction("landgaus").SetLineColor(4)
+ hCaloCharge[chan].Draw()
+
+ leg.append( ROOT.TLegend(0.55,0.55,0.89,0.75) )
+ leg[chan].AddEntry(hCaloCharge[chan].GetFunction("landau"),"Landau MPV = "+str(hCaloCharge[chan].GetFunction("landau").GetParameter(1))[:6]+" #pm "+str(hCaloCharge[chan].GetFunction("landau").GetParError(1))[:6],"L")
+ leg[chan].AddEntry(hCaloCharge[chan].GetFunction("landgaus"),"Landguas MPV = "+str(hCaloCharge[chan].GetFunction("landgaus").GetParameter(1))[:6]+" #pm "+str(hCaloCharge[chan].GetFunction("landgaus").GetParError(1))[:6],"L")
+ leg[chan].AddEntry(hCaloCharge[chan].GetFunction("landgaus"),"Landguas gaussian width = "+str(hCaloCharge[chan].GetFunction("landgaus").GetParameter(4))[:6],"")
+ leg[chan].SetBorderSize(0)
+ leg[chan].Draw()
+c.Print(filename)
+
+leg = []
+c = ROOT.TCanvas()
+c.Divide(2,2)
+for chan in range(4):
+ c.cd(1+chan)
+ hCaloPeak[chan].Fit("landau")
+ hCaloPeak[chan].Draw()
+
+ leg.append( ROOT.TLegend(0.55,0.55,0.89,0.75) )
+ leg[chan].AddEntry(hCaloPeak[chan].GetFunction("landau"),"Landau MPV = "+str(hCaloPeak[chan].GetFunction("landau").GetParameter(1))[:6]+" #pm "+str(hCaloPeak[chan].GetFunction("landau").GetParError(1))[:6],"L")
+ leg[chan].SetBorderSize(0)
+ leg[chan].Draw()
+c.Print(filename)
+
+c = ROOT.TCanvas()
+c.Divide(1,2)
+c.cd(1)
+hCaloThetaX.Draw()
+c.cd(2)
+hCaloThetaY.Draw()
+c.Print(filename)
+
+c = ROOT.TCanvas()
+c.Divide(1,2)
+c.cd(1)
+hTrackPvsPYdiff.GetYaxis().SetRangeUser(hTrackPvsPYdiff.GetMean(2) - hTrackPvsPYdiff.GetStdDev(2), hTrackPvsPYdiff.GetMean(2) + hTrackPvsPYdiff.GetStdDev(2))
+hTrackPvsPYdiff.Draw()
+c.cd(2)
+hTrackPvsPXdiff.GetYaxis().SetRangeUser(hTrackPvsPXdiff.GetMean(2) - hTrackPvsPXdiff.GetStdDev(2), hTrackPvsPXdiff.GetMean(2) + hTrackPvsPXdiff.GetStdDev(2))
+hTrackPvsPXdiff.Draw()
+c.Print(filename)
+
+c = ROOT.TCanvas()
+c.Divide(4,4)
+for chan in range(15):
+ c.cd(1+chan)
+ hRandomCharge[chan].Draw()
+c.Print(filename)
+
+# Must close file at the end
+c.Print(filename+']')
+
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py
new file mode 100755
index 0000000000000000000000000000000000000000..c8ab9bb76af0be2c537fd201c6ffb2387bbbd319
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py
@@ -0,0 +1,127 @@
+#!/usr/bin/env python3
+
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the FASER collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+
+
+def NtupleDumperAlgCfg(flags, OutName, **kwargs):
+ # Initialize GeoModel
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+
+ acc.merge(MagneticFieldSvcCfg(flags))
+ # acc.merge(FaserActsTrackingGeometrySvcCfg(flags))
+ # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+
+ actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
+ actsExtrapolationTool.MaxSteps = 10000
+ actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+
+ NtupleDumperAlg = CompFactory.NtupleDumperAlg("NtupleDumperAlg",**kwargs)
+ NtupleDumperAlg.ExtrapolationTool = actsExtrapolationTool
+ acc.addEventAlgo(NtupleDumperAlg)
+
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output += [f"HIST2 DATAFILE='{OutName}' OPT='RECREATE'"]
+ acc.addService(thistSvc)
+
+ return acc
+
+if __name__ == "__main__":
+
+ import glob
+ import sys
+ import ROOT
+
+ runno=int(sys.argv[1])
+ num=int(sys.argv[2])
+ filesPerJob=int(sys.argv[3])
+ run_config=str(sys.argv[4])
+
+ ptag="p0008"
+
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ # Set up logging and new style config
+ log.setLevel(INFO)
+ Configurable.configurableRun3Behavior = True
+
+ dataDir=f"/eos/experiment/faser/rec/2022/{ptag}/{runno:06d}"
+ files=sorted(glob.glob(f"{dataDir}/Faser-Physics*"))
+ fileListInitial=files[num*filesPerJob:(num+1)*filesPerJob]
+ fileList=[]
+ for fName in fileListInitial:
+ try:
+ fh=ROOT.TFile(fName)
+ fileList.append(fName)
+ except OSError:
+ print("Warning bad file: ",fName)
+
+ log.info(f"Analyzing Run {runno} files {num*filesPerJob} to {(num+1)*filesPerJob} (num={num})")
+ log.info(f"Got {len(fileList)} files out of {len(fileListInitial)}")
+
+ outName=f"Data-tuple-{run_config}-{runno:06d}-{num:05d}-{filesPerJob}.root"
+
+ # Configure
+ ConfigFlags.Input.Files = fileList
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outName, UseFlukaWeights=True, CaloConfig=run_config))
+
+ AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
+ AthenaEventLoopMgr.EventPrintoutInterval=1000
+ acc.addService(AthenaEventLoopMgr)
+
+ # # 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
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.sh b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.sh
new file mode 100755
index 0000000000000000000000000000000000000000..dc42bc521c678a3bed8536e93b6b68ad25549ca7
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+runno=$1
+num=$2
+filesPerJob=$3
+runconfig=$4
+
+WD=$PWD
+
+cd /afs/cern.ch/user/d/dfellers/faser
+source setup.sh
+cd build
+source x86_64-centos7-gcc11-opt/setup.sh
+cd $WD
+echo "Starting analysis"
+analyzeRun.py $runno $num $filesPerJob $runconfig
+cp Data-tuple*.root /eos/project/f/faser-commissioning/DeionsNtuples/$runno/
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/submitAllJobsThatHaveErrorLogs.py b/PhysicsAnalysis/NtupleDumper/scripts/submitAllJobsThatHaveErrorLogs.py
new file mode 100755
index 0000000000000000000000000000000000000000..5f0805b7d876e4573c5cbab096deba0364314726
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/submitAllJobsThatHaveErrorLogs.py
@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+
+import glob
+import os
+import sys
+import pandas as pd
+
+print("NtupleDumping all condor jobs that produced non-empty error logs")
+
+table_runs = pd.read_html('http://aagaard.web.cern.ch/aagaard/FASERruns.html') # load in run tables from website
+df = table_runs[0] # table_runs is a list of all tables on the website, we only want the first table
+
+# make a liist of all runs and batch job numbers that failed and thus have error logs that are not empty
+run_list = []
+allErrorLogs_list = glob.glob('/afs/cern.ch/user/d/dfellers/faser/ntuple-dumper/run/logs/*/*err')
+for error_log in allErrorLogs_list:
+ if os.path.getsize(error_log) != 0:
+ print('Error Log is not empty: ', error_log)
+ run_num = int(error_log.split('/')[-2].split('-')[-1])
+ bath_num = int(error_log.split('.')[-2])
+ run_list.append([run_num,bath_num])
+
+print("list to be re-submitted:", run_list)
+
+ptag="p0008"
+filesPerJob=100
+
+for i,run_info in enumerate(run_list):
+ runno = run_info[0]
+ batch_number = run_info[1]
+
+ runconfig=str(df.at[df.loc[df['Run'] == runno].index[0],'Configuration'].replace(' ','_')) # get config from website run log telling if run is High_gain or Low_gain calo
+
+ print("%i of %i runs processed. Currently processing run %i-%i (%s)"%(i,len(run_list),runno,batch_number,runconfig))
+
+ batchFile=f"batch/Run-{runno:06d}-{batch_number}.sub"
+ fh=open(batchFile,"w")
+ pwd=os.getenv("PWD")
+ fh.write(f"""
+ executable = {pwd}/analyzeRun.sh
+ arguments = {runno} {batch_number} {filesPerJob} {runconfig}
+ output = {pwd}/logs/Run-{runno:06d}/batch.{batch_number}.out
+ error = {pwd}/logs/Run-{runno:06d}/batch.{batch_number}.err
+ log = {pwd}/logs/Run-{runno:06d}/batch.log
+ requirements = (Arch == "X86_64" && OpSysAndVer =?= "CentOS7")
+ getenv = False
+ transfer_output_files = ""
+ +JobFlavour = "workday"
+ queue 1
+ """)
+ fh.close()
+ os.system(f"echo condor_submit {batchFile}")
+ os.system(f"condor_submit {batchFile}")
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/submitAllStableRuns.py b/PhysicsAnalysis/NtupleDumper/scripts/submitAllStableRuns.py
new file mode 100755
index 0000000000000000000000000000000000000000..92f49d3392476a245892a5a322c8cfd5c352fde3
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/submitAllStableRuns.py
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+
+import glob
+import os
+import sys
+import pandas as pd
+
+print("NtupleDumping all stable-beam runs found at http://aagaard.web.cern.ch/aagaard/FASERruns.html")
+
+table_runs = pd.read_html('http://aagaard.web.cern.ch/aagaard/FASERruns.html') # load in run tables from website
+df = table_runs[0] # table_runs is a list of all tables on the website, we only want the first table
+df.columns = [c.replace(' ', '_') for c in df.columns] # rename the columns such that names with spaces are replaced with '_' (needed to access 'Stable_Beam' column)
+df.drop(df[df.Stable_Beam != 'Yes'].index, inplace=True) # drop all runs that are not stable beam runs
+df.drop(df[(df.Configuration != 'Low gain') & (df.Configuration != 'High gain')].index, inplace=True) # drop all runs that are not 'Low gain' or 'High gain'
+
+run_list = df['Run'].tolist()
+
+ptag="p0008"
+filesPerJob=100
+
+for i,runno in enumerate(run_list):
+ runconfig=str(df.at[df.loc[df['Run'] == runno].index[0],'Configuration'].replace(' ','_')) # get config from website run log telling if run is High_gain or Low_gain calo
+
+ print("%i of %i runs processed. Currently processing run %i (%s)"%(i,len(run_list),runno,runconfig))
+
+ os.system(f"mkdir -p logs/Run-{runno:06d}")
+ os.system(f"mkdir -p batch")
+ os.system(f"mkdir -p /eos/project/f/faser-commissioning/DeionsNtuples/{runno}")
+
+ dataDir=f"/eos/experiment/faser/rec/2022/{ptag}/{runno:06d}"
+ files=glob.glob(f"{dataDir}/Faser-Physics*")
+ numFiles=len(files)
+ numJobs=numFiles//filesPerJob+(numFiles%filesPerJob!=0)
+ batchFile=f"batch/Run-{runno:06d}.sub"
+ fh=open(batchFile,"w")
+ pwd=os.getenv("PWD")
+ fh.write(f"""
+ executable = {pwd}/analyzeRun.sh
+ arguments = {runno} $(ProcId) {filesPerJob} {runconfig}
+ output = {pwd}/logs/Run-{runno:06d}/batch.$(ProcId).out
+ error = {pwd}/logs/Run-{runno:06d}/batch.$(ProcId).err
+ log = {pwd}/logs/Run-{runno:06d}/batch.log
+ requirements = (Arch == "X86_64" && OpSysAndVer =?= "CentOS7")
+ getenv = False
+ transfer_output_files = ""
+ +JobFlavour = "workday"
+ queue {numJobs}
+ """)
+ fh.close()
+ os.system(f"echo condor_submit {batchFile}")
+ os.system(f"condor_submit {batchFile}")
+
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/submitRun.py b/PhysicsAnalysis/NtupleDumper/scripts/submitRun.py
new file mode 100755
index 0000000000000000000000000000000000000000..b408759504ab5a1afc1df332154250cfcd7b4325
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/submitRun.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+
+import glob
+import os
+import sys
+import pandas as pd
+
+table_runs = pd.read_html('http://aagaard.web.cern.ch/aagaard/FASERruns.html') # load in run tables from website
+df = table_runs[0] # table_runs is a list of all tables on the website, we only want the first one
+
+ptag="p0008"
+filesPerJob=100
+
+runno=int(sys.argv[1])
+runconfig=str(df.at[df.loc[df['Run'] == runno].index[0],'Configuration'].replace(' ','_')) # get config from website run log telling if run is High_gain or Low_gain calo
+
+os.system(f"mkdir -p logs/Run-{runno:06d}")
+os.system(f"mkdir -p batch")
+os.system(f"mkdir -p /eos/project/f/faser-commissioning/DeionsNtuples/{runno}")
+
+dataDir=f"/eos/experiment/faser/rec/2022/{ptag}/{runno:06d}"
+files=glob.glob(f"{dataDir}/Faser-Physics*")
+numFiles=len(files)
+numJobs=numFiles//filesPerJob+(numFiles%filesPerJob!=0)
+batchFile=f"batch/Run-{runno:06d}.sub"
+fh=open(batchFile,"w")
+pwd=os.getenv("PWD")
+fh.write(f"""
+executable = {pwd}/analyzeRun.sh
+arguments = {runno} $(ProcId) {filesPerJob} {runconfig}
+output = {pwd}/logs/Run-{runno:06d}/batch.$(ProcId).out
+error = {pwd}/logs/Run-{runno:06d}/batch.$(ProcId).err
+log = {pwd}/logs/Run-{runno:06d}/batch.log
+requirements = (Arch == "X86_64" && OpSysAndVer =?= "CentOS7")
+getenv = False
+transfer_output_files = ""
++JobFlavour = "workday"
+queue {numJobs}
+""")
+fh.close()
+os.system(f"echo condor_submit {batchFile}")
+
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d1a07f428c73dfb7d6c7f1bf69ab20b73f12a196
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -0,0 +1,890 @@
+#include "NtupleDumperAlg.h"
+#include "TrkTrack/Track.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "ScintIdentifier/VetoNuID.h"
+#include "ScintIdentifier/VetoID.h"
+#include "ScintIdentifier/TriggerID.h"
+#include "ScintIdentifier/PreshowerID.h"
+#include "FaserCaloIdentifier/EcalID.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "TrackerSpacePoint/FaserSCT_SpacePoint.h"
+#include "Identifier/Identifier.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "xAODTruth/TruthParticle.h"
+#include <cmath>
+#include <TH1F.h>
+
+
+NtupleDumperAlg::NtupleDumperAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator),
+ AthHistogramming(name),
+ m_histSvc("THistSvc/THistSvc", name) {}
+
+
+void NtupleDumperAlg::addBranch(const std::string &name,
+ float* var) {
+ m_tree->Branch(name.c_str(),var,(name+"/F").c_str());
+}
+void NtupleDumperAlg::addBranch(const std::string &name,
+ unsigned int* var) {
+ m_tree->Branch(name.c_str(),var,(name+"/I").c_str());
+}
+
+void NtupleDumperAlg::addWaveBranches(const std::string &name,
+ int nchannels,
+ int first) {
+ for(int ch=0;ch<nchannels;ch++) {
+ std::string base=name+std::to_string(ch)+"_";
+ addBranch(base+"time",&m_wave_localtime[first]);
+ addBranch(base+"peak",&m_wave_peak[first]);
+ addBranch(base+"width",&m_wave_width[first]);
+ addBranch(base+"charge",&m_wave_charge[first]);
+ addBranch(base+"raw_peak",&m_wave_raw_peak[first]);
+ addBranch(base+"raw_charge",&m_wave_raw_charge[first]);
+ addBranch(base+"baseline",&m_wave_baseline_mean[first]);
+ addBranch(base+"baseline_rms",&m_wave_baseline_rms[first]);
+ addBranch(base+"status",&m_wave_status[first]);
+ first++;
+ }
+}
+
+void NtupleDumperAlg::FillWaveBranches(const xAOD::WaveformHitContainer &wave) const {
+ for (auto hit : wave) {
+ if ((hit->hit_status()&2)==0) { // dont store secoondary hits as they can overwrite the primary hit
+ int ch=hit->channel();
+ m_wave_localtime[ch]=hit->localtime()+m_clock_phase;
+ m_wave_peak[ch]=hit->peak();
+ m_wave_width[ch]=hit->width();
+ m_wave_charge[ch]=hit->integral()/50;
+
+ m_wave_raw_peak[ch]=hit->raw_peak();
+ m_wave_raw_charge[ch]=hit->raw_integral()/50;
+ m_wave_baseline_mean[ch]=hit->baseline_mean();
+ m_wave_baseline_rms[ch]=hit->baseline_rms();
+ m_wave_status[ch]=hit->hit_status();
+ }
+ }
+}
+
+StatusCode NtupleDumperAlg::initialize()
+{
+ ATH_CHECK(m_truthEventContainer.initialize());
+ ATH_CHECK(m_truthParticleContainer.initialize());
+ ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_trackSegmentCollection.initialize());
+ ATH_CHECK(m_vetoNuContainer.initialize());
+ ATH_CHECK(m_vetoContainer.initialize());
+ ATH_CHECK(m_triggerContainer.initialize());
+ ATH_CHECK(m_preshowerContainer.initialize());
+ ATH_CHECK(m_ecalContainer.initialize());
+ ATH_CHECK(m_clusterContainer.initialize());
+ ATH_CHECK(m_simDataCollection.initialize());
+ ATH_CHECK(m_FaserTriggerData.initialize());
+ ATH_CHECK(m_ClockWaveformContainer.initialize());
+
+ ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
+ ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID"));
+ ATH_CHECK(detStore()->retrieve(m_vetoHelper, "VetoID"));
+ ATH_CHECK(detStore()->retrieve(m_triggerHelper, "TriggerID"));
+ ATH_CHECK(detStore()->retrieve(m_preshowerHelper, "PreshowerID"));
+ ATH_CHECK(detStore()->retrieve(m_ecalHelper, "EcalID"));
+
+ ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT"));
+ ATH_CHECK(m_extrapolationTool.retrieve());
+ ATH_CHECK(m_trackingGeometryTool.retrieve());
+
+ ATH_CHECK(m_spacePointContainerKey.initialize());
+
+ if (m_useFlukaWeights)
+ {
+ m_baseEventCrossSection = (m_flukaCrossSection * kfemtoBarnsPerMilliBarn)/m_flukaCollisions;
+ }
+ else if (m_useGenieWeights)
+ {
+ m_baseEventCrossSection = 1.0/m_genieLuminosity;
+ }
+ else
+ {
+ m_baseEventCrossSection = 1.0;
+ }
+
+ m_tree = new TTree("nt", "NtupleDumper tree");
+ m_tree->Branch("run", &m_run_number, "run/I");
+ m_tree->Branch("eventID", &m_event_number, "eventID/I");
+ m_tree->Branch("eventTime", &m_event_time, "eventTime/I");
+ m_tree->Branch("BCID", &m_bcid, "BCID/I");
+
+ m_tree->Branch("TBP", &m_tbp, "TBP/I");
+ m_tree->Branch("TAP", &m_tap, "TAP/I");
+ m_tree->Branch("inputBits", &m_inputBits, "inputBits/I");
+ m_tree->Branch("inputBitsNext", &m_inputBitsNext, "inputBitsNext/I");
+
+ addWaveBranches("VetoNu",2,4);
+ addWaveBranches("VetoSt1",2,6);
+ addWaveBranches("VetoSt2",1,14);
+ addWaveBranches("Timing",4,8);
+ addWaveBranches("Preshower",2,12);
+ addWaveBranches("Calo",4,0);
+ addBranch("Calo_total_charge", &m_calo_total);
+ addBranch("Calo_total_raw_charge", &m_calo_rawtotal);
+
+ addBranch("Calo0_Edep", &m_Calo0_Edep);
+ addBranch("Calo1_Edep", &m_Calo1_Edep);
+ addBranch("Calo2_Edep", &m_Calo2_Edep);
+ addBranch("Calo3_Edep", &m_Calo3_Edep);
+ addBranch("Calo_total_Edep", &m_Calo_Total_Edep);
+ addBranch("Preshower12_Edep", &m_Preshower12_Edep);
+ addBranch("Preshower13_Edep", &m_Preshower13_Edep);
+
+ addBranch("nClusters0",&m_station0Clusters);
+ addBranch("nClusters1",&m_station1Clusters);
+ addBranch("nClusters2",&m_station2Clusters);
+ addBranch("nClusters3",&m_station3Clusters);
+
+ addBranch("SpacePoints",&m_nspacepoints);
+ m_tree->Branch("SpacePoint_x", &m_spacepointX);
+ m_tree->Branch("SpacePoint_y", &m_spacepointY);
+ m_tree->Branch("SpacePoint_z", &m_spacepointZ);
+
+ addBranch("TrackSegments",&m_ntracksegs);
+ m_tree->Branch("TrackSegment_Chi2", &m_trackseg_Chi2);
+ m_tree->Branch("TrackSegment_nDoF", &m_trackseg_DoF);
+ m_tree->Branch("TrackSegment_x", &m_trackseg_x);
+ m_tree->Branch("TrackSegment_y", &m_trackseg_y);
+ m_tree->Branch("TrackSegment_z", &m_trackseg_z);
+ m_tree->Branch("TrackSegment_px", &m_trackseg_px);
+ m_tree->Branch("TrackSegment_py", &m_trackseg_py);
+ m_tree->Branch("TrackSegment_pz", &m_trackseg_pz);
+
+ m_tree->Branch("longTracks", &m_longTracks, "longTracks/I");
+ m_tree->Branch("Track_Chi2", &m_Chi2);
+ m_tree->Branch("Track_nDoF", &m_DoF);
+ m_tree->Branch("Track_x0", &m_xup);
+ m_tree->Branch("Track_y0", &m_yup);
+ m_tree->Branch("Track_z0", &m_zup);
+ m_tree->Branch("Track_px0", &m_pxup);
+ m_tree->Branch("Track_py0", &m_pyup);
+ m_tree->Branch("Track_pz0", &m_pzup);
+ m_tree->Branch("Track_p0", &m_pup);
+ m_tree->Branch("Track_x1", &m_xdown);
+ m_tree->Branch("Track_y1", &m_ydown);
+ m_tree->Branch("Track_z1", &m_zdown);
+ m_tree->Branch("Track_px1", &m_pxdown);
+ m_tree->Branch("Track_py1", &m_pydown);
+ m_tree->Branch("Track_pz1", &m_pzdown);
+ m_tree->Branch("Track_p1", &m_pdown);
+ m_tree->Branch("Track_charge", &m_charge);
+ m_tree->Branch("Track_nLayers", &m_nLayers);
+
+ m_tree->Branch("Track_InStation0",&m_nHit0);
+ m_tree->Branch("Track_InStation1",&m_nHit1);
+ m_tree->Branch("Track_InStation2",&m_nHit2);
+ m_tree->Branch("Track_InStation3",&m_nHit3);
+
+ m_tree->Branch("Track_X_atVetoNu", &m_xVetoNu);
+ m_tree->Branch("Track_Y_atVetoNu", &m_yVetoNu);
+ m_tree->Branch("Track_ThetaX_atVetoNu", &m_thetaxVetoNu);
+ m_tree->Branch("Track_ThetaY_atVetoNu", &m_thetayVetoNu);
+
+ m_tree->Branch("Track_X_atVetoStation1", &m_xVetoStation1);
+ m_tree->Branch("Track_Y_atVetoStation1", &m_yVetoStation1);
+ m_tree->Branch("Track_ThetaX_atVetoStation1", &m_thetaxVetoStation1);
+ m_tree->Branch("Track_ThetaY_atVetoStation1", &m_thetayVetoStation1);
+
+ m_tree->Branch("Track_X_atVetoStation2", &m_xVetoStation2);
+ m_tree->Branch("Track_Y_atVetoStation2", &m_yVetoStation2);
+ m_tree->Branch("Track_ThetaX_atVetoStation2", &m_thetaxVetoStation2);
+ m_tree->Branch("Track_ThetaY_atVetoStation2", &m_thetayVetoStation2);
+
+ m_tree->Branch("Track_X_atTrig", &m_xTrig);
+ m_tree->Branch("Track_Y_atTrig", &m_yTrig);
+ m_tree->Branch("Track_ThetaX_atTrig", &m_thetaxTrig);
+ m_tree->Branch("Track_ThetaY_atTrig", &m_thetayTrig);
+
+ m_tree->Branch("Track_X_atPreshower1", &m_xPreshower1);
+ m_tree->Branch("Track_Y_atPreshower1", &m_yPreshower1);
+ m_tree->Branch("Track_ThetaX_atPreshower1", &m_thetaxPreshower1);
+ m_tree->Branch("Track_ThetaY_atPreshower1", &m_thetayPreshower1);
+
+ m_tree->Branch("Track_X_atPreshower2", &m_xPreshower2);
+ m_tree->Branch("Track_Y_atPreshower2", &m_yPreshower2);
+ m_tree->Branch("Track_ThetaX_atPreshower2", &m_thetaxPreshower2);
+ m_tree->Branch("Track_ThetaY_atPreshower2", &m_thetayPreshower2);
+
+ m_tree->Branch("Track_X_atCalo", &m_xCalo);
+ m_tree->Branch("Track_Y_atCalo", &m_yCalo);
+ m_tree->Branch("Track_ThetaX_atCalo", &m_thetaxCalo);
+ m_tree->Branch("Track_ThetaY_atCalo", &m_thetayCalo);
+
+ m_tree->Branch("pTruthLepton", &m_truthLeptonMomentum, "pTruthLepton/D");
+ m_tree->Branch("truthBarcode", &m_truthBarcode, "truthBarcode/I");
+ m_tree->Branch("truthPdg", &m_truthPdg, "truthPdg/I");
+ m_tree->Branch("CrossSection", &m_crossSection, "crossSection/D");
+
+ ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
+
+ // Register histograms
+ m_HistRandomCharge[0] = new TH1F("hRandomCharge0", "Calo ch0 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[1] = new TH1F("hRandomCharge1", "Calo ch1 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[2] = new TH1F("hRandomCharge2", "Calo ch2 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[3] = new TH1F("hRandomCharge3", "Calo ch3 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[4] = new TH1F("hRandomCharge4", "VetoNu ch4 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[5] = new TH1F("hRandomCharge5", "VetoNu ch5 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[6] = new TH1F("hRandomCharge6", "Veto ch6 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[7] = new TH1F("hRandomCharge7", "Veto ch7 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[8] = new TH1F("hRandomCharge8", "Trig ch8 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[9] = new TH1F("hRandomCharge9", "Trig ch9 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[10] = new TH1F("hRandomCharge10", "Trig ch10 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[11] = new TH1F("hRandomCharge11", "Trig ch11 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[12] = new TH1F("hRandomCharge12", "Preshower ch12 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[13] = new TH1F("hRandomCharge13", "Preshower ch13 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+ m_HistRandomCharge[14] = new TH1F("hRandomCharge14", "Veto ch14 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
+
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge0", m_HistRandomCharge[0]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge1", m_HistRandomCharge[1]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge2", m_HistRandomCharge[2]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge3", m_HistRandomCharge[3]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge4", m_HistRandomCharge[4]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge5", m_HistRandomCharge[5]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge6", m_HistRandomCharge[6]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge7", m_HistRandomCharge[7]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge8", m_HistRandomCharge[8]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge9", m_HistRandomCharge[9]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge10", m_HistRandomCharge[10]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge11", m_HistRandomCharge[11]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge12", m_HistRandomCharge[12]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge13", m_HistRandomCharge[13]));
+ ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge14", m_HistRandomCharge[14]));
+
+ m_MIP_sim_Edep_calo = 0.0585; // MIP deposits 0.0585 GeV of energy in calo
+ m_MIP_sim_Edep_preshower = 0.004894; // MIP deposits 0.004894 GeV of energy in a preshower layer
+
+ if (m_doBlinding) {
+ ATH_MSG_INFO("Blinding will be enforced for real data.");
+ } else {
+ ATH_MSG_INFO("Blinding will NOT be enforced for real data.");
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
+{
+ clearTree();
+
+ // check if real data or simulation data
+ bool realData = true;
+ SG::ReadHandle<xAOD::TruthEventContainer> truthEventContainer { m_truthEventContainer, ctx };
+ if (truthEventContainer.isValid() && truthEventContainer->size() > 0)
+ {
+ realData = false;
+ }
+
+ // if real data, store charge in histograms from random events and only fill ntuple from coincidence events
+ if (realData) { //no trigger simulation yet
+ SG::ReadHandle<xAOD::FaserTriggerData> triggerData(m_FaserTriggerData, ctx);
+ m_tap=triggerData->tap();
+ if (m_tap==16) { // random trigger, store charge of scintillators in histograms
+ // Read in Waveform containers
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
+ ATH_CHECK(vetoNuContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
+ ATH_CHECK(vetoContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
+ ATH_CHECK(triggerContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
+ ATH_CHECK(preshowerContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
+ ATH_CHECK(ecalContainer.isValid());
+
+ if (vetoNuContainer.isValid()) {
+ for (auto hit : *vetoNuContainer) {
+ int ch=hit->channel();
+ m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
+ }
+ }
+ if (vetoContainer.isValid()) {
+ for (auto hit : *vetoContainer) {
+ int ch=hit->channel();
+ m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
+ }
+ }
+ if (triggerContainer.isValid()) {
+ for (auto hit : *triggerContainer) {
+ int ch=hit->channel();
+ m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
+ }
+ }
+ if (preshowerContainer.isValid()) {
+ for (auto hit : *preshowerContainer) {
+ int ch=hit->channel();
+ m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
+ }
+ }
+ if (ecalContainer.isValid()) {
+ for (auto hit : *ecalContainer) {
+ int ch=hit->channel();
+ m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
+ }
+ }
+
+ return StatusCode::SUCCESS; // finished with this event
+
+ } else if ( ((m_tap&8)==0) && (((m_tap&4)==0)||((m_tap&2)==0)) && (((m_tap&4)==0)||((m_tap&1)==0)) && (((m_tap&2)==0)||((m_tap&1)==0)) ) { // don't process events that don't trigger coincidence triggers: 1=calo, 2=veotnu|neto1|preshower, 4=TimingLayer, 8=(VetoNu|Veto2)&Preshower
+ return StatusCode::SUCCESS;
+ }
+ m_tbp=triggerData->tbp();
+ m_tap=triggerData->tap();
+ m_inputBits=triggerData->inputBits();
+ m_inputBitsNext=triggerData->inputBitsNextClk();
+ }
+
+ m_run_number = ctx.eventID().run_number();
+ m_event_number = ctx.eventID().event_number();
+ m_event_time = ctx.eventID().time_stamp();
+ m_bcid = ctx.eventID().bunch_crossing_id();
+
+ if (!realData) { // if simulation find MC cross section and primary lepton
+ // Work out effective cross section for MC
+ if (m_useFlukaWeights)
+ {
+ double flukaWeight = truthEventContainer->at(0)->weights()[0];
+ ATH_MSG_ALWAYS("Found fluka weight = " << flukaWeight);
+ m_crossSection = m_baseEventCrossSection * flukaWeight;
+ }
+ else if (m_useGenieWeights)
+ {
+ m_crossSection = m_baseEventCrossSection;
+ }
+ else
+ {
+ //ATH_MSG_WARNING("Monte carlo event with no weighting scheme specified. Setting crossSection (weight) to " << m_baseEventCrossSection << " fb.");
+ m_crossSection = m_baseEventCrossSection;
+ }
+
+ // Find the primary lepton (if any)
+ SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
+ if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0)
+ {
+ for (auto particle : *truthParticleContainer)
+ {
+ if ( particle->absPdgId() == 11 || particle->absPdgId() == 13 || particle->absPdgId() == 15 )
+ {
+ if (particle->status() == 1 && (particle->nParents() == 0 || particle->nParents() == 2) )
+ {
+ m_truthLeptonMomentum = particle->p4().P();
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ if (realData) { // correct waveform time with clock phase
+ SG::ReadHandle<xAOD::WaveformClock> clockHandle(m_ClockWaveformContainer, ctx);
+ ATH_CHECK(clockHandle.isValid());
+
+ if (clockHandle->phase() < -2.0) { // wrap around clock pahse so -pi goes to pi
+ m_clock_phase = ((clockHandle->phase() + 3.14159) / 3.14159) * 12.5;
+ } else {
+ m_clock_phase = (clockHandle->phase() / 3.14159) * 12.5;
+ }
+ }
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
+ ATH_CHECK(vetoNuContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
+ ATH_CHECK(vetoContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
+ ATH_CHECK(triggerContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
+ ATH_CHECK(preshowerContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
+ ATH_CHECK(ecalContainer.isValid());
+
+ FillWaveBranches(*vetoNuContainer);
+ FillWaveBranches(*vetoContainer);
+ FillWaveBranches(*triggerContainer);
+ FillWaveBranches(*preshowerContainer);
+ FillWaveBranches(*ecalContainer);
+
+ m_calo_total=m_wave_charge[0]+m_wave_charge[1]+m_wave_charge[2]+m_wave_charge[3];
+ m_calo_rawtotal=m_wave_raw_charge[0]+m_wave_raw_charge[1]+m_wave_raw_charge[2]+m_wave_raw_charge[3];
+
+ // do calibration of calo channels from pC to GeV deposited
+ if (m_CaloConfig == "High_gain") {
+ m_Calo0_Edep = (m_wave_charge[0] / 23.709) * m_MIP_sim_Edep_calo;
+ m_Calo1_Edep = (m_wave_charge[1] / 24.333) * m_MIP_sim_Edep_calo;
+ m_Calo2_Edep = (m_wave_charge[2] / 24.409) * m_MIP_sim_Edep_calo;
+ m_Calo3_Edep = (m_wave_charge[3] / 25.555) * m_MIP_sim_Edep_calo;
+ } else if (m_CaloConfig == "Low_gain") { // assume low gain calo
+ m_Calo0_Edep = (m_wave_charge[0] / 0.7909) * m_MIP_sim_Edep_calo;
+ m_Calo1_Edep = (m_wave_charge[1] / 0.8197) * m_MIP_sim_Edep_calo;
+ m_Calo2_Edep = (m_wave_charge[2] / 0.8256) * m_MIP_sim_Edep_calo;
+ m_Calo3_Edep = (m_wave_charge[3] / 0.8821) * m_MIP_sim_Edep_calo;
+ } else {
+ ATH_MSG_WARNING("Run config is neither High_gain nor Low_gain, it is " << m_CaloConfig << ", calo calibration will be zero");
+ }
+ m_Calo_Total_Edep = m_Calo0_Edep + m_Calo1_Edep + m_Calo2_Edep + m_Calo3_Edep;
+
+ // do calibration of preshower channels from pC to GeV deposited
+ m_Preshower12_Edep = (m_wave_charge[12] / 5.0) * m_MIP_sim_Edep_preshower; // 5 pC per MIP is rough measurement
+ m_Preshower13_Edep = (m_wave_charge[12] / 5.0) * m_MIP_sim_Edep_preshower;
+
+ if (realData && m_doBlinding) { // enforce blinding such that events with large calo signals are skipped and not in the output root file
+ if ((m_Calo_Total_Edep/0.155) > 10.0) { // only save events with a shower less than a 10 GeV e- (assume 10 GeV electron deposits 15.5% of their energy in calo)
+ return StatusCode::SUCCESS;
+ }
+ }
+
+ SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx };
+ ATH_CHECK(clusterContainer.isValid());
+
+ FaserActsGeometryContext faserGeometryContext = m_trackingGeometryTool->getNominalGeometryContext();
+ auto gctx = faserGeometryContext.context();
+
+ for (auto collection : *clusterContainer)
+ {
+ Identifier id = collection->identify();
+ int station = m_sctHelper->station(id);
+ int clusters = (int) collection->size();
+ switch (station)
+ {
+ case 0:
+ m_station0Clusters += clusters;
+ // following lines commented out depict how to access cluster position
+ //for (auto cluster : *collection) {
+ // if (cluster == nullptr) continue;
+ // auto pos = cluster->globalPosition();
+ // m_station0ClusterX.push_back(pos.x());
+ //}
+ break;
+ case 1:
+ m_station1Clusters += clusters;
+ break;
+ case 2:
+ m_station2Clusters += clusters;
+ break;
+ case 3:
+ m_station3Clusters += clusters;
+ break;
+ default:
+ ATH_MSG_FATAL("Unknown tracker station number " << station);
+ break;
+ }
+ }
+
+ SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey, ctx};
+ ATH_CHECK(spacePointContainer.isValid());
+ for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) {
+ m_nspacepoints += spacePointCollection->size();
+ for (const Tracker::FaserSCT_SpacePoint *spacePoint: *spacePointCollection) {
+ auto pos = spacePoint->globalPosition();
+ m_spacepointX.push_back(pos.x());
+ m_spacepointY.push_back(pos.y());
+ m_spacepointZ.push_back(pos.z());
+ }
+ }
+
+ SG::ReadHandle<TrackCollection> trackSegmentCollection {m_trackSegmentCollection, ctx};
+ ATH_CHECK(trackSegmentCollection.isValid());
+ for (const Trk::Track* trackSeg : *trackSegmentCollection) {
+ if (trackSeg == nullptr) continue;
+ m_ntracksegs += 1;
+ m_trackseg_Chi2.push_back(trackSeg->fitQuality()->chiSquared());
+ m_trackseg_DoF.push_back(trackSeg->fitQuality()->numberDoF());
+ auto SegParameters = trackSeg->trackParameters()->front();
+ const Amg::Vector3D SegPosition = SegParameters->position();
+ const Amg::Vector3D SegMomentum = SegParameters->momentum();
+ m_trackseg_x.push_back(SegPosition.x());
+ m_trackseg_y.push_back(SegPosition.y());
+ m_trackseg_z.push_back(SegPosition.z());
+ m_trackseg_px.push_back(SegMomentum.x());
+ m_trackseg_py.push_back(SegMomentum.y());
+ m_trackseg_pz.push_back(SegMomentum.z());
+ }
+
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+ ATH_CHECK(trackCollection.isValid());
+ const Trk::TrackParameters* candidateParameters {nullptr};
+ const Trk::TrackParameters* candidateDownParameters {nullptr};
+ for (const Trk::Track* track : *trackCollection)
+ {
+ if (track == nullptr) continue;
+ std::set<std::pair<int, int>> layerMap;
+ std::set<int> stationMap;
+
+ // Check for hit in the three downstream stations
+ for (auto measurement : *(track->measurementsOnTrack())) {
+ const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+ if (cluster != nullptr) {
+ Identifier id = cluster->identify();
+ int station = m_sctHelper->station(id);
+ int layer = m_sctHelper->layer(id);
+ stationMap.emplace(station);
+ layerMap.emplace(station, layer);
+ }
+ }
+ if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue;
+
+ int nLayers = std::count_if(layerMap.begin(), layerMap.end(), [](std::pair<int,int> p){return p.first != 0;});
+ const Trk::TrackParameters* upstreamParameters = track->trackParameters()->front();
+ const Trk::TrackParameters* downstreamParameters = track->trackParameters()->back();
+
+ if (candidateParameters == nullptr || upstreamParameters->momentum().mag() > candidateParameters->momentum().mag())
+ {
+ candidateParameters = upstreamParameters;
+ candidateDownParameters = downstreamParameters;
+ }
+
+ if ((candidateParameters == nullptr) || (candidateDownParameters == nullptr)) continue;
+
+ m_nLayers.push_back(nLayers);
+
+ m_Chi2.push_back(track->fitQuality()->chiSquared());
+ m_DoF.push_back(track->fitQuality()->numberDoF());
+
+ m_nHit0.push_back(stationMap.count(0));
+ m_nHit1.push_back(stationMap.count(1));
+ m_nHit2.push_back(stationMap.count(2));
+ m_nHit3.push_back(stationMap.count(3));
+
+ m_charge.push_back( (int) candidateParameters->charge() );
+
+ m_xup.push_back(candidateParameters->position().x());
+ m_yup.push_back(candidateParameters->position().y());
+ m_zup.push_back(candidateParameters->position().z());
+ m_pxup.push_back(candidateParameters->momentum().x());
+ m_pyup.push_back(candidateParameters->momentum().y());
+ m_pzup.push_back(candidateParameters->momentum().z());
+ m_pup.push_back(sqrt( pow(candidateParameters->momentum().x(),2) + pow(candidateParameters->momentum().y(),2) + pow(candidateParameters->momentum().z(),2) ));
+
+ m_xdown.push_back(candidateDownParameters->position().x());
+ m_ydown.push_back(candidateDownParameters->position().y());
+ m_zdown.push_back(candidateDownParameters->position().z());
+ m_pxdown.push_back(candidateDownParameters->momentum().x());
+ m_pydown.push_back(candidateDownParameters->momentum().y());
+ m_pzdown.push_back(candidateDownParameters->momentum().z());
+ m_pdown.push_back(sqrt( pow(candidateDownParameters->momentum().x(),2) + pow(candidateDownParameters->momentum().y(),2) + pow(candidateDownParameters->momentum().z(),2) ));
+
+ // fill extrapolation vectors with filler values that get changed iif the track extrapolation succeeds
+ m_xVetoNu.push_back(-10000);
+ m_yVetoNu.push_back(-10000);
+ m_thetaxVetoNu.push_back(-10000);
+ m_thetayVetoNu.push_back(-10000);
+ m_xVetoStation1.push_back(-10000);
+ m_yVetoStation1.push_back(-10000);
+ m_thetaxVetoStation1.push_back(-10000);
+ m_thetayVetoStation1.push_back(-10000);
+ m_xVetoStation2.push_back(-10000);
+ m_yVetoStation2.push_back(-10000);
+ m_thetaxVetoStation2.push_back(-10000);
+ m_thetayVetoStation2.push_back(-10000);
+ m_xTrig.push_back(-10000);
+ m_yTrig.push_back(-10000);
+ m_thetaxTrig.push_back(-10000);
+ m_thetayTrig.push_back(-10000);
+ m_xPreshower1.push_back(-10000);
+ m_yPreshower1.push_back(-10000);
+ m_thetaxPreshower1.push_back(-10000);
+ m_thetayPreshower1.push_back(-10000);
+ m_xPreshower2.push_back(-10000);
+ m_yPreshower2.push_back(-10000);
+ m_thetaxPreshower2.push_back(-10000);
+ m_thetayPreshower2.push_back(-10000);
+ m_xCalo.push_back(-10000);
+ m_yCalo.push_back(-10000);
+ m_thetaxCalo.push_back(-10000);
+ m_thetayCalo.push_back(-10000);
+
+ // extrapolate track from IFT
+ if (stationMap.count(0) > 0) { // extrapolation crashes if the track does not start in the IFT, as it is too far away to extrapolate
+ Amg::Vector3D position = candidateParameters->position();
+ Amg::Vector3D momentum = candidateParameters->momentum();
+ Acts::BoundVector params = Acts::BoundVector::Zero();
+ params[Acts::eBoundLoc0] = -position.y();
+ params[Acts::eBoundLoc1] = position.x();
+ params[Acts::eBoundPhi] = momentum.phi();
+ params[Acts::eBoundTheta] = momentum.theta();
+ params[Acts::eBoundQOverP] = candidateParameters->charge() / momentum.mag();
+ params[Acts::eBoundTime] = 0;
+ auto startSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position.z()), Acts::Vector3(0, 0, 1));
+ Acts::BoundTrackParameters startParameters(std::move(startSurface), params, candidateParameters->charge());
+
+ auto targetSurface_VetoNu = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -3112.0), Acts::Vector3(0, 0, 1)); // -3112 mm is z position of VetoNu planes touching
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_VetoNu =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_VetoNu, Acts::backward);
+ if (targetParameters_VetoNu != nullptr) {
+ auto targetPosition_VetoNu = targetParameters_VetoNu->position(gctx);
+ auto targetMomentum_VetoNu = targetParameters_VetoNu->momentum();
+ m_xVetoNu[m_longTracks] = targetPosition_VetoNu.x();
+ m_yVetoNu[m_longTracks] = targetPosition_VetoNu.y();
+ m_thetaxVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[0]/targetMomentum_VetoNu[2]);
+ m_thetayVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[1]/targetMomentum_VetoNu[2]);
+ } else {
+ ATH_MSG_INFO("vetoNu null targetParameters");
+ }
+
+ auto targetSurface_Veto1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1769.65), Acts::Vector3(0, 0, 1)); // -1769.65 mm is z position of center of operational layer in Veto station 1
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto1 =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Veto1, Acts::forward);
+ if (targetParameters_Veto1 != nullptr) {
+ auto targetPosition_Veto1 = targetParameters_Veto1->position(gctx);
+ auto targetMomentum_Veto1 = targetParameters_Veto1->momentum();
+ m_xVetoStation1[m_longTracks] = targetPosition_Veto1.x();
+ m_yVetoStation1[m_longTracks] = targetPosition_Veto1.y();
+ m_thetaxVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[0]/targetMomentum_Veto1[2]);
+ m_thetayVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[1]/targetMomentum_Veto1[2]);
+ } else {
+ ATH_MSG_INFO("veto1 null targetParameters");
+ }
+
+ auto targetSurface_Veto2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1609.65), Acts::Vector3(0, 0, 1)); // -1609.65 mm is z position of where planes touch in Veto station 2
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto2 =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Veto2, Acts::forward);
+ if (targetParameters_Veto2 != nullptr) {
+ auto targetPosition_Veto2 = targetParameters_Veto2->position(gctx);
+ auto targetMomentum_Veto2 = targetParameters_Veto2->momentum();
+ m_xVetoStation2[m_longTracks] = targetPosition_Veto2.x();
+ m_yVetoStation2[m_longTracks] = targetPosition_Veto2.y();
+ m_thetaxVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[0]/targetMomentum_Veto2[2]);
+ m_thetayVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[1]/targetMomentum_Veto2[2]);
+ } else {
+ ATH_MSG_INFO("veto2 null targetParameters");
+ }
+
+ auto targetSurface_Trig = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 0.0), Acts::Vector3(0, 0, 1)); // 0 mm is z position of Trig planes overlapping
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Trig =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Trig, Acts::forward); // must extrapolate forward to trig plane if track starts in IFT
+ if (targetParameters_Trig != nullptr) {
+ auto targetPosition_Trig = targetParameters_Trig->position(gctx);
+ auto targetMomentum_Trig = targetParameters_Trig->momentum();
+ m_xTrig[m_longTracks] = targetPosition_Trig.x();
+ m_yTrig[m_longTracks] = targetPosition_Trig.y();
+ m_thetaxTrig[m_longTracks] = atan(targetMomentum_Trig[0]/targetMomentum_Trig[2]);
+ m_thetayTrig[m_longTracks] = atan(targetMomentum_Trig[1]/targetMomentum_Trig[2]);
+ } else {
+ ATH_MSG_INFO("Trig null targetParameters");
+ }
+
+ }
+
+ // extrapolate track from tracking station 3
+ if (stationMap.count(3) > 0) { // extrapolation crashes if the track does not end in the Station 3, as it is too far away to extrapolate
+ Amg::Vector3D positionDown = candidateDownParameters->position();
+ Amg::Vector3D momentumDown = candidateDownParameters->momentum();
+ Acts::BoundVector paramsDown = Acts::BoundVector::Zero();
+ paramsDown[Acts::eBoundLoc0] = -positionDown.y();
+ paramsDown[Acts::eBoundLoc1] = positionDown.x();
+ paramsDown[Acts::eBoundPhi] = momentumDown.phi();
+ paramsDown[Acts::eBoundTheta] = momentumDown.theta();
+ paramsDown[Acts::eBoundQOverP] = candidateDownParameters->charge() / momentumDown.mag();
+ paramsDown[Acts::eBoundTime] = 0;
+ auto startSurfaceDown = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, positionDown.z()), Acts::Vector3(0, 0, 1));
+ Acts::BoundTrackParameters startParametersDown(std::move(startSurfaceDown), paramsDown, candidateDownParameters->charge());
+
+ auto targetSurface_Preshower1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2582.68), Acts::Vector3(0, 0, 1)); // 2582.68 mm is z position of center of upstream preshower layer
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower1 =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Preshower1, Acts::forward);
+ if (targetParameters_Preshower1 != nullptr) {
+ auto targetPosition_Preshower1 = targetParameters_Preshower1->position(gctx);
+ auto targetMomentum_Preshower1 = targetParameters_Preshower1->momentum();
+ m_xPreshower1[m_longTracks] = targetPosition_Preshower1.x();
+ m_yPreshower1[m_longTracks] = targetPosition_Preshower1.y();
+ m_thetaxPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[0]/targetMomentum_Preshower1[2]);
+ m_thetayPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[1]/targetMomentum_Preshower1[2]);
+ } else {
+ ATH_MSG_INFO("Preshower1 null targetParameters");
+ }
+
+ auto targetSurface_Preshower2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2657.68), Acts::Vector3(0, 0, 1)); // 2657.68 mm is z position of center of downstream preshower layer
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower2 =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Preshower2, Acts::forward);
+ if (targetParameters_Preshower2 != nullptr) {
+ auto targetPosition_Preshower2 = targetParameters_Preshower2->position(gctx);
+ auto targetMomentum_Preshower2 = targetParameters_Preshower2->momentum();
+ m_xPreshower2[m_longTracks] = targetPosition_Preshower2.x();
+ m_yPreshower2[m_longTracks] = targetPosition_Preshower2.y();
+ m_thetaxPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[0]/targetMomentum_Preshower2[2]);
+ m_thetayPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[1]/targetMomentum_Preshower2[2]);
+ } else {
+ ATH_MSG_INFO("Preshower2 null targetParameters");
+ }
+
+ auto targetSurface_Calo = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2760.0), Acts::Vector3(0, 0, 1)); // 2760 mm is estimated z position of calorimeter face
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Calo =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Calo, Acts::forward);
+ if (targetParameters_Calo != nullptr) {
+ auto targetPosition_Calo = targetParameters_Calo->position(gctx);
+ auto targetMomentum_Calo = targetParameters_Calo->momentum();
+ m_xCalo[m_longTracks] = targetPosition_Calo.x();
+ m_yCalo[m_longTracks] = targetPosition_Calo.y();
+ m_thetaxCalo[m_longTracks] = atan(targetMomentum_Calo[0]/targetMomentum_Calo[2]) ;
+ m_thetayCalo[m_longTracks] = atan(targetMomentum_Calo[1]/targetMomentum_Calo[2]) ;
+ } else {
+ ATH_MSG_INFO("Calo null targetParameters");
+ }
+ }
+
+ m_longTracks++;
+ }
+
+ /*
+ // Here we apply the signal selection
+ // Very simple/unrealistic to start
+ if (m_vetoUpstream == 0 || m_vetoDownstream == 0 ||
+ m_triggerTotal == 0 ||
+ m_preshower0 == 0 || m_preshower1 == 0 ||
+ // m_ecalTotal == 0 ||
+ candidateParameters == nullptr)
+ return StatusCode::SUCCESS;
+ */
+ m_tree->Fill();
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode NtupleDumperAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
+
+bool NtupleDumperAlg::waveformHitOK(const xAOD::WaveformHit* hit) const
+{
+ if (hit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED) || hit->status_bit(xAOD::WaveformStatus::SECONDARY)) return false;
+ return true;
+}
+
+void
+NtupleDumperAlg::clearTree() const
+{
+ m_run_number = 0;
+ m_event_number = 0;
+ m_event_time = 0;
+ m_bcid = 0;
+
+ m_tbp=0;
+ m_tap=0;
+ m_inputBits=0;
+ m_inputBitsNext=0;
+
+ for(int ii=0;ii<15;ii++) {
+ m_wave_localtime[ii]=0;
+ m_wave_peak[ii]=0;
+ m_wave_width[ii]=0;
+ m_wave_charge[ii]=0;
+
+ m_wave_raw_peak[ii]=0;
+ m_wave_raw_charge[ii]=0;
+ m_wave_baseline_mean[ii]=0;
+ m_wave_baseline_rms[ii]=0;
+ m_wave_status[ii]=0;
+ }
+
+ m_calo_total=0;
+ m_calo_rawtotal=0;
+
+ m_Calo0_Edep=0;
+ m_Calo1_Edep=0;
+ m_Calo2_Edep=0;
+ m_Calo3_Edep=0;
+ m_Calo_Total_Edep=0;
+ m_Preshower12_Edep=0;
+ m_Preshower13_Edep=0;
+
+ m_clock_phase=0;
+
+ m_station0Clusters = 0;
+ m_station1Clusters = 0;
+ m_station2Clusters = 0;
+ m_station3Clusters = 0;
+ m_crossSection = 0;
+
+ m_nspacepoints = 0;
+ m_spacepointX.clear();
+ m_spacepointY.clear();
+ m_spacepointZ.clear();
+
+ m_ntracksegs = 0;
+ m_trackseg_Chi2.clear();
+ m_trackseg_DoF.clear();
+ m_trackseg_x.clear();
+ m_trackseg_y.clear();
+ m_trackseg_z.clear();
+ m_trackseg_px.clear();
+ m_trackseg_py.clear();
+ m_trackseg_pz.clear();
+
+ m_xup.clear();
+ m_yup.clear();
+ m_zup.clear();
+ m_pxup.clear();
+ m_pyup.clear();
+ m_pzup.clear();
+ m_pup.clear();
+
+ m_xdown.clear();
+ m_ydown.clear();
+ m_zdown.clear();
+ m_pxdown.clear();
+ m_pydown.clear();
+ m_pzdown.clear();
+ m_pdown.clear();
+
+ m_Chi2.clear();
+ m_DoF.clear();
+ m_charge.clear();
+ m_nLayers.clear();
+ m_longTracks = 0;
+
+ m_nHit0.clear();
+ m_nHit1.clear();
+ m_nHit2.clear();
+ m_nHit3.clear();
+
+ m_xVetoNu.clear();
+ m_yVetoNu.clear();
+ m_thetaxVetoNu.clear();
+ m_thetayVetoNu.clear();
+
+ m_xVetoStation1.clear();
+ m_yVetoStation1.clear();
+ m_thetaxVetoStation1.clear();
+ m_thetayVetoStation1.clear();
+
+ m_xVetoStation2.clear();
+ m_yVetoStation2.clear();
+ m_thetaxVetoStation2.clear();
+ m_thetayVetoStation2.clear();
+
+ m_xTrig.clear();
+ m_yTrig.clear();
+ m_thetaxTrig.clear();
+ m_thetayTrig.clear();
+
+ m_xPreshower1.clear();
+ m_yPreshower1.clear();
+ m_thetaxPreshower1.clear();
+ m_thetayPreshower1.clear();
+
+ m_xPreshower2.clear();
+ m_yPreshower2.clear();
+ m_thetaxPreshower2.clear();
+ m_thetayPreshower2.clear();
+
+ m_xCalo.clear();
+ m_yCalo.clear();
+ m_thetaxCalo.clear();
+ m_thetayCalo.clear();
+
+ m_truthLeptonMomentum = 0;
+ m_truthBarcode = 0;
+ m_truthPdg = 0;
+}
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..1b26caadeb7290c34a199ef0e941f9203b600bba
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -0,0 +1,217 @@
+#ifndef NTUPLEDUMPER_NTUPLEDUMPERALG_H
+#define NTUPLEDUMPER_NTUPLEDUMPERALG_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthHistogramming.h"
+#include "TrkTrack/TrackCollection.h"
+#include "xAODFaserTrigger/FaserTriggerData.h"
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "xAODFaserWaveform/WaveformHit.h"
+#include "xAODFaserWaveform/WaveformClock.h"
+#include "xAODTruth/TruthEventContainer.h"
+#include "xAODTruth/TruthParticleContainer.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h"
+#include "TrackerSimData/TrackerSimDataCollection.h"
+#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+
+#include <vector>
+
+class TTree;
+class TH1;
+class FaserSCT_ID;
+class VetoNuID;
+class VetoID;
+class TriggerID;
+class PreshowerID;
+class EcalID;
+namespace TrackerDD
+{
+ class SCT_DetectorManager;
+}
+
+class NtupleDumperAlg : public AthReentrantAlgorithm, AthHistogramming {
+public:
+ NtupleDumperAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~NtupleDumperAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+ const ServiceHandle <ITHistSvc> &histSvc() const;
+
+private:
+
+ bool waveformHitOK(const xAOD::WaveformHit* hit) const;
+ void clearTree() const;
+ void addBranch(const std::string &name,float* var);
+ void addBranch(const std::string &name,unsigned int* var);
+ void addWaveBranches(const std::string &name, int nchannels, int first);
+ void FillWaveBranches(const xAOD::WaveformHitContainer &wave) const;
+
+ ServiceHandle <ITHistSvc> m_histSvc;
+
+ SG::ReadHandleKey<xAOD::TruthEventContainer> m_truthEventContainer { this, "EventContainer", "TruthEvents", "Truth event container name." };
+ SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthParticleContainer { this, "ParticleContainer", "TruthParticles", "Truth particle container name." };
+ SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollection {this, "TrackerSimDataCollection", "SCT_SDO_Map"};
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollection", "Input track collection name" };
+ SG::ReadHandleKey<TrackCollection> m_trackSegmentCollection {this, "TrackSegmentCollection", "SegmentFit", "Input track segment collection name"};
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_vetoNuContainer { this, "VetoNuContainer", "VetoNuWaveformHits", "VetoNu hit container name" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_vetoContainer { this, "VetoContainer", "VetoWaveformHits", "Veto hit container name" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_triggerContainer { this, "TriggerContainer", "TriggerWaveformHits", "Trigger hit container name" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_preshowerContainer { this, "PreshowerContainer", "PreshowerWaveformHits", "Preshower hit container name" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_ecalContainer { this, "EcalContainer", "CaloWaveformHits", "Ecal hit container name" };
+ SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainer { this, "ClusterContainer", "SCT_ClusterContainer", "Tracker cluster container name" };
+ SG::ReadHandleKey<FaserSCT_SpacePointContainer> m_spacePointContainerKey { this, "SpacePoints", "SCT_SpacePointContainer", "space point container"};
+
+ SG::ReadHandleKey<xAOD::FaserTriggerData> m_FaserTriggerData { this, "FaserTriggerDataKey", "FaserTriggerData", "ReadHandleKey for xAOD::FaserTriggerData"};
+ SG::ReadHandleKey<xAOD::WaveformClock> m_ClockWaveformContainer { this, "WaveformClockKey", "WaveformClock", "ReadHandleKey for ClockWaveforms Container"};
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
+ ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
+
+ const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
+
+ const FaserSCT_ID* m_sctHelper;
+ const VetoNuID* m_vetoNuHelper;
+ const VetoID* m_vetoHelper;
+ const TriggerID* m_triggerHelper;
+ const PreshowerID* m_preshowerHelper;
+ const EcalID* m_ecalHelper;
+
+ StringProperty m_CaloConfig { this, "CaloConfig", "Low_gain", "Configuration found at http://aagaard.web.cern.ch/aagaard/FASERruns.html (spaces replaced with '_')" };
+ BooleanProperty m_doBlinding { this, "DoBlinding", true, "Blinding will not output events with Calo signal > 10 GeV e-" };
+ BooleanProperty m_useFlukaWeights { this, "UseFlukaWeights", false, "Flag to weight events according to value stored in HepMC::GenEvent" };
+ BooleanProperty m_useGenieWeights { this, "UseGenieWeights", false, "Flag to weight events according to Genie luminosity" };
+ IntegerProperty m_flukaCollisions { this, "FlukaCollisions", 137130000, "Number of proton-proton collisions in FLUKA sample." };
+ DoubleProperty m_flukaCrossSection { this, "FlukaCrossSection", 80.0, "Fluka p-p inelastic cross-section in millibarns." };
+ DoubleProperty m_genieLuminosity { this, "GenieLuminosity", 150.0, "Genie luminosity in inverse fb." };
+
+ double m_baseEventCrossSection {1.0};
+ const double kfemtoBarnsPerMilliBarn {1.0e12};
+
+ mutable TTree* m_tree;
+
+ mutable TH1* m_HistRandomCharge[15];
+
+ mutable unsigned int m_run_number;
+ mutable unsigned int m_event_number;
+ mutable unsigned int m_event_time;
+ mutable unsigned int m_bcid;
+
+ mutable unsigned int m_tbp;
+ mutable unsigned int m_tap;
+ mutable unsigned int m_inputBits;
+ mutable unsigned int m_inputBitsNext;
+
+ mutable float m_wave_localtime[15];
+ mutable float m_wave_peak[15];
+ mutable float m_wave_width[15];
+ mutable float m_wave_charge[15];
+
+ mutable float m_wave_raw_peak[15];
+ mutable float m_wave_raw_charge[15];
+ mutable float m_wave_baseline_mean[15];
+ mutable float m_wave_baseline_rms[15];
+ mutable unsigned int m_wave_status[15];
+
+ mutable float m_calo_total;
+ mutable float m_calo_rawtotal;
+
+ mutable float m_Calo0_Edep;
+ mutable float m_Calo1_Edep;
+ mutable float m_Calo2_Edep;
+ mutable float m_Calo3_Edep;
+ mutable float m_Calo_Total_Edep;
+ mutable float m_Preshower12_Edep;
+ mutable float m_Preshower13_Edep;
+
+ mutable float m_MIP_sim_Edep_calo;
+ mutable float m_MIP_sim_Edep_preshower;
+
+ mutable float m_clock_phase;
+
+ mutable unsigned int m_station0Clusters;
+ mutable unsigned int m_station1Clusters;
+ mutable unsigned int m_station2Clusters;
+ mutable unsigned int m_station3Clusters;
+
+ mutable unsigned int m_nspacepoints;
+ mutable std::vector<double> m_spacepointX;
+ mutable std::vector<double> m_spacepointY;
+ mutable std::vector<double> m_spacepointZ;
+
+ mutable unsigned int m_ntracksegs;
+ mutable std::vector<double> m_trackseg_Chi2;
+ mutable std::vector<double> m_trackseg_DoF;
+ mutable std::vector<double> m_trackseg_x;
+ mutable std::vector<double> m_trackseg_y;
+ mutable std::vector<double> m_trackseg_z;
+ mutable std::vector<double> m_trackseg_px;
+ mutable std::vector<double> m_trackseg_py;
+ mutable std::vector<double> m_trackseg_pz;
+
+ mutable int m_longTracks;
+ mutable std::vector<double> m_Chi2;
+ mutable std::vector<double> m_DoF;
+ mutable std::vector<double> m_xup;
+ mutable std::vector<double> m_yup;
+ mutable std::vector<double> m_zup;
+ mutable std::vector<double> m_pxup;
+ mutable std::vector<double> m_pyup;
+ mutable std::vector<double> m_pzup;
+ mutable std::vector<double> m_pup;
+ mutable std::vector<double> m_xdown;
+ mutable std::vector<double> m_ydown;
+ mutable std::vector<double> m_zdown;
+ mutable std::vector<double> m_pxdown;
+ mutable std::vector<double> m_pydown;
+ mutable std::vector<double> m_pzdown;
+ mutable std::vector<double> m_pdown;
+ mutable std::vector<int> m_charge;
+ mutable std::vector<unsigned int> m_nLayers;
+ mutable std::vector<unsigned int> m_nHit0;
+ mutable std::vector<unsigned int> m_nHit1;
+ mutable std::vector<unsigned int> m_nHit2;
+ mutable std::vector<unsigned int> m_nHit3;
+ mutable std::vector<double> m_xVetoNu;
+ mutable std::vector<double> m_yVetoNu;
+ mutable std::vector<double> m_thetaxVetoNu;
+ mutable std::vector<double> m_thetayVetoNu;
+ mutable std::vector<double> m_xVetoStation1;
+ mutable std::vector<double> m_yVetoStation1;
+ mutable std::vector<double> m_thetaxVetoStation1;
+ mutable std::vector<double> m_thetayVetoStation1;
+ mutable std::vector<double> m_xVetoStation2;
+ mutable std::vector<double> m_yVetoStation2;
+ mutable std::vector<double> m_thetaxVetoStation2;
+ mutable std::vector<double> m_thetayVetoStation2;
+ mutable std::vector<double> m_xTrig;
+ mutable std::vector<double> m_yTrig;
+ mutable std::vector<double> m_thetaxTrig;
+ mutable std::vector<double> m_thetayTrig;
+ mutable std::vector<double> m_xPreshower1;
+ mutable std::vector<double> m_yPreshower1;
+ mutable std::vector<double> m_thetaxPreshower1;
+ mutable std::vector<double> m_thetayPreshower1;
+ mutable std::vector<double> m_xPreshower2;
+ mutable std::vector<double> m_yPreshower2;
+ mutable std::vector<double> m_thetaxPreshower2;
+ mutable std::vector<double> m_thetayPreshower2;
+ mutable std::vector<double> m_xCalo;
+ mutable std::vector<double> m_yCalo;
+ mutable std::vector<double> m_thetaxCalo;
+ mutable std::vector<double> m_thetayCalo;
+
+ mutable double m_truthLeptonMomentum;
+ mutable int m_truthBarcode;
+ mutable int m_truthPdg;
+ mutable double m_crossSection;
+
+};
+
+inline const ServiceHandle <ITHistSvc> &NtupleDumperAlg::histSvc() const {
+ return m_histSvc;
+}
+
+#endif // NTUPLEDUMPER_NTUPLEDUMPERALG_H
diff --git a/PhysicsAnalysis/NtupleDumper/src/component/NtupleDumper_entries.cxx b/PhysicsAnalysis/NtupleDumper/src/component/NtupleDumper_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..15dc244e170829e01af0d848e3cac55abbc83619
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/src/component/NtupleDumper_entries.cxx
@@ -0,0 +1,3 @@
+#include "../NtupleDumperAlg.h"
+
+DECLARE_COMPONENT(NtupleDumperAlg)
diff --git a/Simulation/G4Faser/G4FaserAlg/test/runGeantinoScan.py b/Simulation/G4Faser/G4FaserAlg/test/runGeantinoScan.py
index 7e0a93719a520de7bb819fe270f8666a00c42f5c..46021eb4bc7a37eefbcdfd0b14ab72349c563364 100644
--- a/Simulation/G4Faser/G4FaserAlg/test/runGeantinoScan.py
+++ b/Simulation/G4Faser/G4FaserAlg/test/runGeantinoScan.py
@@ -19,41 +19,34 @@ if __name__ == "__main__":
from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
from G4FaserAlg.G4FaserAlgConfigNew import G4FaserAlgCfg
- from G4FaserAlg.G4FaserAlgConfigNew import G4FaserMaterialStepRecorderAlgCfg
from G4FaserServices.G4FaserServicesConfigNew import G4GeometryNotifierSvcCfg
+ from G4FaserServices.G4FaserUserActionConfigNew import UserActionMaterialStepRecorderSvcCfg
#
# Set up logging and new style config
#
log.setLevel(VERBOSE)
Configurable.configurableRun3Behavior = True
-#
-# Input settings (Generator file)
-#
-# from AthenaConfiguration.TestDefaults import defaultTestFiles
-# ConfigFlags.Input.Files = defaultTestFiles.EVNT
-#
-# Alternatively, these must ALL be explicitly set to run without an input file
-# (if missing, it will try to read metadata from a non-existent file and crash)
-#
+ from AthenaConfiguration.Enums import ProductionStep
+ ConfigFlags.Common.ProductionStep = ProductionStep.Simulation
+ ConfigFlags.Sim.ReleaseGeoModel = False
ConfigFlags.Input.Files = [""]
ConfigFlags.Input.isMC = True
- ConfigFlags.Input.RunNumber = 12345
- ConfigFlags.Input.Collections = [""]
- ConfigFlags.Input.ProjectName = "mc19"
+ ConfigFlags.Input.RunNumber = [12345]
+ ConfigFlags.Input.OverrideRunNumber = True
+ ConfigFlags.Input.LumiBlockNumber = [1]
+ Configurable.configurableRun3Behavior = 1
ConfigFlags.Common.isOnline = False
ConfigFlags.Beam.Type = "collisions"
ConfigFlags.Beam.Energy = 7*TeV # Informational, does not affect simulation
- ConfigFlags.GeoModel.FaserVersion = "FASER-01" # Always needed
- ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" # Always needed; must match FaserVersion
-# Workaround for bug/missing flag; unimportant otherwise
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Always needed
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion
ConfigFlags.addFlag("Input.InitialTimeStamp", 0)
-# Workaround to avoid problematic ISF code
- ConfigFlags.GeoModel.Layout = "Development"
#
# Output settings
#
ConfigFlags.Output.HITSFileName = "MaterialStepCollection.root"
ConfigFlags.GeoModel.GeoExportFile = "faserGeo.db" # Optional dump of geometry for browsing in vp1light
+ ConfigFlags.GeoModel.Align.Dynamic = False
#
# Geometry-related settings
# Do not change!
@@ -61,8 +54,6 @@ if __name__ == "__main__":
detectors = ['Veto', 'Trigger', 'Preshower', 'FaserSCT', 'Dipole', 'Ecal']
from CalypsoConfiguration.DetectorConfigFlags import setupDetectorsFromList
setupDetectorsFromList(ConfigFlags, detectors, toggle_geometry=True)
- ConfigFlags.GeoModel.Align.Dynamic = False
- ConfigFlags.Sim.ReleaseGeoModel = False
#
# All flags should be set before calling lock
#
@@ -79,9 +70,9 @@ if __name__ == "__main__":
pg.McEventKey = "BeamTruthEvent"
pg.randomSeed = 123456
pg.sampler.pid = 999
- pg.sampler.mom = PG.EThetaMPhiSampler(energy=1*TeV, theta=[0, pi/20], phi=[0, 2*pi])
- pg.sampler.pos = PG.PosSampler(x=[-5, 5], y=[-5, 5], z=-2100.0, t=0.0)
- acc.addEventAlgo(pg, "AthBeginSeq") # to run *before* G4
+ pg.sampler.mom = PG.EThetaMPhiSampler(energy=1*TeV, theta=[0, pi/200], phi=[0, 2*pi])
+ pg.sampler.pos = PG.PosSampler(x=[-150, 150], y=[-150, 150], z=-2100.0, t=0.0)
+ acc.addEventAlgo(pg, "AthBeginSeq", primary = True) # to run *before* G4
#
# Only one of these two should be used in a given job
# (MCEventSelectorCfg for generating events with no input file,
@@ -107,7 +98,8 @@ if __name__ == "__main__":
# Here is the configuration of the Geant4 pieces
#
acc.merge(FaserGeometryCfg(ConfigFlags))
- acc.merge(G4FaserMaterialStepRecorderAlgCfg(ConfigFlags))
+ acc.merge(UserActionMaterialStepRecorderSvcCfg(ConfigFlags))
+ acc.merge(G4FaserAlgCfg(ConfigFlags))
acc.addService(G4GeometryNotifierSvcCfg(ConfigFlags, ActivateLVNotifier=True))
#
# Verbosity
@@ -125,4 +117,4 @@ if __name__ == "__main__":
#
# Execute and finish
#
- sys.exit(int(acc.run(maxEvents=20000).isFailure()))
+ sys.exit(int(acc.run(maxEvents=1000000).isFailure()))
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 4bf06b76f66ba0af57e64646cebc12b18b978612..829b0df04de5615d1509b7dbe741c410ef6b41df 100644
--- a/Simulation/ISF/ISF_Core/FaserISF_Services/python/FaserISF_ServicesConfigNew.py
+++ b/Simulation/ISF/ISF_Core/FaserISF_Services/python/FaserISF_ServicesConfigNew.py
@@ -11,7 +11,7 @@ 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 TruthStrategyGroupCfg
+from FaserISF_HepMC_Tools.FaserISF_HepMC_ToolsConfigNew import TruthStrategyGroupCfg, TrenchStrategyGroupCfg
ISF__FaserTruthSvc, ISF__FaserGeoIDSvc, ISF__FaserInputConverter = CompFactory.getComps("ISF::FaserTruthSvc","ISF::FaserGeoIDSvc","ISF::FaserInputConverter")
@@ -65,8 +65,10 @@ def FaserTruthServiceCfg(ConfigFlags, name="FaserISF_TruthService", **kwargs):
# acc2= FaserDipoleTruthStrategyCfg(ConfigFlags)
# kwargs.setdefault('TruthStrategies',[result.popToolsAndMerge(acc), result.popToolsAndMerge(acc2)])
acc = TruthStrategyGroupCfg(ConfigFlags)
- kwargs.setdefault('TruthStrategies', [result.popToolsAndMerge(acc)])
-
+ # FaserNu hack
+ acc2 = TrenchStrategyGroupCfg(ConfigFlags)
+ kwargs.setdefault('TruthStrategies', [result.popToolsAndMerge(acc), result.popToolsAndMerge(acc2)])
+
kwargs.setdefault('SkipIfNoChildren', True)
kwargs.setdefault('SkipIfNoParentBarcode', True)
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 4f739f6c0ca89afab6020eab6f42e096e1f53545..21a86e2026e8a13dee346d43616b677b601f1f7a 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
@@ -162,7 +162,7 @@ def TruthStrategyGroupCfg(ConfigFlags, name="ISF_MCTruthStrategyGroupID", **kwar
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("VertexTypes", [3, 14, 15, 4, 5, 6, 7, 2, 12, 13, 111, 121, 131, 132, 141, 151, 152, 161, 310]) # EM *and* nuclear
kwargs.setdefault("VertexTypeRangeLow" , 201) # All kinds of decay processes
kwargs.setdefault("VertexTypeRangeHigh" , 298) # ...
kwargs.setdefault("Regions", [FaserRegion.fFaserNeutrino,
@@ -173,6 +173,23 @@ def TruthStrategyGroupCfg(ConfigFlags, name="ISF_MCTruthStrategyGroupID", **kwar
result.setPrivateTools(CompFactory.ISF.FaserTruthStrategy(name, **kwargs))
return result
+# FaserNu hack
+def TrenchStrategyGroupCfg(ConfigFlags, name="ISF_TrenchStrategyGroupID", **kwargs):
+
+ import ROOT, cppyy
+ cppyy.load_library('FaserDetDescrDict')
+ FaserRegion = ROOT.FaserDetDescr.FaserRegion
+
+ result = ComponentAccumulator()
+ kwargs.setdefault("ParentMinEkin", 10000.*MeV)
+ kwargs.setdefault("ChildMinEkin" , 10000.*MeV)
+ kwargs.setdefault("VertexTypes", [3, 14, 15, 4, 5, 6, 7, 2, 12, 13, 111, 121, 131, 132, 141, 151, 152, 161, 310]) # EM *and* nuclear
+ kwargs.setdefault("VertexTypeRangeLow" , 201) # All kinds of decay processes
+ kwargs.setdefault("VertexTypeRangeHigh" , 298) # ...
+ kwargs.setdefault("Regions", [FaserRegion.fFaserTrench])
+ result.setPrivateTools(CompFactory.ISF.FaserTruthStrategy(name, **kwargs))
+ return result
+
# def FaserDipoleTruthStrategyCfg(ConfigFlags, name="ISF_FaserDipoleTruthStrategy", **kwargs):
# result = ComponentAccumulator()
diff --git a/Tracker/TrackerEventCnv/TrackerEventTPCnv/src/FaserSCT_SpacePointCnv_p0.cxx b/Tracker/TrackerEventCnv/TrackerEventTPCnv/src/FaserSCT_SpacePointCnv_p0.cxx
index b6cbc9da2508bd734634af3702165b262885ceff..6dc7c62a91fb18c44e701175173d3df82b4c02b0 100644
--- a/Tracker/TrackerEventCnv/TrackerEventTPCnv/src/FaserSCT_SpacePointCnv_p0.cxx
+++ b/Tracker/TrackerEventCnv/TrackerEventTPCnv/src/FaserSCT_SpacePointCnv_p0.cxx
@@ -15,7 +15,7 @@ StatusCode FaserSCT_SpacePointCnv_p0::initialize(MsgStream& log ) {
// ISvcLocator* svcLocator = Gaudi::svcLocator();
// Get the messaging service, print where you are
- log << MSG::INFO << "FaserSCT_SpacePointCnv::initialize()" << endmsg;
+ log << MSG::DEBUG << "FaserSCT_SpacePointCnv::initialize()" << endmsg;
if(m_sctClusContName.initialize()!=StatusCode::SUCCESS)
log << MSG::WARNING<< "FaserSCT_SpacePointCnv failed to initialize the sct cluster container" << endmsg;
diff --git a/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt b/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..85053a1f49bd551ea136c5cd61f6cee40543ed26
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
@@ -0,0 +1,12 @@
+atlas_subdir(PairVertex)
+
+atlas_add_component(
+ PairVertex
+ src/PairVertexAlg.h
+ src/PairVertexAlg.cxx
+ src/component/PairVertex_entries.cxx
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives
+)
+
+atlas_install_python_modules(python/*.py)
+# atlas_install_scripts(test/*.py)
diff --git a/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py b/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..f98af1629a39bfc590ef3ab0929da4b849ad5ce6
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py
@@ -0,0 +1,90 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+
+def PairVertexAlgCfg(flags, **kwargs):
+ acc = FaserSCT_GeometryCfg(flags)
+ acc.merge(MagneticFieldSvcCfg(flags))
+ # acc.merge(FaserActsTrackingGeometrySvcCfg(flags))
+ # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+
+ actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
+ actsExtrapolationTool.MaxSteps = 1000
+ actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+
+ PairVertexAlg = CompFactory.Tracker.PairVertexAlg("PairVertexAlg",**kwargs)
+ PairVertexAlg.ExtrapolationTool = actsExtrapolationTool
+ acc.addEventAlgo(PairVertexAlg)
+ return acc
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ ConfigFlags.Input.Files = [
+ '/eos/experiment/faser/sim/mdc/foresee/110004/rec/r0007/FaserMC-MDC_FS_Amm_316MeV_2Em6-110004-00000-s0003-r0007-xAOD.root'
+ ]
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(PairVertexAlgCfg(ConfigFlags))
+
+ # # 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
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/PairVertex/python/__init__.py b/Tracker/TrackerRecAlgs/PairVertex/python/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..99058166c7252ddbff852b0a09ec8d49a15e54fb
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
@@ -0,0 +1,163 @@
+#include "PairVertexAlg.h"
+
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+#include "GeoPrimitives/CLHEPtoEigenConverter.h"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+#include <cmath>
+
+
+namespace Tracker {
+
+PairVertexAlg::PairVertexAlg(const std::string &name, ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator) {}
+
+StatusCode PairVertexAlg::initialize()
+{
+ ATH_CHECK(m_mcEventCollectionKey.initialize());
+ ATH_CHECK(m_trackCollectionKey.initialize());
+ ATH_CHECK(m_extrapolationTool.retrieve());
+ ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+ ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT"));
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PairVertexAlg::execute(const EventContext &ctx) const
+{
+ m_numberOfEvents++;
+ const Acts::GeometryContext gctx =
+ m_extrapolationTool->trackingGeometryTool()->getNominalGeometryContext().context();
+
+ std::vector<double> z_positions {};
+
+ SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey, ctx};
+ ATH_CHECK(mcEvents.isValid());
+ if (mcEvents->size() != 1)
+ {
+ ATH_MSG_ERROR("There should be exactly one event in the McEventCollection.");
+ return StatusCode::FAILURE;
+ }
+
+ SG::ReadHandle<TrackCollection> tracks { m_trackCollectionKey, ctx };
+ ATH_CHECK(tracks.isValid());
+
+ const Trk::TrackParameters* positive {nullptr};
+ const Trk::TrackParameters* negative {nullptr};
+
+ for (auto trk : *tracks)
+ {
+ const Trk::TrackParameters* upstream {nullptr};
+ if (trk == nullptr)
+ {
+ ATH_MSG_WARNING("Null pointer from TrackContainer.");
+ m_invalidTrackContainerEntries++;
+ continue;
+ }
+ auto parameters = trk->trackParameters();
+ if (parameters == nullptr || parameters->empty())
+ {
+ m_numberOfNoParameterTracks++;
+ ATH_MSG_WARNING("Track without parameters found.");
+ return StatusCode::SUCCESS;
+ }
+ for (auto state : *parameters)
+ {
+ if (state == nullptr)
+ {
+ m_numberOfNullParameters++;
+ ATH_MSG_WARNING("Null track parameters returned.");
+ continue;
+ }
+ if (!state->hasSurface())
+ {
+ m_numberOfParametersWithoutSurface++;
+ ATH_MSG_WARNING("Track state without surface found.");
+ continue;
+ }
+ if (!upstream || upstream->associatedSurface().center().z() > state->associatedSurface().center().z())
+ upstream = state;
+ }
+ if (!upstream)
+ {
+ m_numberOfNoUpstreamTracks++;
+ ATH_MSG_WARNING("Unable to find track parameters on any surface");
+ continue;
+ }
+ auto momentum = upstream->momentum();
+ double charge = upstream->charge();
+ if (charge > 0 || momentum.mag() > positive->momentum().mag())
+ {
+ positive = upstream;
+ }
+ else if (charge < 0 || momentum.mag() > negative->momentum().mag())
+ {
+ negative = upstream;
+ }
+ }
+
+ if (positive != nullptr) m_numberOfPositiveTracks++;
+ if (negative != nullptr) m_numberOfNegativeTracks++;
+
+ if (positive != nullptr && negative != nullptr) m_numberOfOppositeSignPairs++;
+
+ // for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range())
+ // {
+ // if ((std::abs(particle->pdg_id()) != 13)) continue;
+ // const HepMC::FourVector& vertex = particle->production_vertex()->position();
+ // if (vertex.z() > 0) continue;
+ // const HepMC::FourVector& momentum = particle->momentum();
+ // double phi = momentum.phi();
+ // double theta = momentum.theta();
+ // double charge = particle->pdg_id() > 0 ? -1 : 1;
+ // double abs_momentum = momentum.rho() * m_MeV2GeV;
+ // double qop = charge / abs_momentum;
+ // // The coordinate system of the Acts::PlaneSurface is defined as
+ // // T = Z = normal, U = X x T = -Y, V = T x U = x
+ // // Acts::BoundVector pars;
+ // // pars << -vertex.y(), vertex.x(), phi, theta, qop, vertex.t();
+ // Acts::BoundVector pars = Acts::BoundVector::Zero();
+ // pars[Acts::eBoundLoc0] = -vertex.y();
+ // pars[Acts::eBoundLoc1] = vertex.x();
+ // pars[Acts::eBoundPhi] = phi;
+ // pars[Acts::eBoundTheta] = theta;
+ // pars[Acts::eBoundQOverP] = qop;
+ // pars[Acts::eBoundTime] = vertex.t();
+
+ // auto startSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ // Acts::Vector3(0, 0, vertex.z()), Acts::Vector3(0, 0, 1));
+ // auto targetSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ // Acts::Vector3(0, 0, z_mean), Acts::Vector3(0, 0, 1));
+ // Acts::BoundTrackParameters startParameters(
+ // std::move(startSurface), pars, charge);
+ // ATH_MSG_DEBUG("vertex: " << vertex.x() << ", " << vertex.y() << ", " << vertex.z());
+ // ATH_MSG_DEBUG("vertex momentum: " << momentum.x() * m_MeV2GeV << ", " << momentum.y() * m_MeV2GeV << ", " << momentum.z() * m_MeV2GeV);
+ // std::unique_ptr<const Acts::BoundTrackParameters> targetParameters =
+ // m_extrapolationTool->propagate(ctx, startParameters, *targetSurface);
+ // if (targetParameters)
+ // {
+ // Acts::Vector3 targetPosition = targetParameters->position(gctx);
+ // Acts::Vector3 targetMomentum = targetParameters->momentum();
+ // ATH_MSG_DEBUG("vertex: " << vertex.x() << ", " << vertex.y() << ", " << vertex.z());
+ // ATH_MSG_DEBUG("origin: " << targetPosition.x() << ", " << targetPosition.y() << ", " << targetPosition.z());
+ // ATH_MSG_DEBUG("vertex momentum: " << momentum.x() * m_MeV2GeV << ", " << momentum.y() * m_MeV2GeV << ", " << momentum.z() * m_MeV2GeV);
+ // ATH_MSG_DEBUG("origin momentum: " << targetMomentum.x() << ", " << targetMomentum.y() << ", " << targetMomentum.z());
+ // }
+ // }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PairVertexAlg::finalize()
+{
+ ATH_MSG_INFO("Found " << m_numberOfOppositeSignPairs << " opposite sign pairs in " << m_numberOfEvents << " total events.");
+ ATH_MSG_INFO(m_numberOfPositiveTracks << " events had a positive track, and " << m_numberOfNegativeTracks << " had a negative track.");
+ ATH_MSG_INFO(m_numberOfNoUpstreamTracks << " tracks could not locate an upstream position.");
+ ATH_MSG_INFO(m_numberOfNoParameterTracks << " tracks in track pairs had no track parameters.");
+ ATH_MSG_INFO(m_numberOfNullParameters << " track parameters were null and " << m_numberOfParametersWithoutSurface << " had no associated surface.");
+ ATH_MSG_INFO(m_invalidTrackContainerEntries << " invalid TrackContainer entries found.");
+ return StatusCode::SUCCESS;
+}
+
+} // Tracker
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..d15fddf6e288828d143218fa6d2f686a50076ba0
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
@@ -0,0 +1,51 @@
+/*
+Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+*/
+
+#pragma once
+
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
+#include "GeneratorObjects/McEventCollection.h"
+#include "TrkTrack/TrackCollection.h"
+
+class FaserSCT_ID;
+namespace TrackerDD
+{
+ class SCT_DetectorManager;
+}
+
+namespace Tracker
+{
+ class PairVertexAlg : public AthReentrantAlgorithm
+ {
+ public:
+ PairVertexAlg(const std::string& name, ISvcLocator* pSvcLocator);
+
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext& ctx) const override;
+ virtual StatusCode finalize() override;
+
+ private:
+ double m_MeV2GeV = 1e-3;
+ const FaserSCT_ID* m_idHelper {nullptr};
+ const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
+ SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey { this, "McEventCollection", "TruthEvent" };
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollectionKey { this, "TrackCollection", "CKFTrackCollection" };
+
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
+
+ mutable std::atomic<size_t> m_numberOfEvents{0};
+ mutable std::atomic<size_t> m_numberOfPositiveTracks{0};
+ mutable std::atomic<size_t> m_numberOfNegativeTracks{0};
+ mutable std::atomic<size_t> m_numberOfOppositeSignPairs{0};
+ mutable std::atomic<size_t> m_numberOfNoParameterTracks{0};
+ mutable std::atomic<size_t> m_invalidTrackContainerEntries{0};
+ mutable std::atomic<size_t> m_numberOfNullParameters{0};
+ mutable std::atomic<size_t> m_numberOfParametersWithoutSurface{0};
+ mutable std::atomic<size_t> m_numberOfNoUpstreamTracks{0};
+ };
+
+}
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/component/PairVertex_entries.cxx b/Tracker/TrackerRecAlgs/PairVertex/src/component/PairVertex_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..db08d4e764775c0a28da9e53538e2e5daae46932
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/component/PairVertex_entries.cxx
@@ -0,0 +1,3 @@
+#include "../PairVertexAlg.h"
+
+DECLARE_COMPONENT(Tracker::PairVertexAlg)
\ No newline at end of file
diff --git a/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsLayerBuilder.h b/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsLayerBuilder.h
index 132fd1ca2069375eceb43531071b3ff5ebd2dbdf..52aa75240b5f548433a165b19ec837827c42107a 100644
--- a/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsLayerBuilder.h
+++ b/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsLayerBuilder.h
@@ -51,9 +51,10 @@ public:
std::string configurationName = "undefined";
FaserActsDetectorElement::Subdetector subdetector
= FaserActsDetectorElement::Subdetector::SCT;
- const TrackerDD::SCT_DetectorManager* mng;
- const GeoVDetectorManager* vetomng;
- const GeoVDetectorManager* triggermng;
+ const TrackerDD::SCT_DetectorManager* mng=nullptr;
+ const GeoVDetectorManager* vetomng=nullptr;
+ const GeoVDetectorManager* triggermng=nullptr;
+ const GeoVDetectorManager* dipolemng=nullptr;
std::shared_ptr<const Acts::LayerCreator> layerCreator = nullptr;
std::shared_ptr<std::vector<std::shared_ptr<const FaserActsDetectorElement>>> elementStore;
std::shared_ptr<std::map<Identifier, Acts::GeometryIdentifier>> identifierMap;
@@ -109,15 +110,16 @@ public:
private:
- FaserActs::CuboidVolumeBuilder::VolumeConfig buildScintVolume(double , double , std::string );
+ FaserActs::CuboidVolumeBuilder::VolumeConfig buildScintVolume(const Acts::GeometryContext& ,double, double, double ,double, double, double , std::string );
+FaserActs::CuboidVolumeBuilder::VolumeConfig buildDipoleVolume(const Acts::GeometryContext& ,double , double , double , double , std::string );
double m_ModuleWidth;
double m_ModuleLength;
/// configruation object
Config m_cfg;
- Acts::Vector3 m_worldDimensions = { 400.0_mm, 400.0_mm, 8000.0_mm };
+ Acts::Vector3 m_worldDimensions = { 600.0_mm, 600.0_mm, 8000.0_mm };
Acts::Vector3 m_worldCenter = {0.0, 0.0, 0.0};
- Acts::Vector3 m_trackerDimensions = { 400.0_mm, 400.0_mm, 100.0_mm };
+ Acts::Vector3 m_trackerDimensions = { 600.0_mm, 600.0_mm, 100.0_mm };
/// Private access to the logger
const Acts::Logger&
diff --git a/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py b/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py
index 99e81e8cd5546cc7db45f34f27978989e0ec0fcd..92c1b6039a754f68dd5c562575a0dc3b9cc2e9d5 100644
--- a/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py
+++ b/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py
@@ -26,10 +26,6 @@ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
def ActsMaterialMappingCfg(configFlags, name = "FaserActsMaterialMapping", **kwargs):
result = ComponentAccumulator()
- MaterialStepConverterTool = ActsMaterialStepConverterToolCfg()
- kwargs["MaterialStepConverterTool"] = MaterialStepConverterTool.getPrimary()
- result.merge(MaterialStepConverterTool)
-
ActsSurfaceMappingTool = ActsSurfaceMappingToolCfg(configFlags)
kwargs["SurfaceMappingTool"] = ActsSurfaceMappingTool.getPrimary()
result.merge(ActsSurfaceMappingTool)
@@ -64,11 +60,10 @@ if "__main__" == __name__:
## Just enable ID for the moment.
ConfigFlags.Input.isMC = True
ConfigFlags.Beam.Type = "collisions"
- ConfigFlags.GeoModel.FaserVersion = "FASER-01"
- ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" # Always needed; must match FaserVersion
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion
ConfigFlags.TrackingGeometry.MaterialSource = "geometry-maps.json"
- ConfigFlags.Concurrency.NumThreads = 1
- ConfigFlags.Concurrency.NumConcurrentEvents = 1
+ ConfigFlags.GeoModel.Align.Dynamic = False
ConfigFlags.lock()
ConfigFlags.dump()
@@ -101,4 +96,4 @@ if "__main__" == __name__:
log.info("CONFIG DONE")
- cfg.run(80000)
+ cfg.run(-1)
diff --git a/Tracking/Acts/FaserActsGeometry/python/FaserActsWriteTrackingGeometryConfig.py b/Tracking/Acts/FaserActsGeometry/python/FaserActsWriteTrackingGeometryConfig.py
index c760be8845ae6a9afe61c262c0fb064fdb049aca..77e6921bfa63b17816a019964ec94ff80875db5c 100644
--- a/Tracking/Acts/FaserActsGeometry/python/FaserActsWriteTrackingGeometryConfig.py
+++ b/Tracking/Acts/FaserActsGeometry/python/FaserActsWriteTrackingGeometryConfig.py
@@ -31,7 +31,7 @@ def FaserActsWriteTrackingGeometryBasicCfg(flags, **kwargs):
actsWriteTrackingGeometry.ObjWriterTool=FaserActsObjWriterTool("FaserActsObjWriterTool",OutputDirectory="./", SubDetectors=["Station_0","Station_1","Station_2","Station_3"])
actsWriteTrackingGeometry.MaterialJsonWriterTool= FaserActsMaterialJsonWriterTool(OutputFile = "geometry-maps.json",
processSensitives = False,
- processnonmaterial = True)
+ processnonmaterial = False)
acc.addEventAlgo(actsWriteTrackingGeometry)
return acc
diff --git a/Tracking/Acts/FaserActsGeometry/src/CuboidVolumeBuilder.cxx b/Tracking/Acts/FaserActsGeometry/src/CuboidVolumeBuilder.cxx
index df96a311c852ff26268bb992d211679851a7fba6..e382d4cd31a023ef9a9fb4bb8afa744920bb500e 100644
--- a/Tracking/Acts/FaserActsGeometry/src/CuboidVolumeBuilder.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/CuboidVolumeBuilder.cxx
@@ -138,8 +138,6 @@ std::shared_ptr<Acts::TrackingVolume> CuboidVolumeBuilder::buildVolume(
cfg.layerCfg.push_back(lCfg);
}
- else{
- }
// Gather the layers
Acts::LayerVector layVec;
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
index 8ed2d78936f6c6fd7c36bb38225ed0a3046bcae7..c5f6a68ab01b8da95f3d62149d23bdf41684263d 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
@@ -12,11 +12,13 @@
#include "FaserActsGeometry/CuboidVolumeBuilder.h"
#include "ActsInterop/IdentityHelper.h"
#include "GeoModelKernel/GeoBox.h"
+#include "GeoModelKernel/GeoTube.h"
// ACTS
#include "Acts/Material/ProtoSurfaceMaterial.hpp"
#include "Acts/Surfaces/CylinderSurface.hpp"
#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/DiscSurface.hpp"
#include "Acts/Surfaces/RectangleBounds.hpp"
#include "Acts/Geometry/GenericApproachDescriptor.hpp"
#include "Acts/Geometry/ApproachDescriptor.hpp"
@@ -89,41 +91,63 @@ FaserActs::CuboidVolumeBuilder::Config FaserActsLayerBuilder::buildVolume(const
//veto station
//Veto A, VetoRadiator, Veto B
auto vetoManager = static_cast<const GeoVDetectorManager*>(m_cfg.vetomng);
+ if(vetoManager!=nullptr){
for(unsigned int i=0; i<vetoManager->getNumTreeTops(); i++){
auto vol0 = vetoManager->getTreeTop(i)->getLogVol();
//get the shape params and translation
const GeoBox* shape = dynamic_cast<const GeoBox*>(vol0->getShape());
auto trans = vetoManager->getTreeTop(i)->getX();
if(vetoManager->getTreeTop(i)->getNChildVols()==0){
- volumeConfigs.emplace_back(buildScintVolume(trans.translation().z(), shape->getZHalfLength(),vol0->getName()));
+ volumeConfigs.emplace_back(buildScintVolume(gctx, trans.translation().x(),trans.translation().y(),trans.translation().z(), shape->getXHalfLength(),shape->getYHalfLength(),shape->getZHalfLength(),vol0->getName()));
}
else{
for(size_t j =0 ;j< vetoManager->getTreeTop(i)->getNChildVols();j++){
auto childtrans=vetoManager->getTreeTop(i)->getXToChildVol(j);
const GeoBox* childshape = dynamic_cast<const GeoBox*>(vetoManager->getTreeTop(i)->getChildVol(j)->getLogVol()->getShape());
- volumeConfigs.emplace_back(buildScintVolume(trans.translation().z() + childtrans.translation().z(), childshape->getZHalfLength(),vol0->getName()+std::string("_")+std::to_string(j)));
+ volumeConfigs.emplace_back(buildScintVolume(gctx, trans.translation().x() + childtrans.translation().x(),trans.translation().y() + childtrans.translation().y(),trans.translation().z() + childtrans.translation().z(), childshape->getXHalfLength(), childshape->getXHalfLength(), childshape->getZHalfLength(),vol0->getName()+std::string("_")+std::to_string(j)));
}
}
}
+ }
//build trigger stations
auto triggerManager = static_cast<const GeoVDetectorManager*>(m_cfg.triggermng);
+ if(triggerManager!=nullptr){
for(unsigned int i=0; i<triggerManager->getNumTreeTops(); i++){
auto vol0 = triggerManager->getTreeTop(i)->getLogVol();
//get the shape params and translation
const GeoBox* shape = dynamic_cast<const GeoBox*>(vol0->getShape());
auto trans = triggerManager->getTreeTop(i)->getX();
if(triggerManager->getTreeTop(i)->getNChildVols()==0){
- volumeConfigs.emplace_back(buildScintVolume(trans.translation().z(), shape->getZHalfLength(),vol0->getName()));
+ volumeConfigs.emplace_back(buildScintVolume(gctx, trans.translation().x(), trans.translation().y(),trans.translation().z(), shape->getXHalfLength(), shape->getYHalfLength(), shape->getZHalfLength(),vol0->getName()));
}
else{
for(size_t j =0 ;j< triggerManager->getTreeTop(i)->getNChildVols();j++){
auto childtrans=triggerManager->getTreeTop(i)->getXToChildVol(j);
const GeoBox* childshape = dynamic_cast<const GeoBox*>(triggerManager->getTreeTop(i)->getChildVol(j)->getLogVol()->getShape());
- volumeConfigs.emplace_back(buildScintVolume(trans.translation().z() + childtrans.translation().z(), childshape->getZHalfLength(),vol0->getName()+std::string("_")+std::to_string(j)));
+ volumeConfigs.emplace_back(buildScintVolume(gctx, trans.translation().x() + childtrans.translation().x(), trans.translation().y() + childtrans.translation().y(), trans.translation().z() + childtrans.translation().z(), childshape->getXHalfLength(), childshape->getYHalfLength(), childshape->getZHalfLength(),vol0->getName()+std::string("_")+std::to_string(j)));
}
}
}
+ }
+ //build dipole magnets
+ auto dipoleManager = static_cast<const GeoVDetectorManager*>(m_cfg.dipolemng);
+ for(unsigned int i=0; i<dipoleManager->getNumTreeTops(); i++){
+ auto vol0 = dipoleManager->getTreeTop(i)->getLogVol();
+ //get the shape params and translation
+ const GeoTube* shape = dynamic_cast<const GeoTube*>(vol0->getShape());
+ auto trans = dipoleManager->getTreeTop(i)->getX();
+ if(dipoleManager->getTreeTop(i)->getNChildVols()==0){
+ volumeConfigs.emplace_back(buildDipoleVolume(gctx, trans.translation().z(), shape->getZHalfLength(), shape->getRMin(), shape->getRMax(), vol0->getName()+std::string("_")+std::to_string(i)));
+ }
+ else{
+ for(size_t j =0 ;j< dipoleManager->getTreeTop(i)->getNChildVols();j++){
+ auto childtrans=dipoleManager->getTreeTop(i)->getXToChildVol(j);
+ const GeoTube* childshape = dynamic_cast<const GeoTube*>(dipoleManager->getTreeTop(i)->getChildVol(j)->getLogVol()->getShape());
+ volumeConfigs.emplace_back(buildDipoleVolume(gctx, trans.translation().z() + childtrans.translation().z(), childshape->getZHalfLength(), childshape->getRMin(), childshape->getRMax(), vol0->getName()+std::string("_")+std::to_string(j)));
+ }
+ }
+ }
//sort the geometry by the Z position, neccessary to have the correct boundary
auto sortFunc = [](FaserActs::CuboidVolumeBuilder::VolumeConfig& v1, FaserActs::CuboidVolumeBuilder::VolumeConfig& v2){return v1.position.z()<v2.position.z();};
@@ -242,7 +266,8 @@ FaserActsLayerBuilder::buildLayers(const Acts::GeometryContext& gctx,
}
-FaserActs::CuboidVolumeBuilder::VolumeConfig FaserActsLayerBuilder::buildScintVolume(double zpos, double zhalflength, std::string name){
+
+FaserActs::CuboidVolumeBuilder::VolumeConfig FaserActsLayerBuilder::buildDipoleVolume(const Acts::GeometryContext& gctx, double zpos, double zhalflength, double rmin, double rmax, std::string name){
FaserActs::CuboidVolumeBuilder::VolumeConfig volumeConfig;
@@ -253,7 +278,6 @@ FaserActs::CuboidVolumeBuilder::VolumeConfig FaserActsLayerBuilder::buildScintVo
FaserActs::CuboidVolumeBuilder::SurfaceConfig surfacecfg;
surfacecfg.position = Acts::Vector3(0, 0, zpos);
surfacecfg.thickness=zhalflength;
- layercfg.active = true;
std::vector<std::shared_ptr<const Acts::Surface>> surfaces;
surfaces.clear();
@@ -264,37 +288,52 @@ FaserActs::CuboidVolumeBuilder::VolumeConfig FaserActsLayerBuilder::buildScintVo
layercfg.active = true;
//bounds are hard coded
- auto rBounds = std::make_shared<const Acts::RectangleBounds>( 130.*1_mm, 130.*1_mm) ;
- surfacecfg.rBounds=rBounds;
layercfg.surfaceCfg = surfacecfg;
+ Transform3 transformCenter(Translation3(0., 0., (zpos)*1_mm));
- Transform3 transformInner(Translation3(0., 0., (zpos)*1_mm));
+ std::shared_ptr<Acts::DiscSurface> centerBoundary =
+ Acts::Surface::makeShared<Acts::DiscSurface>(
+ transformCenter, rmin, rmax);
- std::shared_ptr<Acts::PlaneSurface> innerBoundary
- = Acts::Surface::makeShared<Acts::PlaneSurface>(transformInner, rBounds);
+ Transform3 transformInner(Translation3(0., 0., (zpos-zhalflength)*1_mm));
+
+ std::shared_ptr<Acts::DiscSurface> innerBoundary =
+ Acts::Surface::makeShared<Acts::DiscSurface>(
+ transformInner, rmin, rmax);
+
+ Transform3 transformOuter(Translation3(0., 0., (zpos+zhalflength)*1_mm));
+ std::shared_ptr<Acts::DiscSurface> outerBoundary =
+ Acts::Surface::makeShared<Acts::DiscSurface>(
+ transformOuter, rmin, rmax);
- size_t matBinsX = m_cfg.MaterialBins.first;
+
+// size_t matBinsX = m_cfg.MaterialBins.first;
size_t matBinsY = m_cfg.MaterialBins.second;
+// set bin size to 1 at phi and MaterialBins-Y at r direction
+ Acts::BinUtility materialBinUtil(1, -M_PI, M_PI, Acts::closed, Acts::binPhi);
+ materialBinUtil +=
+ Acts::BinUtility(matBinsY, rmin, rmax, Acts::open, Acts::binR, transformInner);
+
+ std::shared_ptr<const Acts::ProtoSurfaceMaterial> materialProxy =
+ std::make_shared<const Acts::ProtoSurfaceMaterial>(materialBinUtil);
- Acts::BinUtility materialBinUtil(
- matBinsX, -130*1_mm, 130.*1_mm, Acts::open, Acts::binX);
- materialBinUtil += Acts::BinUtility(
- matBinsY, -130*1_mm, 130.*1_mm, Acts::open, Acts::binY, transformInner);
- std::shared_ptr<const Acts::ProtoSurfaceMaterial> materialProxy
- = std::make_shared<const Acts::ProtoSurfaceMaterial>(materialBinUtil);
innerBoundary->assignSurfaceMaterial(materialProxy);
std::vector<std::shared_ptr<const Acts::Surface>> aSurfaces;
aSurfaces.push_back(std::move(innerBoundary));
+// aSurfaces.push_back(std::move(centerBoundary));
+// aSurfaces.push_back(std::move(outerBoundary));
+ Acts::ProtoLayer pl(gctx, aSurfaces);
layercfg.surfaces = aSurfaces;
- layerConfigs.push_back(layercfg);
+ layercfg.approachDescriptor = new Acts::GenericApproachDescriptor(aSurfaces);
+ layerConfigs.push_back(layercfg);
volumeConfig.position = Acts::Vector3(0, 0, zpos);
- Acts::Vector3 length= { 300.0*1_mm, 300.0*1_mm, zhalflength*2.*1_mm };
+ Acts::Vector3 length= { 500.0*1_mm, 500.0*1_mm, zhalflength*2.*1_mm };
volumeConfig.length = length;
//volumeConfig.layerCfg = {};
volumeConfig.layerCfg = layerConfigs;
@@ -304,4 +343,79 @@ FaserActs::CuboidVolumeBuilder::VolumeConfig FaserActsLayerBuilder::buildScintVo
}
+FaserActs::CuboidVolumeBuilder::VolumeConfig FaserActsLayerBuilder::buildScintVolume(const Acts::GeometryContext& gctx, double xpos, double ypos, double zpos, double xhalflength, double yhalflength, double zhalflength, std::string name){
+
+ FaserActs::CuboidVolumeBuilder::VolumeConfig volumeConfig;
+
+ std::vector<FaserActs::CuboidVolumeBuilder::LayerConfig> layerConfigs;
+ layerConfigs.clear();
+
+ FaserActs::CuboidVolumeBuilder::LayerConfig layercfg;
+ FaserActs::CuboidVolumeBuilder::SurfaceConfig surfacecfg;
+ surfacecfg.position = Acts::Vector3(xpos, ypos, zpos);
+ surfacecfg.thickness=zhalflength;
+
+ std::vector<std::shared_ptr<const Acts::Surface>> surfaces;
+ surfaces.clear();
+
+ layercfg.binsX = 1;
+ layercfg.binsY = 1;
+
+ layercfg.active = true;
+
+ //bounds are hard coded
+ layercfg.surfaceCfg = surfacecfg;
+
+ auto rBounds = std::make_shared<const Acts::RectangleBounds>( xhalflength*1_mm, yhalflength*1_mm) ;
+// surfacecfg.rBounds=rBounds;
+
+
+ Transform3 transformInner(Translation3(xpos*1_mm, ypos*1_mm, (zpos-zhalflength)*1_mm));
+
+ std::shared_ptr<Acts::PlaneSurface> innerBoundary =
+ Acts::Surface::makeShared<Acts::PlaneSurface>(
+ transformInner, rBounds);
+
+ Transform3 transformCenter(Translation3(xpos*1_mm, ypos*1_mm, zpos*1_mm));
+ std::shared_ptr<Acts::PlaneSurface> centerBoundary =
+ Acts::Surface::makeShared<Acts::PlaneSurface>(
+ transformCenter, rBounds);
+
+ Transform3 transformOuter(Translation3(xpos*1_mm, ypos*1_mm, (zpos+zhalflength)*1_mm));
+
+ std::shared_ptr<Acts::PlaneSurface> outerBoundary =
+ Acts::Surface::makeShared<Acts::PlaneSurface>(
+ transformOuter, rBounds);
+
+
+// set bin size to 2
+ Acts::BinUtility materialBinUtil(2, (0.-xhalflength)*1_mm, xhalflength*1_mm, Acts::open, Acts::binX);
+ materialBinUtil +=
+ Acts::BinUtility(2, (0.-yhalflength)*1_mm, yhalflength*1_mm, Acts::open, Acts::binY, transformInner);
+
+ std::shared_ptr<const Acts::ProtoSurfaceMaterial> materialProxy =
+ std::make_shared<const Acts::ProtoSurfaceMaterial>(materialBinUtil);
+
+
+ innerBoundary->assignSurfaceMaterial(materialProxy);
+
+ std::vector<std::shared_ptr<const Acts::Surface>> aSurfaces;
+ aSurfaces.push_back(std::move(innerBoundary));
+ aSurfaces.push_back(std::move(centerBoundary));
+ aSurfaces.push_back(std::move(outerBoundary));
+ Acts::ProtoLayer pl(gctx, aSurfaces);
+
+ layercfg.surfaces = aSurfaces;
+ layercfg.approachDescriptor = new Acts::GenericApproachDescriptor(aSurfaces);
+
+ layerConfigs.push_back(layercfg);
+ volumeConfig.position = Acts::Vector3(0, 0, zpos);
+ Acts::Vector3 length= { 500.0*1_mm, 500.0*1_mm, zhalflength*2.*1_mm };
+ volumeConfig.length = length;
+ //volumeConfig.layerCfg = {};
+ volumeConfig.layerCfg = layerConfigs;
+ volumeConfig.name = name;
+
+ return volumeConfig;
+}
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.cxx
index 9edaa72ff056bbeab111cce25c127f9f7382e9a8..2595b42dc613db98997d4ca5c8ded263e65a9963 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.cxx
@@ -79,7 +79,8 @@ StatusCode FaserActsMaterialMapping::execute(const EventContext &ctx) const {
const Trk::MaterialStep* mat = *iter;
Trk::MaterialStep* newmat = const_cast<Trk::MaterialStep*>(mat);
//remove the magnets, need to be updated
- if((newmat->hitZ()<-1600&&newmat->hitZ()>-1920)||(newmat->hitZ()>-200&&newmat->hitZ()<100)||(newmat->hitZ()>1000&&newmat->hitZ()<1300)||(newmat->hitZ()>2000&&newmat->hitZ()<2500))
+// if((newmat->hitZ()<-1600&&newmat->hitZ()>-1920)||(newmat->hitZ()>-200&&newmat->hitZ()<100)||(newmat->hitZ()>1000&&newmat->hitZ()<1300)||(newmat->hitZ()>2000&&newmat->hitZ()<2500))
+ if(newmat->hitZ()>-1920&&newmat->hitZ()<2500)
newmaterialStepCollection->push_back(newmat);
}
if(newmaterialStepCollection->size()<1)
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx
index c6f476e7c1292f81e11916f24aa25d1b7a993ded..a45c7c06d37214248a19b189ca40886339c9f1b4 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx
@@ -63,9 +63,10 @@ FaserActsTrackingGeometrySvc::initialize()
//get all the subdetectors
const GeoModelExperiment* theExpt = nullptr;
- ATH_CHECK( m_detStore->retrieve( theExpt ));
+ ATH_CHECK( m_detStore->retrieve( theExpt ,"FASER"));
const GeoVDetectorManager *vetoManager = theExpt->getManager("Veto");
const GeoVDetectorManager *triggerManager = theExpt->getManager("Trigger");
+ const GeoVDetectorManager *dipoleManager = theExpt->getManager("Dipole");
Acts::LayerArrayCreator::Config lacCfg;
auto layerArrayCreator = std::make_shared<const Acts::LayerArrayCreator>(
@@ -104,7 +105,7 @@ FaserActsTrackingGeometrySvc::initialize()
// SCT
tgbConfig.trackingVolumeBuilders.push_back([&](
const auto& gctx, const auto& /*inner*/, const auto&) {
- auto tv = makeVolumeBuilder(gctx, p_SCTManager, vetoManager, triggerManager);
+ auto tv = makeVolumeBuilder(gctx, p_SCTManager, vetoManager, triggerManager, dipoleManager);
return tv->trackingVolume(gctx);
});
}
@@ -151,7 +152,7 @@ FaserActsTrackingGeometrySvc::trackingGeometry() {
}
std::shared_ptr<const Acts::ITrackingVolumeBuilder>
-FaserActsTrackingGeometrySvc::makeVolumeBuilder(const Acts::GeometryContext& gctx, const TrackerDD::SCT_DetectorManager* manager, const GeoVDetectorManager* vetoManager, const GeoVDetectorManager* triggerManager)
+FaserActsTrackingGeometrySvc::makeVolumeBuilder(const Acts::GeometryContext& gctx, const TrackerDD::SCT_DetectorManager* manager, const GeoVDetectorManager* vetoManager, const GeoVDetectorManager* triggerManager, const GeoVDetectorManager* dipoleManager)
{
std::string managerName = manager->getName();
@@ -161,8 +162,18 @@ FaserActsTrackingGeometrySvc::makeVolumeBuilder(const Acts::GeometryContext& gct
cfg.subdetector = FaserActsDetectorElement::Subdetector::SCT;
cfg.mng = static_cast<const TrackerDD::SCT_DetectorManager*>(manager);
//get veto and trigger manager
- cfg.vetomng = static_cast<const GeoVDetectorManager*>(vetoManager);
- cfg.triggermng = static_cast<const GeoVDetectorManager*>(triggerManager);
+ if(vetoManager!=nullptr)
+ cfg.vetomng = static_cast<const GeoVDetectorManager*>(vetoManager);
+ else
+ ATH_MSG_WARNING(name() << " can not find veto stations");
+ if(triggerManager!=nullptr)
+ cfg.triggermng = static_cast<const GeoVDetectorManager*>(triggerManager);
+ else
+ ATH_MSG_WARNING(name() << " can not find trigger stations");
+ if(dipoleManager!=nullptr)
+ cfg.dipolemng = static_cast<const GeoVDetectorManager*>(dipoleManager);
+ else
+ ATH_MSG_WARNING(name() << " can not find dipoles");
cfg.elementStore = m_elementStore;
cfg.identifierMap = m_identifierMap;
std::vector<size_t> matBins(m_MaterialBins);
@@ -181,8 +192,6 @@ FaserActsTrackingGeometrySvc::populateAlignmentStore(FaserActsAlignmentStore *st
// populate the alignment store with all detector elements
m_trackingGeometry->visitSurfaces(
[store](const Acts::Surface* srf) {
- FaserActsGeometryContext constructionContext;
- constructionContext.construction = true;
const Acts::DetectorElementBase* detElem = srf->associatedDetectorElement();
if(detElem){
const auto* gmde = dynamic_cast<const FaserActsDetectorElement*>(detElem);
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.h
index 83ef0efaefff6bf315f05bbc7f10bed8b8bd8665..3b191202b5959fc0ee606f7d2a21c4a82155044e 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.h
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.h
@@ -71,7 +71,7 @@ public:
private:
std::shared_ptr<const Acts::ITrackingVolumeBuilder>
- makeVolumeBuilder(const Acts::GeometryContext& gctx, const TrackerDD::SCT_DetectorManager* manager, const GeoVDetectorManager* vetoManager, const GeoVDetectorManager* triggerManager);
+ makeVolumeBuilder(const Acts::GeometryContext& gctx, const TrackerDD::SCT_DetectorManager* manager, const GeoVDetectorManager* vetoManager, const GeoVDetectorManager* triggerManager, const GeoVDetectorManager* dipoleManager);
ServiceHandle<StoreGateSvc> m_detStore;
const TrackerDD::SCT_DetectorManager* p_SCTManager;
diff --git a/Tracking/Acts/FaserActsGeometry/test/FaserActsWriteTrackingGeometry.py b/Tracking/Acts/FaserActsGeometry/test/FaserActsWriteTrackingGeometry.py
index 21dfba7e51488c1f432807e43c6bdbce5759b1e0..ed998b3b4607cf633e2825b5cb868c08039308d6 100644
--- a/Tracking/Acts/FaserActsGeometry/test/FaserActsWriteTrackingGeometry.py
+++ b/Tracking/Acts/FaserActsGeometry/test/FaserActsWriteTrackingGeometry.py
@@ -17,7 +17,7 @@ log.setLevel(DEBUG)
Configurable.configurableRun3Behavior = True
# Configure
-ConfigFlags.Input.Files = ["my03.HITS.pool.root"]
+ConfigFlags.Input.Files = ["myevt4.HITS.pool.root"]
#ConfigFlags.Output.RDOFileName = "myRDO_sp.pool.root"
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion
ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Always needed
diff --git a/faser-common b/faser-common
index 69a90ec95da88a00097fb809bede6c2bae8c02d6..89ce6a07128eb2ebc367b6b68f29c9c88220e3e6 160000
--- a/faser-common
+++ b/faser-common
@@ -1 +1 @@
-Subproject commit 69a90ec95da88a00097fb809bede6c2bae8c02d6
+Subproject commit 89ce6a07128eb2ebc367b6b68f29c9c88220e3e6