diff --git a/Generators/GeneratorUtils/python/ShiftLOS.py b/Generators/GeneratorUtils/python/ShiftLOS.py
index 4327b3f3afa7c559b97c01ac94b1f730157cb0fc..44a94c7a2572bdef501c4f8954a79c3fa8e2a544 100644
--- a/Generators/GeneratorUtils/python/ShiftLOS.py
+++ b/Generators/GeneratorUtils/python/ShiftLOS.py
@@ -15,8 +15,8 @@ class ShiftLOS(PyAthena.Alg):
super(ShiftLOS,self).__init__(name=name)
self.InputMCEventKey = InputMCEventKey
self.OutputMCEventKey = OutputMCEventKey
- self.xcross = xcross * 1e-6
- self.ycross = ycross * 1e-6
+ self.xcross = xcross
+ self.ycross = ycross
self.xshift = xshift
self.yshift = yshift
self.distance = 480 * m # Assumes 480m is 0 of FASER coordinate system
@@ -163,7 +163,7 @@ if __name__ == "__main__":
acc = ComponentAccumulator()
alg = ShiftLOS("ShiftLOS", InputMCEventKey=args.InputMCEventKey, OutputMCEventKey=args.OutputMCEventKey,
- xcross = args.xcross, ycross = args.ycross, xshift = args.xshift, yshift = args.yshift)
+ xcross = args.xcross * 1e-6, ycross = args.ycross * 1e-6, xshift = args.xshift, yshift = args.yshift)
alg.OutputLevel = INFO
acc.addEventAlgo(alg)
cfg.merge(acc)
diff --git a/Tracker/TrackerConditions/FaserSCT_ConditionsTools/CMakeLists.txt b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/CMakeLists.txt
index b19cee24e9260527698b0040339c7560a97feb49..1450b0a07133ff3ea8a5517244cbdb9a54b0262a 100644
--- a/Tracker/TrackerConditions/FaserSCT_ConditionsTools/CMakeLists.txt
+++ b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/CMakeLists.txt
@@ -19,6 +19,8 @@ atlas_add_component ( FaserSCT_ConditionsTools
atlas_add_library( FaserSCT_ConditionsToolsLib
+ FaserSCT_ConditionsTools/*.h
+ src/*.h
src/*.cxx
PUBLIC_HEADERS FaserSCT_ConditionsTools
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
diff --git a/Tracker/TrackerConditions/FaserSCT_ConditionsTools/FaserSCT_ConditionsTools/ISCT_NoisyStripTool.h b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/FaserSCT_ConditionsTools/ISCT_NoisyStripTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..9fcf6b9a54f891e52edcedb729a3ff8510af319a
--- /dev/null
+++ b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/FaserSCT_ConditionsTools/ISCT_NoisyStripTool.h
@@ -0,0 +1,23 @@
+/*
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS and FAsER collaborations
+*/
+
+#ifndef ISCT_NOISY_STRIP_TOOL
+#define ISCT_NOISY_STRIP_TOOL
+
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/EventContext.h"
+#include <vector>
+
+
+class ISCT_NoisyStripTool: virtual public IAlgTool {
+public:
+ virtual ~ISCT_NoisyStripTool() = default;
+
+ DeclareInterfaceID(ISCT_NoisyStripTool, 1, 0);
+
+ virtual std::map<std::pair<int,int>, double> getNoisyStrips(const EventContext& ctx) const = 0;
+ virtual std::map<std::pair<int,int>, double> getNoisyStrips(void) const = 0;
+};
+
+#endif // ISCT_NOISY_STRIP_TOOL
diff --git a/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/FaserSCT_NoisyStripTool.cxx b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/FaserSCT_NoisyStripTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0e2b5c444a8839e0421a073fd79b0737a89ed635
--- /dev/null
+++ b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/FaserSCT_NoisyStripTool.cxx
@@ -0,0 +1,113 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS and FASER collaborations
+*/
+
+#include "FaserSCT_NoisyStripTool.h"
+
+
+FaserSCT_NoisyStripTool::FaserSCT_NoisyStripTool (const std::string& type,
+ const std::string& name, const IInterface* parent) :
+ base_class(type, name, parent) {}
+
+
+StatusCode FaserSCT_NoisyStripTool::initialize() {
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode FaserSCT_NoisyStripTool::finalize() {
+ return StatusCode::SUCCESS;
+}
+
+
+std::map<std::pair<int, int>, double>
+FaserSCT_NoisyStripTool::getNoisyStrips(const EventContext& /*ctx*/) const {
+ // TODO fix hard-coded definition of noisy strip and read from NoisyPixels.xml database
+ std::map<std::pair<int, int>, double> noisyStrips {};
+ // noisyStrips.insert({std::make_pair(10, 150), 0.25279282895176935});
+ // noisyStrips.insert({std::make_pair(20, 155), 0.9950721341449819});
+ // noisyStrips.insert({std::make_pair(38, 41), 1.0});
+ // noisyStrips.insert({std::make_pair(48, 643), 0.6209110977322898});
+ // noisyStrips.insert({std::make_pair(49, 69), 1.0});
+ // noisyStrips.insert({std::make_pair(49, 508), 0.1250027872544429});
+ // noisyStrips.insert({std::make_pair(49, 660), 1.0});
+ // noisyStrips.insert({std::make_pair(61, 350), 0.3669587709322809});
+ // noisyStrips.insert({std::make_pair(61, 351), 0.32956496532655477});
+ // noisyStrips.insert({std::make_pair(64, 287), 1.0});
+ // noisyStrips.insert({std::make_pair(65, 323), 0.5987691484380226});
+ // noisyStrips.insert({std::make_pair(67, 147), 1.0});
+ // noisyStrips.insert({std::make_pair(67, 207), 1.0});
+ // noisyStrips.insert({std::make_pair(67, 346), 0.6463977523580173});
+ // noisyStrips.insert({std::make_pair(83, 114), 0.9968113809173412});
+ // noisyStrips.insert({std::make_pair(86, 544), 0.46609583695676415});
+ // noisyStrips.insert({std::make_pair(96, 79), 1.0});
+ // noisyStrips.insert({std::make_pair(96, 183), 1.0});
+ // noisyStrips.insert({std::make_pair(97, 550), 1.0});
+ // noisyStrips.insert({std::make_pair(100, 215), 1.0});
+ // noisyStrips.insert({std::make_pair(100, 610), 1.0});
+ // noisyStrips.insert({std::make_pair(130, 722), 1.0});
+ // noisyStrips.insert({std::make_pair(132, 297), 0.8765803732691151});
+ // noisyStrips.insert({std::make_pair(144, 597), 1.0});
+ // noisyStrips.insert({std::make_pair(144, 665), 1.0});
+ // noisyStrips.insert({std::make_pair(145, 9), 1.0});
+ // noisyStrips.insert({std::make_pair(145, 492), 1.0});
+ // noisyStrips.insert({std::make_pair(145, 566), 1.0});
+ // noisyStrips.insert({std::make_pair(146, 267), 1.0});
+ // noisyStrips.insert({std::make_pair(147, 393), 0.4758623765246282});
+ // noisyStrips.insert({std::make_pair(147, 394), 0.5172252324570206});
+ // noisyStrips.insert({std::make_pair(147, 395), 0.5058532343300556});
+ // noisyStrips.insert({std::make_pair(147, 396), 0.5272816464869445});
+ // noisyStrips.insert({std::make_pair(147, 397), 0.4782036702566504});
+ // noisyStrips.insert({std::make_pair(147, 398), 0.5092202376970589});
+ // noisyStrips.insert({std::make_pair(147, 399), 0.4811247129127924});
+ // noisyStrips.insert({std::make_pair(147, 600), 1.0});
+ // noisyStrips.insert({std::make_pair(151, 417), 0.6434767097018753});
+ // noisyStrips.insert({std::make_pair(152, 698), 0.2651013445715433});
+ // noisyStrips.insert({std::make_pair(152, 699), 0.3080250629919504});
+ // noisyStrips.insert({std::make_pair(153, 396), 0.4862532610876982});
+ // noisyStrips.insert({std::make_pair(153, 397), 0.44509108747519344});
+ // noisyStrips.insert({std::make_pair(154, 620), 0.4075634936562089});
+ // noisyStrips.insert({std::make_pair(154, 621), 0.5001449372310299});
+ // noisyStrips.insert({std::make_pair(155, 146), 0.4566637679220461});
+ // noisyStrips.insert({std::make_pair(155, 147), 0.4941021695988584});
+ // noisyStrips.insert({std::make_pair(158, 737), 0.5486208665016612});
+ // noisyStrips.insert({std::make_pair(158, 738), 0.5591901353490758});
+ // noisyStrips.insert({std::make_pair(158, 739), 0.5590786451713604});
+ // noisyStrips.insert({std::make_pair(160, 7), 1.0});
+ // noisyStrips.insert({std::make_pair(161, 763), 1.0});
+ // noisyStrips.insert({std::make_pair(164, 613), 1.0});
+ // noisyStrips.insert({std::make_pair(167, 175), 1.0});
+ // noisyStrips.insert({std::make_pair(170, 90), 0.48484848484848486});
+ // noisyStrips.insert({std::make_pair(170, 91), 0.4570874305973644});
+ // noisyStrips.insert({std::make_pair(173, 18), 1.0});
+ // noisyStrips.insert({std::make_pair(173, 484), 1.0});
+ // noisyStrips.insert({std::make_pair(174, 230), 1.0});
+ // noisyStrips.insert({std::make_pair(174, 530), 1.0});
+ // noisyStrips.insert({std::make_pair(175, 683), 1.0});
+ // noisyStrips.insert({std::make_pair(176, 418), 1.0});
+ // noisyStrips.insert({std::make_pair(177, 149), 1.0});
+ // noisyStrips.insert({std::make_pair(177, 345), 1.0});
+ // noisyStrips.insert({std::make_pair(178, 214), 1.0});
+ // noisyStrips.insert({std::make_pair(178, 508), 0.5192097576203537});
+ // noisyStrips.insert({std::make_pair(178, 673), 0.6028496889424042});
+ // noisyStrips.insert({std::make_pair(179, 651), 0.999977701964457});
+ // noisyStrips.insert({std::make_pair(182, 525), 0.5044707561263853});
+ // noisyStrips.insert({std::make_pair(182, 526), 0.5083506143108792});
+ // noisyStrips.insert({std::make_pair(185, 493), 0.42738644725399694});
+ // noisyStrips.insert({std::make_pair(185, 494), 0.43757664949717934});
+ // noisyStrips.insert({std::make_pair(187, 427), 0.9203737150757019});
+ // noisyStrips.insert({std::make_pair(188, 696), 0.6201752625593685});
+ // noisyStrips.insert({std::make_pair(188, 752), 1.0});
+ // noisyStrips.insert({std::make_pair(189, 249), 0.1250027872544429});
+ // noisyStrips.insert({std::make_pair(189, 338), 0.25925925925925924});
+ // noisyStrips.insert({std::make_pair(191, 170), 1.0});
+ // noisyStrips.insert({std::make_pair(191, 406), 1.0});
+ return noisyStrips;
+}
+
+
+std::map<std::pair<int, int>, double>
+FaserSCT_NoisyStripTool::getNoisyStrips() const {
+ const EventContext& ctx{Gaudi::Hive::currentContext()};
+ return getNoisyStrips(ctx);
+}
diff --git a/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/FaserSCT_NoisyStripTool.h b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/FaserSCT_NoisyStripTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..c554620885b0cf28beb5c3d93c6d6fa31377e3a4
--- /dev/null
+++ b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/FaserSCT_NoisyStripTool.h
@@ -0,0 +1,25 @@
+/*
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS and CERN collaborations
+*/
+
+#ifndef FASERSCT_NOISY_STRIP_TOOL
+#define FASERSCT_NOISY_STRIP_TOOL
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserSCT_ConditionsTools/ISCT_NoisyStripTool.h"
+#include "GaudiKernel/ICondSvc.h"
+#include "GaudiKernel/EventContext.h"
+
+
+class FaserSCT_NoisyStripTool: public extends<AthAlgTool, ISCT_NoisyStripTool> {
+public:
+ FaserSCT_NoisyStripTool(const std::string& type, const std::string& name, const IInterface* parent);
+ virtual ~FaserSCT_NoisyStripTool() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+
+ virtual std::map<std::pair<int, int>, double> getNoisyStrips(const EventContext& ctx) const override;
+ virtual std::map<std::pair<int, int>, double> getNoisyStrips() const override;
+};
+
+#endif // FASERSCT_NOISY_STRIP_TOOL
diff --git a/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/components/FaserSCT_ConditionsTools_entries.cxx b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/components/FaserSCT_ConditionsTools_entries.cxx
index cae3df3c65e74cecd79c4729d2dc57cc8e7e4cf3..adad7ee37337a3d347bedc9c1f5b28cf731fb837 100644
--- a/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/components/FaserSCT_ConditionsTools_entries.cxx
+++ b/Tracker/TrackerConditions/FaserSCT_ConditionsTools/src/components/FaserSCT_ConditionsTools_entries.cxx
@@ -20,6 +20,7 @@
// #include "../SCT_StripVetoTool.h"
// #include "../SCT_TdaqEnabledTool.h"
#include "../FaserSCT_CableMappingTool.h"
+#include "../FaserSCT_NoisyStripTool.h"
// DECLARE_COMPONENT( SCT_ByteStreamErrorsTool )
// DECLARE_COMPONENT( SCT_ChargeTrappingTool )
@@ -42,4 +43,5 @@ DECLARE_COMPONENT( FaserSCT_ReadCalibChipDataTool )
DECLARE_COMPONENT( FaserSCT_SiliconConditionsTool )
// DECLARE_COMPONENT( SCT_StripVetoTool )
// DECLARE_COMPONENT( SCT_TdaqEnabledTool )
-DECLARE_COMPONENT( FaserSCT_CableMappingTool )
\ No newline at end of file
+DECLARE_COMPONENT( FaserSCT_CableMappingTool )
+DECLARE_COMPONENT( FaserSCT_NoisyStripTool )
diff --git a/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerByteStreamCnv.cxx b/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerByteStreamCnv.cxx
index bd6ff38f7fc46fe44c62e5ccb339744e46cad90c..87bb2437c36c0d0c4f2cbe0cbfb4f3a1270c1e1a 100644
--- a/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerByteStreamCnv.cxx
+++ b/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerByteStreamCnv.cxx
@@ -28,6 +28,7 @@ TrackerByteStreamCnv::TrackerByteStreamCnv(ISvcLocator* svcloc)
, m_tool("TrackerDataDecoderTool")
, m_mappingTool("FaserSCT_CableMappingTool")
, m_rdpSvc("FaserROBDataProviderSvc", "TrackerByteStreamCnv")
+ , m_noisyStripTool("FaserSCT_NoisyStripTool")
, m_detStoreSvc("StoreGateSvc/DetectorStore", "TrackerByteStreamCnv")
{
}
@@ -48,6 +49,7 @@ StatusCode TrackerByteStreamCnv::initialize()
CHECK(m_rdpSvc.retrieve());
CHECK(m_tool.retrieve());
CHECK(m_mappingTool.retrieve());
+ CHECK(m_noisyStripTool.retrieve());
ATH_CHECK(m_detStoreSvc.retrieve());
ATH_CHECK(m_detStoreSvc->retrieve(m_sctID, "FaserSCT_ID"));
@@ -98,10 +100,13 @@ StatusCode TrackerByteStreamCnv::createObj(IOpaqueAddress* pAddr, DataObject*& p
auto mapping = m_mappingTool->getCableMapping();
ATH_MSG_DEBUG("Cable mapping contains " << mapping.size() << " entries");
-
+
+ auto noisyStrips = m_noisyStripTool->getNoisyStrips();
+ ATH_MSG_DEBUG(noisyStrips.size() << " noisy strips");
+
// Convert raw data into this container
- CHECK( m_tool->convert(re, cont, key, mapping) );
+ CHECK( m_tool->convert(re, cont, key, mapping, noisyStrips) );
pObj = SG::asStorable(cont);
diff --git a/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerByteStreamCnv.h b/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerByteStreamCnv.h
index de0095b637bc99a9b0e30e2b34fba5842bcdbcf6..86262fbc244530aecac352d968a13d54d0943728 100644
--- a/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerByteStreamCnv.h
+++ b/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerByteStreamCnv.h
@@ -15,6 +15,7 @@
#include "AthenaBaseComps/AthMessaging.h"
#include "FaserByteStreamCnvSvcBase/FaserByteStreamAddress.h"
#include "FaserSCT_ConditionsTools/ISCT_CableMappingTool.h"
+#include "FaserSCT_ConditionsTools/ISCT_NoisyStripTool.h"
class TrackerDataDecoderTool;
class IFaserROBDataProviderSvc;
@@ -40,10 +41,11 @@ public:
static const CLID& classID();
private:
- ToolHandle<TrackerDataDecoderTool> m_tool;
- ToolHandle<ISCT_CableMappingTool> m_mappingTool;
ServiceHandle<IFaserROBDataProviderSvc> m_rdpSvc;
ServiceHandle<StoreGateSvc> m_detStoreSvc;
+ ToolHandle<TrackerDataDecoderTool> m_tool;
+ ToolHandle<ISCT_CableMappingTool> m_mappingTool;
+ ToolHandle<ISCT_NoisyStripTool> m_noisyStripTool;
const FaserSCT_ID* m_sctID{nullptr};
};
diff --git a/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerDataDecoderTool.cxx b/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerDataDecoderTool.cxx
index 9971d0292356a1e8437f5734c2703a0bac10ee2a..612978a068a931bf4a85b267537d89b626ba53e5 100644
--- a/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerDataDecoderTool.cxx
+++ b/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerDataDecoderTool.cxx
@@ -45,7 +45,7 @@ TrackerDataDecoderTool::initialize()
ATH_CHECK(detStore()->retrieve(m_sctID, "FaserSCT_ID"));
- auto first_wafer = m_sctID->wafer_begin();
+ // auto first_wafer = m_sctID->wafer_begin();
// m_trb0Station = m_sctID->station(*first_wafer);
m_sctContext = m_sctID->wafer_context();
@@ -81,7 +81,9 @@ StatusCode
TrackerDataDecoderTool::convert(const DAQFormats::EventFull* re,
FaserSCT_RDO_Container* container,
std::string key,
- const std::map<int, std::pair<int, int> >& cableMapping)
+ const std::map<int, std::pair<int, int> >& cableMapping,
+ const std::map<std::pair<int, int>, double>& noisyStrips)
+
{
ATH_MSG_DEBUG("TrackerDataDecoderTool::convert()");
@@ -231,6 +233,12 @@ TrackerDataDecoderTool::convert(const DAQFormats::EventFull* re,
ATH_MSG_ERROR("Invalid strip number on side: " << stripOnSide);
continue;
}
+ // check if strip is noisy
+ auto it = noisyStrips.find(std::make_pair(waferHash, stripOnSide));
+ if (it != noisyStrips.end() && it->second >= m_occupancyCut) {
+ ATH_MSG_VERBOSE("Mask wafer " << waferHash << ", strip " << stripOnSide << " with an occupancy of " << it->second);
+ continue;
+ }
Identifier digitID {m_sctID->strip_id(id, stripOnSide)};
int errors{0};
int groupSize{1};
diff --git a/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerDataDecoderTool.h b/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerDataDecoderTool.h
index 32b42ced60b494d93661aa4f3ea8dbcbdc6ce63b..3bc337276a3f8eaba28fa51d90f7eb4f903592c7 100644
--- a/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerDataDecoderTool.h
+++ b/Tracker/TrackerEventCnv/TrackerByteStream/src/TrackerDataDecoderTool.h
@@ -31,7 +31,9 @@ class TrackerDataDecoderTool : public AthAlgTool {
virtual StatusCode initialize();
virtual StatusCode finalize();
- StatusCode convert(const DAQFormats::EventFull* re, FaserSCT_RDO_Container* cont, std::string key, const std::map<int, std::pair<int, int> >& cableMapping);
+ StatusCode convert(const DAQFormats::EventFull* re, FaserSCT_RDO_Container* cont, std::string key,
+ const std::map<int, std::pair<int, int> >& cableMapping,
+ const std::map<std::pair<int, int>, double>& noisyStrips);
private:
const FaserSCT_ID* m_sctID{nullptr};
@@ -41,7 +43,8 @@ private:
"ModuleMap",
{7, 2, 5, 0, 3, 6, 1, 4},
"Mapping from online to offline module numbers" };
-
+ Gaudi::Property<double> m_occupancyCut {this, "OccupancyCut", 0.1, "Mask strips with an occupancy larger than the OccupancyCut"};
+
// Gaudi::Property<uint32_t> m_trb0Station { this, "Trb0StationNumber", 1, "Station number for TRB #0" };
};
diff --git a/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsConfig.py b/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsConfig.py
index 36a2008a833453937a2adf9c4c1c2f3d81525362..7fa948ab3072fbbd11340eafb6c6c59cd92595dc 100644
--- a/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsConfig.py
+++ b/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsConfig.py
@@ -1,5 +1,5 @@
"""
-Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/Tracker/TrackerRecAlgs/MCEvents/python/MCEventsConfig.py b/Tracker/TrackerRecAlgs/MCEvents/python/MCEventsConfig.py
index 87fabb5738421fa804d7b6eced90c6a074a3829c..215d7a0ba1b878434beeb126931b32d832699155 100644
--- a/Tracker/TrackerRecAlgs/MCEvents/python/MCEventsConfig.py
+++ b/Tracker/TrackerRecAlgs/MCEvents/python/MCEventsConfig.py
@@ -1,3 +1,6 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
from AthenaConfiguration.ComponentFactory import CompFactory
from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/MyExtrapolationExampleConfig.py b/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/MyExtrapolationExampleConfig.py
index 6c2738946d33457377968f29831b007461af969e..d87ae7f5973b3f69dc4206fbb0cc5830f64344ae 100644
--- a/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/MyExtrapolationExampleConfig.py
+++ b/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/MyExtrapolationExampleConfig.py
@@ -1,3 +1,7 @@
+"""
+ 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
diff --git a/Tracker/TrackerRecAlgs/NoisyStripFinder/python/NoisyStripFinderConfig.py b/Tracker/TrackerRecAlgs/NoisyStripFinder/python/NoisyStripFinderConfig.py
index 5e94a1fb0e944666f4f245e8f47bc1100fcd6929..b408cb0e0d8c255833b5e51d47dafe461012a48b 100644
--- a/Tracker/TrackerRecAlgs/NoisyStripFinder/python/NoisyStripFinderConfig.py
+++ b/Tracker/TrackerRecAlgs/NoisyStripFinder/python/NoisyStripFinderConfig.py
@@ -1,5 +1,5 @@
"""
-Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/Tracker/TrackerRecAlgs/NoisyStripFinder/test/NoisyStripFinderDbg.py b/Tracker/TrackerRecAlgs/NoisyStripFinder/test/NoisyStripFinderDbg.py
new file mode 100644
index 0000000000000000000000000000000000000000..b3c0697cfcef6539190aeeebf26e16b4e0d3b100
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/NoisyStripFinder/test/NoisyStripFinderDbg.py
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+"""
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+"""
+
+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 CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from NoisyStripFinder.NoisyStripFinderConfig import NoisyStripFinderCfg
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+run = 1792
+ConfigFlags.Input.Files = [f"/home/tboeckh/Documents/data/TI12/Faser-Physics-00{run}-00000.raw"]
+ConfigFlags.Output.ESDFileName = f"run00{run}-00000.ESD.pool.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+ConfigFlags.IOVDb.DatabaseInstance = "OFLP200"
+ConfigFlags.Input.ProjectName = "data21"
+ConfigFlags.Input.isMC = False
+ConfigFlags.GeoModel.FaserVersion = "FASER-01"
+ConfigFlags.Common.isOnline = False
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+ConfigFlags.Detector.GeometryFaserSCT = True
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolWriteCfg(ConfigFlags))
+acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags))
+acc.merge(FaserSCT_ClusterizationCfg(
+ ConfigFlags,
+ name="LevelClustering",
+ DataObjectName="SCT_LEVELMODE_RDOs",
+ ClusterToolTimingPattern="X1X"))
+acc.merge(NoisyStripFinderCfg(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
+
+sc = acc.run(maxEvents=-1)
+sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/TrackCounts/python/TrackCountsConfig.py b/Tracker/TrackerRecAlgs/TrackCounts/python/TrackCountsConfig.py
index 934579044bcf3cd9d731f4199ed374cda045aebf..da136d0381b211477497550e2072cae6fec029d9 100644
--- a/Tracker/TrackerRecAlgs/TrackCounts/python/TrackCountsConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackCounts/python/TrackCountsConfig.py
@@ -1,5 +1,5 @@
"""
-Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/Tracker/TrackerRecAlgs/TrackSeedPerformanceWriter/python/TrackSeedPerformanceWriterConfig.py b/Tracker/TrackerRecAlgs/TrackSeedPerformanceWriter/python/TrackSeedPerformanceWriterConfig.py
index f8b0eca720b42ed668d7918bac86569007f4d68f..7e064cf0c8a8247426d4f7d613d33e6b8c88626e 100644
--- a/Tracker/TrackerRecAlgs/TrackSeedPerformanceWriter/python/TrackSeedPerformanceWriterConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackSeedPerformanceWriter/python/TrackSeedPerformanceWriterConfig.py
@@ -1,5 +1,5 @@
"""
-Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/Tracker/TrackerRecAlgs/TrackerClusterFit/python/TrackerClusterFitConfig.py b/Tracker/TrackerRecAlgs/TrackerClusterFit/python/TrackerClusterFitConfig.py
index 19fb699dae47f023b5a20665461f2cfe3a3510d6..61d1927154c731009a78be77a8c3b2431ae00871 100644
--- a/Tracker/TrackerRecAlgs/TrackerClusterFit/python/TrackerClusterFitConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackerClusterFit/python/TrackerClusterFitConfig.py
@@ -1,6 +1,5 @@
"""Define methods to construct configured SCT Cluster Fit tools and algorithms
-
-Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentFactory import CompFactory
from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
diff --git a/Tracker/TrackerRecAlgs/TrackerData/python/TrackerSegmentFitDataConfig.py b/Tracker/TrackerRecAlgs/TrackerData/python/TrackerSegmentFitDataConfig.py
index 74497753353d66363d6c09ba12cb9a4c899351ad..cb46c5d7a8b748375ffcdca1178b2aa8ef3bdd3e 100644
--- a/Tracker/TrackerRecAlgs/TrackerData/python/TrackerSegmentFitDataConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackerData/python/TrackerSegmentFitDataConfig.py
@@ -1,5 +1,5 @@
"""
-Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx
index 6902a70138cdad5f8957f0f7e0886d36e6d8213e..e78f7a8c5b72a0ad06d4b360388b3a35f995392b 100755
--- a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx
+++ b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx
@@ -1,6 +1,6 @@
/*
- Copyright (C) 2002-2021 ChRN for the benefit of the ATLAS collaboration
- */
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+*/
/***************************************************************************
-------------------
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
index 1ac1f94d55b030b729a1fc1f000142f275d127fb..f0de1582f2ccf47270203d4d1d6eb9a6bca341e9 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
@@ -1,6 +1,5 @@
"""Define methods to construct configured SCT Cluster Fit tools and algorithms
-
-Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentFactory import CompFactory
from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
index f12ccaa7d336bce3f9be4f0758d4ca06e1b878b9..2481c6f2eeb0109ac08c55a6e133e341e6645268 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
@@ -408,6 +408,7 @@ SegmentFitAlg::checkFit(const Tracker::SegmentFitAlg::fitInfo& seed) const
return false;
}
if (m_tanThetaCut > 0.0 && pow(seed.fitParams(2),2.0) + pow(seed.fitParams(3), 2.0) > m_tanThetaCut * m_tanThetaCut) return false;
+ if (m_tanThetaXZCut > 0.0 && pow(seed.fitParams(2), 2.0) > m_tanThetaXZCut * m_tanThetaXZCut) return false;
for (auto i : seed.candidates)
{
double z = seed.clusters[i]->cluster.globalPosition().z();
@@ -518,7 +519,7 @@ SegmentFitAlg::selectFits(const FitMap& fits) const
// smaller than m_minClustersPerFit
info.remove_if([&](std::shared_ptr<fitInfo> p) {return
((p->clusterMask & selected->clusterMask).count() > m_sharedHitFraction * p->clusterMask.count())
- || (p->clusterMask.count() < m_minClustersPerFit) ;}
+ || ((selectedFits.size() >= 4) && (p->clusterMask.count() < m_minClustersPerFit)) ;}
);
}
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
index bce869119f8a4927ddc272d885aa3cffe6c46ea1..068e1eeab0f2280d38227fd6d3c376c0ffb49e18 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
@@ -413,6 +413,7 @@ class SegmentFitAlg : public AthReentrantAlgorithm
DoubleProperty m_waferTolerance { this, "WaferTolerance", 3 * Gaudi::Units::mm, "Wafer boundary tolerance for extrapolated fit." };
DoubleProperty m_yResidualCut { this, "ResidualCut", 3 * Gaudi::Units::mm, "Cluster y tolerance for compatibility with extrapolated fit." };
DoubleProperty m_tanThetaCut { this, "TanThetaCut", 0.0, "Maximum accepted tangent of track angle from beam axis; 0 means no cut."};
+ DoubleProperty m_tanThetaXZCut { this, "TanThetaXZCut", 0.0, "Maximum accepted tangent of track angle from beam axis in x direction; 0 means no cut."};
DoubleProperty m_reducedChi2Cut { this, "ReducedChi2Cut", 10.0, "Maximum accepted chi^2 per degree of freedom for final fits; 0 means no cut." };
DoubleProperty m_sharedHitFraction { this, "SharedHitFraction", -1., "Fraction of hits which are allowed to be shared between two fits." };
UnsignedIntegerProperty m_minClustersPerFit { this, "MinClustersPerFit", 4, "Minimum number of clusters a fit has to have." };
diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/python/TrackerTruthConfig.py b/Tracker/TrackerRecAlgs/TrackerTruth/python/TrackerTruthConfig.py
index 475e7ba6f60844cdea8c4b5da0ef4f43da70d2f5..92f56a8edc34c173636f16fb9162c8ab17d8850f 100644
--- a/Tracker/TrackerRecAlgs/TrackerTruth/python/TrackerTruthConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackerTruth/python/TrackerTruthConfig.py
@@ -1,5 +1,5 @@
"""
-Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py b/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py
index 6466848fb515373c2b33bd8c0978d4f162e053bb..095edc7648474af1510d8c1f1174f8c955eccaa8 100644
--- a/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py
+++ b/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py
@@ -1,4 +1,6 @@
-# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS and FASER collaborations
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
###############################################################
#
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
index df260a0f6e8da85b268940795f33775f54b60d5b..f8803a279be6565afe9597c701d5bf3f221fef94 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
@@ -39,9 +39,11 @@ FaserActs::CuboidVolumeBuilder::Config FaserActsLayerBuilder::buildVolume(const
FaserActs::CuboidVolumeBuilder::Config cvbConfig;
//std::vector<Acts::CuboidVolumeBuilder::VolumeConfig> volumeConfigs = {};
std::vector<FaserActs::CuboidVolumeBuilder::VolumeConfig> volumeConfigs = {};
-
- for (int iStation=1; iStation<4; iStation++) {
+ // iStation starts from 1 for the main detector (Faser-01 geometry) and from 0 if we include the IFT (Faser-02 geometry)
+ auto siDetMng = static_cast<const TrackerDD::SCT_DetectorManager*>(m_cfg.mng);
+ int nStations = siDetMng->numerology().numStations();
+ for (int iStation=4-nStations; iStation<4; iStation++) {
//Acts::CuboidVolumeBuilder::VolumeConfig volumeConfig;
FaserActs::CuboidVolumeBuilder::VolumeConfig volumeConfig;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
index c5c0c5e454ab0b43622689bda4f7c3691ee1abca..22f075abb55848488419bc7129574a6a6a6741fe 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
+++ b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
@@ -27,6 +27,7 @@ atlas_add_library( FaserActsKalmanFilterLib
atlas_add_component(FaserActsKalmanFilter
FaserActsKalmanFilter/ActsTrackSeedTool.h
FaserActsKalmanFilter/CircleFit.h
+ FaserActsKalmanFilter/CircleFitTrackSeedTool.h
FaserActsKalmanFilter/CKF2.h
FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h
FaserActsKalmanFilter/EffPlotTool.h
@@ -34,12 +35,14 @@ atlas_add_component(FaserActsKalmanFilter
FaserActsKalmanFilter/FaserActsGeometryContainers.h
FaserActsKalmanFilter/FaserActsKalmanFilterAlg.h
FaserActsKalmanFilter/FaserActsRecMultiTrajectory.h
+ FaserActsKalmanFilter/GhostBusters.h
FaserActsKalmanFilter/MyTrackSeedTool.h
FaserActsKalmanFilter/IdentifierLink.h
FaserActsKalmanFilter/IndexSourceLink.h
FaserActsKalmanFilter/ITrackFinderTool.h
FaserActsKalmanFilter/ITrackSeedTool.h
FaserActsKalmanFilter/KalmanFitterTool.h
+ FaserActsKalmanFilter/LinearFit.h
# FaserActsKalmanFilter/ClusterTrackSeedTool.h
# FaserActsKalmanFilter/TruthTrackFinderTool.h
FaserActsKalmanFilter/Measurement.h
@@ -49,6 +52,7 @@ atlas_add_component(FaserActsKalmanFilter
FaserActsKalmanFilter/ResPlotTool.h
FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h
+ FaserActsKalmanFilter/SeedingAlg.h
# FaserActsKalmanFilter/SegmentFitClusterTrackFinderTool.h
# FaserActsKalmanFilter/SegmentFitTrackFinderTool.h
FaserActsKalmanFilter/SimWriterTool.h
@@ -64,17 +68,21 @@ atlas_add_component(FaserActsKalmanFilter
# FaserActsKalmanFilter/TruthSeededTrackFinderTool.h
FaserActsKalmanFilter/ThreeStationTrackSeedTool.h
src/ActsTrackSeedTool.cxx
+ src/CircleFit.cxx
+ src/CircleFitTrackSeedTool.cxx
src/CKF2.cxx
# src/ClusterTrackSeedTool.cxx
src/CombinatorialKalmanFilterAlg.cxx
src/EffPlotTool.cxx
src/FaserActsKalmanFilterAlg.cxx
+ src/GhostBusters.cxx
src/MyTrackSeedTool.cxx
src/KalmanFitterTool.cxx
# src/MultiTrackFinderTool.cxx
src/PerformanceWriterTool.cxx
src/PlotHelpers.cxx
src/ResPlotTool.cxx
+ src/SeedingAlg.cxx
src/RootTrajectoryStatesWriterTool.cxx
src/RootTrajectorySummaryWriterTool.cxx
# src/SegmentFitClusterTrackFinderTool.cxx
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ActsTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ActsTrackSeedTool.h
index fe129a47d2867405ece6b807bf56803f814a2882..95311206b6300722344935b535c4f0d963728392 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ActsTrackSeedTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ActsTrackSeedTool.h
@@ -28,6 +28,7 @@ public:
const std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters() const override;
const std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters() const override;
const std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints() const override;
+ double targetZPosition() const override;
private:
std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>> m_initialTrackParameters;
@@ -38,6 +39,7 @@ private:
std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> m_clusters {};
std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> m_seedClusters {};
std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> m_spacePoints {};
+ double m_targetZPosition {0.};
const FaserSCT_ID* m_idHelper {nullptr};
const TrackerDD::SCT_DetectorManager* m_detManager {nullptr};
@@ -107,4 +109,8 @@ ActsTrackSeedTool::spacePoints() const {
return m_spacePoints;
}
+inline double ActsTrackSeedTool::targetZPosition() const {
+ return m_targetZPosition;
+}
+
#endif // FASERACTSKALMANFILTER_ACTSTRACKSEEDTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
index 7e3b0e76af52738963d42f7cd716a1c32644ca47..6209dda840a81c780b8ec4b71d745165e9bc1ce3 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
@@ -131,7 +131,6 @@ private:
Gaudi::Property<double> m_maxSteps {this, "maxSteps", 10000};
Gaudi::Property<double> m_chi2Max {this, "chi2Max", 15};
Gaudi::Property<unsigned long> m_nMax {this, "nMax", 10};
- Gaudi::Property<size_t> m_nTrajectories {this, "nTrajectories", 2};
SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
ToolHandle<ITrackSeedTool> m_trackSeedTool {this, "TrackSeed", "ClusterTrackSeedTool"};
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
@@ -140,6 +139,7 @@ private:
ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"};
ToolHandle<KalmanFitterTool> m_kalmanFitterTool1 {this, "KalmanFitterTool1", "KalmanFitterTool"};
ToolHandle<KalmanFitterTool> m_kalmanFitterTool2 {this, "KalmanFitterTool2", "KalmanFitterTool"};
+ Gaudi::Property<bool> m_isMC {this, "isMC", false};
std::unique_ptr<Trk::Track> makeTrack(Acts::GeometryContext& tgContext, TrackFitterResult& fitResult) const;
std::unique_ptr<Trk::Track> makeTrack(const Acts::GeometryContext &geoCtx, const FaserActsRecMultiTrajectory &traj) const;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFit.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFit.h
index 664b7005ddfdb31998dc3a5665fdd871a9d8ce31..780a476a722707f280b427a1678f16fcfb5b1cdf 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFit.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFit.h
@@ -5,6 +5,25 @@
#include <cmath>
#include <vector>
+namespace CircleFit {
+
+
+class CircleData {
+public:
+ CircleData(const std::vector<const Tracker::FaserSCT_SpacePoint *> &spacePoints);
+ CircleData(const std::vector<Amg::Vector3D> &spacePoints);
+ double meanX() const;
+ double meanY() const;
+ double x(int i) const;
+ double y(int i) const;
+ int size() const;
+
+private:
+ std::vector<double> m_x{};
+ std::vector<double> m_y{};
+ int m_size = 0;
+};
+
struct Circle {
public:
@@ -19,170 +38,9 @@ public:
double j = 0;
};
+double sigma(const CircleData &data, const Circle &circle);
+Circle circleFit(const CircleData &data);
-class CircleData {
-public:
- CircleData(const std::vector<const Tracker::FaserSCT_SpacePoint*> &spacePoints) {
- for (auto sp : spacePoints) {
- m_x.push_back(sp->globalPosition().z());
- m_y.push_back(sp->globalPosition().y());
- }
- m_size = spacePoints.size();
- }
-
- double meanX() const {
- double mean = 0;
- for (double i : m_x) mean += i;
- return mean / m_size;
- }
-
- double meanY() const {
- double mean = 0;
- for (double i : m_y) mean += i;
- return mean / m_size;
- }
-
- double x(int i) const {
- return m_x[i];
- }
-
- double y(int i) const {
- return m_y[i];
- }
-
- int size() const {
- return m_size;
- }
-
-private:
- std::vector<double> m_x {};
- std::vector<double> m_y {};
- int m_size = 0;
-};
-
-
-double Sigma (CircleData& data, Circle& circle) {
- double sum=0.,dx,dy;
-
- for (int i=0; i<data.size(); i++) {
- dx = data.x(i) - circle.cx;
- dy = data.y(i) - circle.cy;
- sum += (sqrt(dx*dx+dy*dy) - circle.r) * (sqrt(dx*dx+dy*dy) - circle.r);
- }
- return sqrt(sum/data.size());
-}
-
-
-Circle CircleFit(CircleData& data)
-/*
- Circle fit to a given set of data points (in 2D)
-
- This is an algebraic fit, due to Taubin, based on the journal article
-
- G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
- Space Curves Defined By Implicit Equations, With
- Applications To Edge And Range Image Segmentation",
- IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
-
- Input: data - the class of data (contains the given points):
-
- data.size() - the number of data points
- data.X[] - the array of X-coordinates
- data.Y[] - the array of Y-coordinates
-
- Output:
- circle - parameters of the fitting circle:
-
- circle.a - the X-coordinate of the center of the fitting circle
- circle.b - the Y-coordinate of the center of the fitting circle
- circle.r - the radius of the fitting circle
- circle.s - the root mean square error (the estimate of sigma)
- circle.j - the total number of iterations
-
- The method is based on the minimization of the function
-
- sum [(x-a)^2 + (y-b)^2 - R^2]^2
- F = -------------------------------
- sum [(x-a)^2 + (y-b)^2]
-
- This method is more balanced than the simple Kasa fit.
-
- It works well whether data points are sampled along an entire circle or
- along a small arc.
-
- It still has a small bias and its statistical accuracy is slightly
- lower than that of the geometric fit (minimizing geometric distances),
- but slightly higher than that of the very similar Pratt fit.
- Besides, the Taubin fit is slightly simpler than the Pratt fit
-
- It provides a very good initial guess for a subsequent geometric fit.
-
- Nikolai Chernov (September 2012)
-
-*/
-{
- int i,iter,IterMAX=99;
-
- double Xi,Yi,Zi;
- double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
- double A0,A1,A2,A22,A3,A33;
- double Dy,xnew,x,ynew,y;
- double DET,Xcenter,Ycenter;
-
- Circle circle;
-
- Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
- for (i=0; i<data.size(); i++) {
- Xi = data.x(i) - data.meanX();
- Yi = data.y(i) - data.meanY();
- Zi = Xi*Xi + Yi*Yi;
-
- Mxy += Xi*Yi;
- Mxx += Xi*Xi;
- Myy += Yi*Yi;
- Mxz += Xi*Zi;
- Myz += Yi*Zi;
- Mzz += Zi*Zi;
- }
- Mxx /= data.size();
- Myy /= data.size();
- Mxy /= data.size();
- Mxz /= data.size();
- Myz /= data.size();
- Mzz /= data.size();
-
- Mz = Mxx + Myy;
- Cov_xy = Mxx*Myy - Mxy*Mxy;
- Var_z = Mzz - Mz*Mz;
- A3 = 4*Mz;
- A2 = -3*Mz*Mz - Mzz;
- A1 = Var_z*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;
- A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;
- A22 = A2 + A2;
- A33 = A3 + A3 + A3;
-
- for (x=0.,y=A0,iter=0; iter<IterMAX; iter++) {
- Dy = A1 + x*(A22 + A33*x);
- xnew = x - y/Dy;
- if ((xnew == x)||(!std::isfinite(xnew))) break;
- ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
- if (abs(ynew)>=abs(y)) break;
- x = xnew; y = ynew;
- }
-
- DET = x*x - x*Mz + Cov_xy;
- Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/2;
- Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/2;
-
- circle.cx = Xcenter + data.meanX();
- circle.cy = Ycenter + data.meanY();
- circle.r = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
- circle.s = Sigma(data,circle);
- circle.i = 0;
- circle.j = iter;
-
- return circle;
-}
-
+} // namespace CircleFit
-#endif // FASERACTSKALMANFILTER_CIRCLEFIT_H
\ No newline at end of file
+#endif // FASERACTSKALMANFILTER_CIRCLEFIT_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFitTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFitTrackSeedTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..368abb1f228e45e69d0c38da02e63dc65c15fb4d
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFitTrackSeedTool.h
@@ -0,0 +1,165 @@
+#ifndef FASERACTSKALMANFILTER_CIRCLEFITTRACKSEEDTOOL_H
+#define FASERACTSKALMANFILTER_CIRCLEFITTRACKSEEDTOOL_H
+
+#include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "Gaudi/Property.h"
+#include "GaudiKernel/IInterface.h"
+#include "GaudiKernel/StatusCode.h"
+
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+#include "FaserActsKalmanFilter/ITrackSeedTool.h"
+#include "TrkTrack/TrackCollection.h"
+// #include "TrackerSimData/TrackerSimDataCollection.h"
+// #include "GeneratorObjects/McEventCollection.h"
+#include <array>
+#include <memory>
+#include <string>
+#include <vector>
+#include <boost/dynamic_bitset.hpp>
+
+typedef boost::dynamic_bitset<> ClusterSet;
+
+class FaserSCT_ID;
+namespace TrackerDD { class SCT_DetectorManager; }
+
+class CircleFitTrackSeedTool : public extends<AthAlgTool, ITrackSeedTool> {
+public:
+ CircleFitTrackSeedTool(const std::string& type, const std::string& name, const IInterface* parent);
+ virtual ~CircleFitTrackSeedTool() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+ virtual StatusCode run() override;
+
+ const std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>> initialTrackParameters() const override;
+ const std::shared_ptr<const Acts::Surface> initialSurface() const override;
+ const std::shared_ptr<std::vector<IndexSourceLink>> sourceLinks() const override;
+ const std::shared_ptr<IdentifierLink> idLinks() const override;
+ const std::shared_ptr<std::vector<Measurement>> measurements() const override;
+ const std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters() const override;
+ const std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters() const override;
+ const std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints() const override;
+ double targetZPosition() const override;
+
+private:
+ struct Segment {
+ public:
+ Segment(const Trk::Track* track, const FaserSCT_ID *idHelper);
+ int station;
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters;
+ std::vector<const Tracker::FaserSCT_SpacePoint*> spacePoints;
+ ClusterSet clusterSet;
+ Acts::Vector3 position;
+ std::vector<Acts::Vector3> fakePositions;
+ Acts::Vector3 momentum;
+ };
+
+ struct Seed {
+ Seed(const std::vector<Segment> &segments);
+
+ ClusterSet clusterSet;
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters;
+ std::vector<const Tracker::FaserSCT_SpacePoint*> spacePoints;
+ std::vector<Acts::Vector3> positions;
+ std::vector<Acts::Vector3> fakePositions;
+
+ double c0, c1, cx, cy, r, chi2, momentum, charge, minZ;
+ Acts::Vector3 direction;
+ size_t size;
+ Acts::CurvilinearTrackParameters get_params(double origin, Acts::BoundSymMatrix cov) const;
+
+ private:
+ void getChi2();
+ double getX(double z);
+ double getY(double z);
+ void fit();
+ void fakeFit(double B=0.55);
+ void linearFit(const std::vector<Acts::Vector2> &points);
+ double m_dx, m_dy;
+ double m_MeV2GeV = 1e-3;
+ double m_sigma_x = 0.8;
+ double m_sigma_y = 0.016;
+ };
+
+ void go(const std::array<std::vector<Segment>, 4> &v, std::vector<Segment> &combination,
+ std::vector<Seed> &seeds, int offset, int k);
+
+ static std::map<Identifier, Index> s_indexMap;
+ static std::map<Identifier, const Tracker::FaserSCT_SpacePoint*> s_spacePointMap;
+
+ std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>> m_initialTrackParameters;
+ std::shared_ptr<const Acts::Surface> m_initialSurface;
+ std::shared_ptr<std::vector<IndexSourceLink>> m_sourceLinks {};
+ std::shared_ptr<IdentifierLink> m_idLinks {};
+ std::shared_ptr<std::vector<Measurement>> m_measurements {};
+ std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> m_clusters {};
+ std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> m_seedClusters {};
+ std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> m_spacePoints {};
+ double m_targetZPosition {0};
+
+ const FaserSCT_ID* m_idHelper {nullptr};
+ const TrackerDD::SCT_DetectorManager* m_detManager {nullptr};
+
+ ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool { this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "SegmentFit", "Input track collection name" };
+ SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainerKey { this, "ClusterContainer", "SCT_ClusterContainer"};
+ SG::ReadHandleKey<FaserSCT_SpacePointContainer> m_spacePointContainerKey { this, "SpacePoints", "SCT_SpacePointContainer"};
+ // SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey { this, "McEventCollection", "TruthEvent"};
+ // SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollectionKey { this, "TrackerSimDataCollection", "SCT_SDO_Map"};
+
+ Gaudi::Property<double> m_std_cluster {this, "std_cluster", 0.04};
+ // covariance of the initial parameters
+ Gaudi::Property<double> m_covLoc0 {this, "covLoc0", 1};
+ Gaudi::Property<double> m_covLoc1 {this, "covLoc1", 1};
+ Gaudi::Property<double> m_covPhi {this, "covPhi", 1};
+ Gaudi::Property<double> m_covTheta {this, "covTheta", 1};
+ Gaudi::Property<double> m_covQOverP {this, "covQOverP", 1};
+ Gaudi::Property<double> m_covTime {this, "covTime", 1};
+};
+
+inline const std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>>
+CircleFitTrackSeedTool::initialTrackParameters() const {
+ return m_initialTrackParameters;
+}
+
+inline const std::shared_ptr<const Acts::Surface>
+CircleFitTrackSeedTool::initialSurface() const {
+ return m_initialSurface;
+}
+
+inline const std::shared_ptr<std::vector<IndexSourceLink>>
+CircleFitTrackSeedTool::sourceLinks() const {
+ return m_sourceLinks;
+}
+
+inline const std::shared_ptr<IdentifierLink>
+CircleFitTrackSeedTool::idLinks() const {
+ return m_idLinks;
+}
+
+inline const std::shared_ptr<std::vector<Measurement>>
+CircleFitTrackSeedTool::measurements() const {
+ return m_measurements;
+}
+
+inline const std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>>
+CircleFitTrackSeedTool::clusters() const {
+ return m_clusters;
+}
+
+inline const std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>>
+CircleFitTrackSeedTool::seedClusters() const {
+ return m_seedClusters;
+}
+
+inline const std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>>
+CircleFitTrackSeedTool::spacePoints() const {
+ return m_spacePoints;
+}
+
+inline double CircleFitTrackSeedTool::targetZPosition() const {
+ return m_targetZPosition;
+}
+
+#endif // FASERACTSKALMANFILTER_CIRCLEFITTRACKSEEDTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h
index ba983f245fdce513858c27cb0a814b48313eac0b..28058887c9b6f05d384962206c8cca1a64a9d66c 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h
@@ -114,6 +114,7 @@ public:
Gaudi::Property<double> m_maxSteps {this, "maxSteps", 1000};
Gaudi::Property<double> m_chi2Max {this, "chi2Max", 15};
Gaudi::Property<unsigned long> m_nMax {this, "nMax", 10};
+ Gaudi::Property<bool> m_isMC {this, "isMC", false};
SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
ToolHandle<ITrackSeedTool> m_trackSeedTool {this, "TrackSeed", "ClusterTrackSeedTool"};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h
new file mode 100644
index 0000000000000000000000000000000000000000..f8ac55d6cbb62c20ac2de7ce9c8b5741a2bfe8cb
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h
@@ -0,0 +1,75 @@
+#ifndef FASERACTSKALMANFILTER_GHOSTBUSTERS_H
+#define FASERACTSKALMANFILTER_GHOSTBUSTERS_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthHistogramming.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrackerSimData/TrackerSimDataCollection.h"
+#include "FaserActsKalmanFilter/TrackClassification.h"
+
+
+class TTree;
+class FaserSCT_ID;
+
+class GhostBusters : public AthReentrantAlgorithm, AthHistogramming {
+public:
+ GhostBusters(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~GhostBusters() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+ const ServiceHandle <ITHistSvc> &histSvc() const;
+
+private:
+ struct Segment {
+ public:
+ Segment(const Trk::Track *track, const FaserSCT_ID *idHelper);
+ // Segment(const Trk::Track *track, const TrackerSimDataCollection &simDataCollection,
+ // const FaserSCT_ID *idHelper);
+ inline int station() const { return m_station; }
+ inline Amg::Vector3D position() const { return m_position; }
+ inline double x() const { return m_position.x(); }
+ inline double y() const { return m_position.y(); }
+ inline double z() const { return m_position.z(); }
+ inline double chi2() const { return m_chi2; }
+ inline size_t size() const { return m_size; }
+ inline bool isGhost() const { return m_isGhost; }
+ // inline size_t majorityHits() const { return m_particleHitCounts.empty() ? -1 : m_particleHitCounts.front().hitCount; }
+ inline const Trk::Track *track() const { return m_track; }
+ inline void setGhost() { m_isGhost = true; }
+ private:
+ // std::vector<ParticleHitCount> m_particleHitCounts {};
+ size_t m_size;
+ Amg::Vector3D m_position;
+ int m_station;
+ double m_chi2;
+ bool m_isGhost = false;
+ const Trk::Track *m_track;
+ };
+
+ ServiceHandle <ITHistSvc> m_histSvc;
+ SG::WriteHandleKey<TrackCollection> m_outputTrackCollection { this, "OutputCollection", "Segments", "Output track collection name" };
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "SegmentFit", "Input track collection name" };
+ // SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollectionKey { this, "TrackerSimDataCollection", "SCT_SDO_Map"};
+ DoubleProperty m_xTolerance {this, "xTolerance", 0.5};
+ DoubleProperty m_yTolerance {this, "yTolerance", 0.25};
+
+ const FaserSCT_ID *m_idHelper;
+
+ mutable TTree* m_tree;
+ mutable unsigned int m_event_number;
+ mutable double m_x;
+ mutable double m_y;
+ mutable double m_z;
+ mutable double m_chi2;
+ mutable unsigned int m_station;
+ // mutable unsigned int m_majorityHits;
+ mutable unsigned int m_size;
+ mutable bool m_isGhost;
+};
+
+inline const ServiceHandle <ITHistSvc> &GhostBusters::histSvc() const {
+ return m_histSvc;
+}
+
+#endif // FASERACTSKALMANFILTER_GHOSTBUSTERS_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ITrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ITrackSeedTool.h
index 8345823daedd83b2ce1220b843303492f3101a57..0b664f78ec6e60f1d75f504a343221970d0b938d 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ITrackSeedTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ITrackSeedTool.h
@@ -23,6 +23,7 @@ public:
virtual const std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters() const = 0;
virtual const std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters() const = 0;
virtual const std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints() const = 0;
+ virtual double targetZPosition() const = 0;
};
#endif // FASERACTSKALMANFILTER_ITRACKSEEDTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
index 9c26faeb3f354618133ae478fca9a4f3cea947bd..6b8f8471c1c15866b58b45bf0b2c9ede581f0312 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
@@ -42,14 +42,14 @@ public:
virtual Acts::MagneticFieldContext getMagneticFieldContext(const EventContext& ctx) const;
std::unique_ptr<Trk::Track> fit(const EventContext& ctx, const Trk::Track& inputTrack,
std::vector<FaserActsRecMultiTrajectory> &trajectories,
- const Acts::BoundVector& inputVector) const;
+ const Acts::BoundVector& inputVector, bool isMC, double origin) const;
private:
const FaserSCT_ID* m_idHelper {nullptr};
std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
getMeasurementsFromTrack(const Trk::Track &track) const;
// Acts::BoundTrackParameters getParametersFromTrack(const Acts::BoundVector& params, const Trk::TrackParameters *inputParameters) const;
- Acts::BoundTrackParameters getParametersFromTrack(const Trk::TrackParameters *inputParameters, const Acts::BoundVector& inputVector) const;
+ Acts::BoundTrackParameters getParametersFromTrack(const Trk::TrackParameters *inputParameters, const Acts::BoundVector& inputVector, double origin) const;
std::shared_ptr<TrackFitterFunction> m_fit;
std::unique_ptr<const Acts::Logger> m_logger;
Gaudi::Property<std::string> m_actsLogging {this, "ActsLogging", "VERBOSE"};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/LinearFit.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/LinearFit.h
new file mode 100644
index 0000000000000000000000000000000000000000..f367110a4b7a05cb9bc610628995cfa6df53ddb0
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/LinearFit.h
@@ -0,0 +1,26 @@
+#ifndef FASERACTSKALMANFILTER_LINEARFIT_H
+#define FASERACTSKALMANFILTER_LINEARFIT_H
+
+#include "Acts/Definitions/Algebra.hpp"
+
+namespace LinearFit {
+
+using Vector2 = Acts::Vector2;
+
+std::pair <Vector2, Vector2 > linearFit(const std::vector<Vector2> &points) {
+ size_t nPoints = points.size();
+ Eigen::Matrix< Vector2::Scalar, Eigen::Dynamic, Eigen::Dynamic > centers(nPoints, 2);
+ for (size_t i = 0; i < nPoints; ++i) centers.row(i) = points[i];
+
+ Vector2 origin = centers.colwise().mean();
+ Eigen::MatrixXd centered = centers.rowwise() - origin.transpose();
+ Eigen::MatrixXd cov = centered.adjoint() * centered;
+ Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
+ Vector2 axis = eig.eigenvectors().col(1).normalized();
+
+ return std::make_pair(origin, axis);
+}
+
+} // namespace LinearFit
+
+#endif // FASERACTSKALMANFILTER_LINEARFIT_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/MyTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/MyTrackSeedTool.h
index bd29bbf5f020fe5ce139c8cbd12eb52154511110..abc4e987434c536b37a20862bef87fc5fe8cbbac 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/MyTrackSeedTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/MyTrackSeedTool.h
@@ -35,6 +35,7 @@ public:
const std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters() const override;
const std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters() const override;
const std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints() const override;
+ double targetZPosition() const override;
private:
std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>> m_initialTrackParameters;
@@ -45,6 +46,7 @@ private:
std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> m_clusters {};
std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> m_seedClusters {};
std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> m_spacePoints {};
+ double m_targetZPosition {0.};
const FaserSCT_ID* m_idHelper {nullptr};
const TrackerDD::SCT_DetectorManager* m_detManager {nullptr};
@@ -114,4 +116,8 @@ MyTrackSeedTool::spacePoints() const {
return m_spacePoints;
}
+inline double MyTrackSeedTool::targetZPosition() const {
+ return m_targetZPosition;
+}
+
#endif // FASERACTSKALMANFILTER_MYTRACKSEEDTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
index ee8ef5d25a9f808d4ed8ae87787b437188969b49..6f0885d33ec1200edbd28761352a36840fd11e47 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
@@ -29,7 +29,7 @@ public:
StatusCode initialize() override;
StatusCode finalize() override;
- StatusCode write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories) const;
+ StatusCode write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories, bool isMC) const;
private:
SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey {this, "McEventCollection", "TruthEvent"};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h
index 33b5ee882b8a5f4358cb8b65203704f4bbd9c9a3..90cf6cf84ca8a95b41397bd00dea30d85135967f 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h
@@ -25,7 +25,7 @@ class RootTrajectorySummaryWriterTool : public AthAlgTool {
public:
RootTrajectorySummaryWriterTool(const std::string& type, const std::string& name, const IInterface* parent);
~RootTrajectorySummaryWriterTool() override = default;
- StatusCode write(const Acts::GeometryContext& geoContext, const TrajectoriesContainer& trajectories) const;
+ StatusCode write( const Acts::GeometryContext& geoContext, const TrajectoriesContainer& trajectories, bool isMC) const;
StatusCode initialize() override;
StatusCode finalize() override;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SeedingAlg.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SeedingAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..49fd1c901ab2f0159cd033401ddde589b61262e6
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SeedingAlg.h
@@ -0,0 +1,23 @@
+#ifndef FASERACTSKALMANFILTER_SEEDINGALG_H
+#define FASERACTSKALMANFILTER_SEEDINGALG_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "FaserActsKalmanFilter/ITrackSeedTool.h"
+#include <boost/dynamic_bitset.hpp>
+
+
+class SeedingAlg : public AthAlgorithm {
+public:
+ SeedingAlg(const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~SeedingAlg() = default;
+
+ StatusCode initialize() override;
+ StatusCode execute() override;
+ StatusCode finalize() override;
+
+private:
+ ToolHandle<ITrackSeedTool> m_trackSeedTool {this, "TrackSeed", "CircleFitTrackSeedTool"};
+};
+
+#endif // FASERACTSKALMANFILTER_SEEDINGALG_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ThreeStationTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ThreeStationTrackSeedTool.h
index 3cf450414f20733b873ceda542f7c4ab44cdd573..22cbf9af8420d5100635bbc0f6a00f760e111097 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ThreeStationTrackSeedTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ThreeStationTrackSeedTool.h
@@ -33,6 +33,7 @@ public:
const std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters() const override;
const std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters() const override;
const std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints() const override;
+ double targetZPosition() const override;
private:
@@ -56,6 +57,7 @@ private:
std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> m_clusters {};
std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> m_seedClusters {};
std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> m_spacePoints {};
+ double m_targetZPosition {0.};
const FaserSCT_ID* m_idHelper {nullptr};
const TrackerDD::SCT_DetectorManager* m_detManager {nullptr};
@@ -125,4 +127,8 @@ ThreeStationTrackSeedTool::spacePoints() const {
return m_spacePoints;
}
+inline double ThreeStationTrackSeedTool::targetZPosition() const {
+ return m_targetZPosition;
+}
+
#endif // FASERACTSKALMANFILTER_THREESTATIONTRACKSEEDTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/TrackClassification.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/TrackClassification.h
index 790cdb8f95643d270fe12e6e7a38ba0b6739382e..2ca594c7dfd08c322c08ae193239fa8cafbc5ae9 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/TrackClassification.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/TrackClassification.h
@@ -15,4 +15,9 @@ void identifyContributingParticles(
const FaserActsRecMultiTrajectory& trajectories, size_t tip,
std::vector<ParticleHitCount>& particleHitCounts);
+void identifyContributingParticles(
+ const TrackerSimDataCollection& simDataCollection,
+ const std::vector<const Tracker::FaserSCT_Cluster*> clusters,
+ std::vector<ParticleHitCount>& particleHitCounts);
+
#endif // FASERACTSKALMANFILTER_TRACKCLASSIFICATION_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
index 728cc83e1cf04a37a67486bcaee6662e5dfb73ad..ea2dc0cb9eb08c309bf0453216230fe0b2074b4f 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
@@ -1,4 +1,6 @@
-# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS and FASER collaborations
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
@@ -34,10 +36,11 @@ def CKF2Cfg(flags, **kwargs):
acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
acc.merge(acts_tracking_geometry_svc )
-# track_seed_tool = CompFactory.ClusterTrackSeedTool()
- track_seed_tool = CompFactory.ThreeStationTrackSeedTool()
+ # track_seed_tool = CompFactory.ClusterTrackSeedTool()
# track_seed_tool = CompFactory.ActsTrackSeedTool()
# track_seed_tool = CompFactory.MyTrackSeedTool()
+ # track_seed_tool = CompFactory.ThreeStationTrackSeedTool()
+ track_seed_tool = CompFactory.CircleFitTrackSeedTool()
sigma_loc0 = 1.9e-2
sigma_loc1 = 9e-1
sigma_phi = 3.3e-2
@@ -53,7 +56,7 @@ def CKF2Cfg(flags, **kwargs):
track_seed_tool.covTheta = initial_variance_inflation[3] * sigma_theta * sigma_theta
track_seed_tool.covQOverP = initial_variance_inflation[4] * sigma_qop * sigma_qop
track_seed_tool.std_cluster = 0.0231
- track_seed_tool.origin = 0
+ track_seed_tool.TrackCollection = "Segments"
trajectory_states_writer_tool = CompFactory.RootTrajectoryStatesWriterTool()
trajectory_states_writer_tool.noDiagnostics = kwargs["noDiagnostics"]
@@ -87,8 +90,9 @@ def CKF2Cfg(flags, **kwargs):
kalman_fitter1.noDiagnostics = kwargs["noDiagnostics"]
kalman_fitter1.ActsLogging = "INFO"
kalman_fitter1.SummaryWriter = True
- kalman_fitter1.StatesWriter = True
+ kalman_fitter1.StatesWriter = False
kalman_fitter1.SeedCovarianceScale = 10
+ kalman_fitter1.isMC = flags.Input.isMC
kalman_fitter1.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool1
kalman_fitter1.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool1
ckf.KalmanFitterTool1 = kalman_fitter1
@@ -97,8 +101,9 @@ def CKF2Cfg(flags, **kwargs):
kalman_fitter2.noDiagnostics = kwargs["noDiagnostics"]
kalman_fitter2.ActsLogging = "INFO"
kalman_fitter2.SummaryWriter = True
- kalman_fitter2.StatesWriter = True
- kalman_fitter2.SeedCovarianceScale = 1
+ kalman_fitter2.StatesWriter = False
+ kalman_fitter2.SeedCovarianceScale = 10
+ kalman_fitter2.isMC = flags.Input.isMC
kalman_fitter2.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool2
kalman_fitter2.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool2
ckf.KalmanFitterTool2 = kalman_fitter2
@@ -107,10 +112,14 @@ def CKF2Cfg(flags, **kwargs):
ckf.ActsLogging = "INFO"
ckf.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool
ckf.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool
- ckf.PerformanceWriter = trajectory_performance_writer_tool
+ ckf.PerformanceWriterTool = trajectory_performance_writer_tool
+ ckf.isMC = flags.Input.isMC
+ ckf.SummaryWriter = True
+ ckf.StatesWriter = False
+ ckf.PerformanceWriter = False
ckf.nMax = 10
ckf.chi2Max = 25
acc.addEventAlgo(ckf)
- acc.merge(CKF2_OutputCfg(flags))
+ # acc.merge(CKF2_OutputCfg(flags))
return acc
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/FaserActsKalmanFilterConfig.py b/Tracking/Acts/FaserActsKalmanFilter/python/FaserActsKalmanFilterConfig.py
index cc44658845f4202f7295faf93de273cee776c552..a2ad7705968c54771f72345dfb56c82e9f03a8c4 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/python/FaserActsKalmanFilterConfig.py
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/FaserActsKalmanFilterConfig.py
@@ -1,10 +1,11 @@
-# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS and FASER collaborations
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
-from math import pi
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
-from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+# from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
# FaserActsKalmanFilterAlg,FaserActsTrackingGeometrySvc, FaserActsTrackingGeometryTool, FaserActsExtrapolationTool,FaserActsAlignmentCondAlg, THistSvc = CompFactory.getComps("FaserActsKalmanFilterAlg","FaserActsTrackingGeometrySvc", "FaserActsTrackingGeometryTool", "FaserActsExtrapolationTool","FaserActsAlignmentCondAlg", "THistSvc")
#
@@ -90,7 +91,7 @@ def FaserActsKalmanFilterCfg(flags, **kwargs):
sigma_loc1 = 9e-1
sigma_phi = 3.3e-2
sigma_theta = 2.9e-4
- sigma_p_rel = 0.1
+ # sigma_p_rel = 0.1
p = 1
sigma_qop = 0.1 * p / (p * p)
initial_variance_inflation = [100, 100, 1000, 1000, 1000]
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py b/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..ebe56d179393fcf3f51cb8df510271fbf5a85ebe
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py
@@ -0,0 +1,35 @@
+"""
+ 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 FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometrySvcCfg
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+
+def GhostBusters_OutputCfg(flags):
+ acc = ComponentAccumulator()
+ itemList = ["xAOD::EventInfo#*",
+ "xAOD::EventAuxInfo#*",
+ "TrackCollection#*",
+ ]
+ acc.merge(OutputStreamCfg(flags, "ESD", itemList))
+ ostream = acc.getEventAlgo("OutputStreamESD")
+ ostream.TakeItemsFromInput = True
+ return acc
+
+
+def GhostBustersCfg(flags, **kwargs):
+ acc = FaserSCT_GeometryCfg(flags)
+ acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
+ acc.merge(acts_tracking_geometry_svc )
+ acc.addEventAlgo(CompFactory.GhostBusters(**kwargs))
+ acc.merge(GhostBusters_OutputCfg(flags))
+
+ # thistSvc = CompFactory.THistSvc()
+ # thistSvc.Output += ["HIST2 DATAFILE='GhostBusters.root' OPT='RECREATE'"]
+ # acc.addService(thistSvc)
+
+ return acc
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/SeedingConfig.py b/Tracking/Acts/FaserActsKalmanFilter/python/SeedingConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..e6a84e5b48f6886a0935e6665f4827db55eba958
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/SeedingConfig.py
@@ -0,0 +1,16 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometrySvcCfg
+
+
+def SeedingCfg(flags, **kwargs):
+ acc = FaserSCT_GeometryCfg(flags)
+ acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
+ acc.merge(acts_tracking_geometry_svc )
+ seedingAlg = CompFactory.SeedingAlg(**kwargs)
+ acc.addEventAlgo(seedingAlg)
+ return acc
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/TI12CKF2Config.py b/Tracking/Acts/FaserActsKalmanFilter/python/TI12CKF2Config.py
new file mode 100644
index 0000000000000000000000000000000000000000..748900ed86612abebc50a0f64557223c922d97c8
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/TI12CKF2Config.py
@@ -0,0 +1,125 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS and FASER collaborations
+
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometrySvcCfg
+
+
+def FaserActsAlignmentCondAlgCfg(flags, **kwargs):
+ acc = ComponentAccumulator()
+ acc.addCondAlgo(CompFactory.FaserActsAlignmentCondAlg(name="FaserActsAlignmentCondAlg", **kwargs))
+ return acc
+
+
+def CKF2_OutputCfg(flags):
+ acc = ComponentAccumulator()
+ itemList = ["xAOD::EventInfo#*",
+ "xAOD::EventAuxInfo#*",
+ "TrackCollection#*",
+ ]
+ acc.merge(OutputStreamCfg(flags, "ESD", itemList))
+ ostream = acc.getEventAlgo("OutputStreamESD")
+ ostream.TakeItemsFromInput = True
+ return acc
+
+
+def TI12CKF2Cfg(flags, **kwargs):
+ # acc = ComponentAccumulator()
+ acc = FaserSCT_GeometryCfg(flags)
+ acc.merge(MagneticFieldSvcCfg(flags))
+ # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+ acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
+ acc.merge(acts_tracking_geometry_svc )
+
+ # track_seed_tool = CompFactory.ClusterTrackSeedTool()
+ # track_seed_tool = CompFactory.ActsTrackSeedTool()
+ # track_seed_tool = CompFactory.MyTrackSeedTool()
+ # track_seed_tool = CompFactory.ThreeStationTrackSeedTool()
+ track_seed_tool = CompFactory.CircleFitTrackSeedTool()
+ sigma_loc0 = 1.9e-2
+ sigma_loc1 = 9e-1
+ sigma_phi = 3.3e-2
+ sigma_theta = 2.9e-4
+ p = 1000
+ sigma_p = 0.1 * p
+ sigma_qop = sigma_p / (p * p)
+ # initial_variance_inflation = [1000, 1000, 100, 100, 10000]
+ initial_variance_inflation = [100, 100, 100, 100, 1000]
+ track_seed_tool.covLoc0 = initial_variance_inflation[0] * sigma_loc1 * sigma_loc1
+ track_seed_tool.covLoc1 = initial_variance_inflation[1] * sigma_loc0 * sigma_loc0
+ track_seed_tool.covPhi = initial_variance_inflation[2] * sigma_phi * sigma_phi
+ track_seed_tool.covTheta = initial_variance_inflation[3] * sigma_theta * sigma_theta
+ track_seed_tool.covQOverP = initial_variance_inflation[4] * sigma_qop * sigma_qop
+ track_seed_tool.std_cluster = 0.0231
+ origin = -1900
+ track_seed_tool.origin = origin
+ track_seed_tool.TrackCollection = "Segments"
+
+ trajectory_states_writer_tool = CompFactory.RootTrajectoryStatesWriterTool()
+ trajectory_states_writer_tool.noDiagnostics = kwargs["noDiagnostics"]
+ trajectory_states_writer_tool1 = CompFactory.RootTrajectoryStatesWriterTool()
+ trajectory_states_writer_tool1.noDiagnostics = kwargs["noDiagnostics"]
+ trajectory_states_writer_tool1.FilePath = "track_states_ckf1.root"
+ trajectory_states_writer_tool2 = CompFactory.RootTrajectoryStatesWriterTool()
+ trajectory_states_writer_tool2.FilePath = "track_states_ckf2.root"
+ trajectory_states_writer_tool2.noDiagnostics = kwargs["noDiagnostics"]
+
+ trajectory_summary_writer_tool = CompFactory.RootTrajectorySummaryWriterTool(**kwargs)
+ trajectory_summary_writer_tool.noDiagnostics = kwargs["noDiagnostics"]
+ trajectory_summary_writer_tool1 = CompFactory.RootTrajectorySummaryWriterTool()
+ trajectory_summary_writer_tool1.FilePath = "track_summary_ckf1.root"
+ trajectory_summary_writer_tool1.noDiagnostics = kwargs["noDiagnostics"]
+ trajectory_summary_writer_tool2 = CompFactory.RootTrajectorySummaryWriterTool(**kwargs)
+ trajectory_summary_writer_tool2.FilePath = "track_summary_ckf2.root"
+ trajectory_summary_writer_tool2.noDiagnostics = kwargs["noDiagnostics"]
+
+ actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
+ actsExtrapolationTool.MaxSteps = 1000
+ actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+
+ # trajectory_performance_writer_tool = CompFactory.PerformanceWriterTool("PerformanceWriterTool", **kwargs)
+ # trajectory_performance_writer_tool.ExtrapolationTool = actsExtrapolationTool
+ # trajectory_performance_writer_tool.noDiagnostics = kwargs["noDiagnostics"]
+
+
+ ckf = CompFactory.CKF2(**kwargs)
+ kalman_fitter1 = CompFactory.KalmanFitterTool(name="fitterTool1", **kwargs)
+ kalman_fitter1.noDiagnostics = kwargs["noDiagnostics"]
+ kalman_fitter1.ActsLogging = "INFO"
+ kalman_fitter1.SummaryWriter = True
+ kalman_fitter1.StatesWriter = True
+ kalman_fitter1.SeedCovarianceScale = 10
+ kalman_fitter1.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool1
+ kalman_fitter1.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool1
+ kalman_fitter1.origin = origin
+ ckf.KalmanFitterTool1 = kalman_fitter1
+
+ kalman_fitter2 = CompFactory.KalmanFitterTool(name="fitterTool2", **kwargs)
+ kalman_fitter2.noDiagnostics = kwargs["noDiagnostics"]
+ kalman_fitter2.ActsLogging = "INFO"
+ kalman_fitter2.SummaryWriter = True
+ kalman_fitter2.StatesWriter = True
+ kalman_fitter2.SeedCovarianceScale = 1
+ kalman_fitter2.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool2
+ kalman_fitter2.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool2
+ kalman_fitter2.origin = origin
+ ckf.KalmanFitterTool2 = kalman_fitter2
+
+ ckf.TrackSeed = track_seed_tool
+ ckf.ActsLogging = "INFO"
+ ckf.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool
+ ckf.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool
+ ckf.SummaryWriter = True
+ ckf.StatesWriter = True
+ ckf.PerformanceWriter = False
+ # ckf.PerformanceWriterTool = trajectory_performance_writer_tool
+ ckf.origin = origin
+
+ ckf.nMax = 10
+ ckf.chi2Max = 100000
+ acc.addEventAlgo(ckf)
+ acc.merge(CKF2_OutputCfg(flags))
+ return acc
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index 68aceccd626edab5d1e0316a1364055cb061d160..a3945c24662f12a1c7d800a226b1141c9c0f197b 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -31,7 +31,6 @@
#include "Acts/EventData/Measurement.hpp"
-#include "FaserActsKalmanFilter/CircleFit.h"
size_t CKF2::TrajectoryInfo::nClusters {0};
@@ -98,18 +97,14 @@ StatusCode CKF2::execute() {
m_trackSeedTool->initialTrackParameters();
std::shared_ptr<std::vector<IndexSourceLink>> sourceLinks =
m_trackSeedTool->sourceLinks();
- std::shared_ptr<IdentifierLink> idLinks = m_trackSeedTool->idLinks();
- std::vector<Identifier> ids {};
- for (const auto& idLink : *idLinks) {
- ids.push_back(idLink.second);
- }
+ double origin = m_trackSeedTool->targetZPosition();
+
std::shared_ptr<std::vector<Measurement>> measurements = m_trackSeedTool->measurements();
std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters = m_trackSeedTool->clusters();
std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints = m_trackSeedTool->spacePoints();
std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters = m_trackSeedTool->seedClusters();
TrajectoryInfo::nClusters = sourceLinks->size();
-
TrajectoriesContainer trajectories;
trajectories.reserve(initialParameters->size());
@@ -124,7 +119,7 @@ StatusCode CKF2::execute() {
rotation.col(0) = Acts::Vector3(0, 0, -1);
rotation.col(1) = Acts::Vector3(0, 1, 0);
rotation.col(2) = Acts::Vector3(1, 0, 0);
- Acts::Translation3 trans(0., 0., 0.);
+ Acts::Translation3 trans(0., 0., origin);
Acts::Transform3 trafo(rotation * trans);
initialSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(trafo);
@@ -150,93 +145,7 @@ StatusCode CKF2::execute() {
auto results = (*m_fit)(tmp, *initialParameters, options);
- /*
- for (std::size_t iseed = 0; iseed < results.size(); ++iseed) {
- // The result for this seed
- auto &result = results[iseed];
- if (result.ok()) {
- // Get the track finding output object
- const auto &trackFindingOutput = result.value();
-
- std::unique_ptr<Trk::Track> track = makeTrack(geoctx, result);
-
- if (!trackFindingOutput.fittedParameters.empty()) {
- const std::unordered_map<size_t, Acts::BoundTrackParameters> ¶ms = trackFindingOutput.fittedParameters;
- for (const auto &surface_params: params) {
- ATH_MSG_VERBOSE("Fitted parameters for track " << iseed << ", surface " << surface_params.first);
- ATH_MSG_VERBOSE(" position: " << surface_params.second.position(geoctx).transpose());
- ATH_MSG_VERBOSE(" momentum: " << surface_params.second.momentum().transpose());
- // const auto& [currentTip, tipState] = trackFindingOutput.activeTips.back();
- // ATH_MSG_VERBOSE(" #measurements: " << tipState.nMeasurements);
- }
- } else {
- ATH_MSG_DEBUG("No fitted parameters for track " << iseed);
- }
-
- m_kalmanFitterTool1->fit(ctx, *track);
-
- // Create a Trajectories result struct
- trajectories.emplace_back(trackFindingOutput.fittedStates,
- trackFindingOutput.lastMeasurementIndices,
- trackFindingOutput.fittedParameters);
- } else {
- ATH_MSG_WARNING("Track finding failed for seed " << iseed << " with error" << result.error());
- // Track finding failed. Add an empty result so the output container has
- // the same number of entries as the input.
- trajectories.push_back(FaserActsRecMultiTrajectory());
- }
- }
- */
-
-
- // std::vector<TrackQuality> trackQuality;
- // selectTracks(results, trackQuality);
- //
- // TrackFinderResult selectedResults;
- // for (size_t i = 0; i < std::min(trackQuality.size(), (size_t)m_nTrajectories); ++i) {
- // selectedResults.push_back(Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(trackQuality[i].track));
- // }
-
- /*
- TrackFinderResult minTrackRequirements {};
- for (auto& result : results) {
- if (not result.ok()) {
- continue;
- }
- auto& ckfResult = result.value();
- auto traj = FaserActsRecMultiTrajectory(ckfResult.fittedStates, ckfResult.lastMeasurementIndices, ckfResult.fittedParameters);
- const auto& mj = traj.multiTrajectory();
- const auto& trackTips = traj.tips();
- size_t maxMeasurements = 0;
- for (const auto& trackTip : trackTips) {
- auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
- size_t nMeasurements = trajState.nMeasurements;
- if (nMeasurements > maxMeasurements) {
- maxMeasurements = nMeasurements;
- }
- }
- if (maxMeasurements >= m_minNumberMeasurements) {
- minTrackRequirements.push_back(Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(ckfResult));
- }
- }
-
- TrackFinderResult selectedResults {};
- if (minTrackRequirements.size() > 2) {
- std::pair<std::size_t, std::size_t> trackIndices = solveAmbiguity(results);
- selectedResults.push_back(
- Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(results.at(trackIndices.first).value()));
- selectedResults.push_back(
- Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(results.at(trackIndices.second).value()));
- } else {
- for (auto& result : minTrackRequirements) {
- selectedResults.push_back(
- Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(result.value()));
- }
- }
- */
- // computeSharedHits(sourceLinks.get(), selectedResults);
-
- // results contains a MultiTrajectory for each track seed with a trajectory of each branch of the CKF.
+ // results contain a MultiTrajectory for each track seed with a trajectory of each branch of the CKF.
// To simplify the ambiguity solving a list of MultiTrajectories is created, each containing only a single track.
std::list<TrajectoryInfo> allTrajectories;
for (auto &result : results) {
@@ -281,9 +190,9 @@ StatusCode CKF2::execute() {
std::unique_ptr<Trk::Track> track = makeTrack(geoctx, traj);
if (track) {
// outputTracks->push_back(std::move(track));
- std::unique_ptr<Trk::Track> track2 = m_kalmanFitterTool1->fit(ctx, *track, trajectories, Acts::BoundVector::Zero());
+ std::unique_ptr<Trk::Track> track2 = m_kalmanFitterTool1->fit(ctx, *track, trajectories, Acts::BoundVector::Zero(), m_isMC, origin);
if (track2) {
- std::unique_ptr<Trk::Track> track3 = m_kalmanFitterTool2->fit(ctx, *track2, trajectories, Acts::BoundVector::Zero());
+ std::unique_ptr<Trk::Track> track3 = m_kalmanFitterTool2->fit(ctx, *track2, trajectories, Acts::BoundVector::Zero(), m_isMC, origin);
outputTracks->push_back(std::move(track3));
}
}
@@ -291,10 +200,10 @@ StatusCode CKF2::execute() {
// run the performance writer
if (m_statesWriter && !m_noDiagnostics) {
- ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, selectedTrajectories));
+ ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, selectedTrajectories, m_isMC));
}
if (m_summaryWriter && !m_noDiagnostics) {
- ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, selectedTrajectories));
+ ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, selectedTrajectories, m_isMC));
}
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool->write(geoctx, selectedTrajectories));
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFit.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFit.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9f1fab11cf10edac1a497c7af498e737b5bd84c6
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFit.cxx
@@ -0,0 +1,169 @@
+#include "FaserActsKalmanFilter/CircleFit.h"
+
+namespace CircleFit {
+
+CircleData::CircleData(const std::vector<const Tracker::FaserSCT_SpacePoint *> &spacePoints) {
+ for (auto sp: spacePoints) {
+ m_x.push_back(sp->globalPosition().z());
+ m_y.push_back(sp->globalPosition().y());
+ }
+ m_size = spacePoints.size();
+}
+
+CircleData::CircleData(const std::vector<Amg::Vector3D> &spacePoints) {
+ for (auto sp: spacePoints) {
+ m_x.push_back(sp.z());
+ m_y.push_back(sp.y());
+ }
+ m_size = spacePoints.size();
+}
+
+double CircleData::meanX() const {
+ double mean = 0;
+ for (double i: m_x) mean += i;
+ return mean / m_size;
+}
+
+double CircleData::meanY() const {
+ double mean = 0;
+ for (double i: m_y) mean += i;
+ return mean / m_size;
+}
+
+double CircleData::x(int i) const {
+ return m_x[i];
+}
+
+double CircleData::y(int i) const {
+ return m_y[i];
+}
+
+int CircleData::size() const {
+ return m_size;
+}
+
+
+double sigma(const CircleData &data, const Circle &circle) {
+ double sum=0.,dx,dy;
+ for (int i=0; i<data.size(); i++) {
+ dx = data.x(i) - circle.cx;
+ dy = data.y(i) - circle.cy;
+ sum += (sqrt(dx*dx+dy*dy) - circle.r) * (sqrt(dx*dx+dy*dy) - circle.r);
+ }
+ return sqrt(sum/data.size());
+}
+
+
+Circle circleFit(const CircleData &data)
+/*
+ Circle fit to a given set of data points (in 2D)
+
+ This is an algebraic fit, due to Taubin, based on the journal article
+
+ G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
+ Space Curves Defined By Implicit Equations, With
+ Applications To Edge And Range Image Segmentation",
+ IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
+
+ Input: data - the class of data (contains the given points):
+
+ data.size() - the number of data points
+ data.X[] - the array of X-coordinates
+ data.Y[] - the array of Y-coordinates
+
+ Output:
+ circle - parameters of the fitting circle:
+
+ circle.a - the X-coordinate of the center of the fitting circle
+ circle.b - the Y-coordinate of the center of the fitting circle
+ circle.r - the radius of the fitting circle
+ circle.s - the root mean square error (the estimate of sigma)
+ circle.j - the total number of iterations
+
+ The method is based on the minimization of the function
+
+ sum [(x-a)^2 + (y-b)^2 - R^2]^2
+ F = -------------------------------
+ sum [(x-a)^2 + (y-b)^2]
+
+ This method is more balanced than the simple Kasa fit.
+
+ It works well whether data points are sampled along an entire circle or
+ along a small arc.
+
+ It still has a small bias and its statistical accuracy is slightly
+ lower than that of the geometric fit (minimizing geometric distances),
+ but slightly higher than that of the very similar Pratt fit.
+ Besides, the Taubin fit is slightly simpler than the Pratt fit
+
+ It provides a very good initial guess for a subsequent geometric fit.
+
+ Nikolai Chernov (September 2012)
+
+*/
+{
+ int i, iter, IterMAX = 999;
+
+ double Xi, Yi, Zi;
+ double Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Cov_xy, Var_z;
+ double A0, A1, A2, A22, A3, A33;
+ double Dy, xnew, x, ynew, y;
+ double DET, Xcenter, Ycenter;
+
+ Circle circle;
+
+ Mxx = Myy = Mxy = Mxz = Myz = Mzz = 0.;
+ for (i = 0; i < data.size(); i++) {
+ Xi = data.x(i) - data.meanX();
+ Yi = data.y(i) - data.meanY();
+ Zi = Xi * Xi + Yi * Yi;
+
+ Mxy += Xi * Yi;
+ Mxx += Xi * Xi;
+ Myy += Yi * Yi;
+ Mxz += Xi * Zi;
+ Myz += Yi * Zi;
+ Mzz += Zi * Zi;
+ }
+ Mxx /= data.size();
+ Myy /= data.size();
+ Mxy /= data.size();
+ Mxz /= data.size();
+ Myz /= data.size();
+ Mzz /= data.size();
+
+ Mz = Mxx + Myy;
+ Cov_xy = Mxx * Myy - Mxy * Mxy;
+ Var_z = Mzz - Mz * Mz;
+ A3 = 4 * Mz;
+ A2 = -3 * Mz * Mz - Mzz;
+ A1 = Var_z * Mz + 4 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz;
+ A0 = Mxz * (Mxz * Myy - Myz * Mxy) + Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy;
+ A22 = A2 + A2;
+ A33 = A3 + A3 + A3;
+
+ for (x = 0., y = A0, iter = 0; iter < IterMAX; iter++) {
+ Dy = A1 + x * (A22 + A33 * x);
+ xnew = x - y / Dy;
+ if ((xnew == x) || (!std::isfinite(xnew))) break;
+ ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
+ if (abs(ynew) >= abs(y)) break;
+ x = xnew;
+ y = ynew;
+ }
+
+ DET = x * x - x * Mz + Cov_xy;
+ Xcenter = (Mxz * (Myy - x) - Myz * Mxy) / DET / 2;
+ Ycenter = (Myz * (Mxx - x) - Mxz * Mxy) / DET / 2;
+
+ circle.cx = Xcenter + data.meanX();
+ circle.cy = Ycenter + data.meanY();
+ circle.r = sqrt(Xcenter * Xcenter + Ycenter * Ycenter + Mz);
+ circle.s = sigma(data, circle);
+ circle.i = 0;
+ circle.j = iter;
+
+ return circle;
+}
+
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..18499e4964e49ebf4ef041a870125bdae2d2d694
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
@@ -0,0 +1,346 @@
+#include "FaserActsKalmanFilter/CircleFitTrackSeedTool.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterCollection.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "FaserActsKalmanFilter/CircleFit.h"
+#include "FaserActsKalmanFilter/LinearFit.h"
+#include "FaserActsKalmanFilter/TrackClassification.h"
+#include <array>
+#include <algorithm>
+
+
+std::map<Identifier,Index> CircleFitTrackSeedTool::s_indexMap {};
+std::map<Identifier, const Tracker::FaserSCT_SpacePoint*> CircleFitTrackSeedTool::s_spacePointMap {};
+
+CircleFitTrackSeedTool::CircleFitTrackSeedTool(
+ const std::string& type, const std::string& name, const IInterface* parent)
+ : base_class(type, name, parent) {}
+
+
+StatusCode CircleFitTrackSeedTool::initialize() {
+ ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+ ATH_CHECK(detStore()->retrieve(m_detManager, "SCT"));
+ ATH_CHECK(m_trackingGeometryTool.retrieve());
+ ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_clusterContainerKey.initialize());
+ ATH_CHECK(m_spacePointContainerKey.initialize());
+ // ATH_CHECK(m_simDataCollectionKey.initialize());
+ // ATH_CHECK(m_mcEventCollectionKey.initialize());
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode CircleFitTrackSeedTool::run() {
+ SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer {m_clusterContainerKey};
+ ATH_CHECK(clusterContainer.isValid());
+
+ SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey};
+ ATH_CHECK(spacePointContainer.isValid());
+
+ // SG::ReadHandle<TrackerSimDataCollection> simData {m_simDataCollectionKey};
+ // ATH_CHECK(simData.isValid());
+
+ // SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey};
+ // ATH_CHECK(mcEvents.isValid());
+
+ using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
+ std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
+ // std::map<int, const HepMC::GenParticle*> particles {};
+ // for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
+ // particles[particle->barcode()] = particle;
+ // }
+
+ // std::map<Identifier, const Tracker::FaserSCT_SpacePoint*> s_spacePointMap {};
+ std::vector<const Tracker::FaserSCT_SpacePoint*> spacePoints {};
+ for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) {
+ for (const Tracker::FaserSCT_SpacePoint *spacePoint: *spacePointCollection) {
+ spacePoints.push_back(spacePoint);
+ Identifier id1 = spacePoint->cluster1()->identify();
+ CircleFitTrackSeedTool::s_spacePointMap[id1] = spacePoint;
+ Identifier id2 = spacePoint->cluster2()->identify();
+ CircleFitTrackSeedTool::s_spacePointMap[id2] = spacePoint;
+ }
+ }
+
+
+ using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, 1>;
+ std::array<Acts::BoundIndices, 1> Indices = {Acts::eBoundLoc0};
+ std::vector<IndexSourceLink> sourceLinks {};
+ std::vector<Measurement> measurements {};
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters {};
+ for (const Tracker::FaserSCT_ClusterCollection* clusterCollection : *clusterContainer) {
+ for (const Tracker::FaserSCT_Cluster* cluster : *clusterCollection) {
+ clusters.push_back(cluster);
+ Identifier id = cluster->detectorElement()->identify();
+ CircleFitTrackSeedTool::s_indexMap[cluster->identify()] = measurements.size();
+ if (identifierMap->count(id) != 0) {
+ Acts::GeometryIdentifier geoId = identifierMap->at(id);
+ IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
+ // create measurement
+ const auto& par = cluster->localPosition();
+ Eigen::Matrix<double, 1, 1> pos {par.x(),};
+ Eigen::Matrix<double, 1, 1> cov {0.0231 * 0.0231,};
+ ThisMeasurement meas(sourceLink, Indices, pos, cov);
+ sourceLinks.push_back(sourceLink);
+ measurements.emplace_back(std::move(meas));
+ }
+ }
+ }
+
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection};
+ ATH_CHECK(trackCollection.isValid());
+
+ std::array<std::vector<Segment>, 4> segments {};
+ for (const Trk::Track* track : *trackCollection) {
+ auto s = Segment(track, m_idHelper);
+ segments[s.station].push_back(s);
+ }
+
+ std::vector<Segment> combination {};
+ std::vector<Seed> seeds {};
+ // create seeds from four stations
+ go(segments, combination, seeds, 0, 4);
+ if (seeds.size() < 2) {
+ // create seeds from three stations
+ go(segments, combination, seeds, 0, 3);
+ }
+ // create seeds from two stations
+ if (seeds.size() < 2) {
+ go(segments, combination, seeds, 0, 2);
+ }
+ if (seeds.size() < 2) {
+ go(segments, combination, seeds, 0, 1);
+ }
+
+ std::list<Seed> allSeeds;
+ for (const Seed &seed : seeds) allSeeds.push_back(seed);
+
+ /*
+ allSeeds.sort([](const Seed &left, const Seed &right) {
+ if (left.size > right.size) return true;
+ if (left.size < right.size) return false;
+ if (left.chi2 < right.chi2) return true;
+ else return false;
+ });
+
+ std::vector<Seed> selectedSeeds {};
+ while (not allSeeds.empty()) {
+ Seed selected = allSeeds.front();
+ selectedSeeds.push_back(selected);
+ allSeeds.remove_if([&](const Seed &p) {
+ return (p.size <= 12) || ((p.clusterSet & selected.clusterSet).count() > 6);
+ });
+ }
+ */
+
+ std::vector<Seed> selectedSeeds {};
+ for (const Seed &seed : allSeeds) {
+ selectedSeeds.push_back(seed);
+ }
+
+ Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
+ cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = m_covLoc0;
+ cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = m_covLoc1;
+ cov(Acts::eBoundPhi, Acts::eBoundPhi) = m_covPhi;
+ cov(Acts::eBoundTheta, Acts::eBoundTheta) = m_covTheta;
+ cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = m_covQOverP;
+ cov(Acts::eBoundTime, Acts::eBoundTime) = m_covTime;
+
+ auto minSeed = std::min_element(
+ selectedSeeds.begin(), selectedSeeds.end(), [](const Seed &lhs, const Seed &rhs) {
+ return lhs.minZ < rhs.minZ;
+ });
+ double origin = !selectedSeeds.empty() ? minSeed->minZ : 0;
+ m_targetZPosition = origin;
+ std::vector<Acts::CurvilinearTrackParameters> initParams {};
+ for (const Seed &seed : selectedSeeds) {
+ initParams.push_back(seed.get_params(origin, cov));
+ }
+
+ m_initialTrackParameters = std::make_shared<std::vector<Acts::CurvilinearTrackParameters>>(initParams);
+ m_sourceLinks = std::make_shared<std::vector<IndexSourceLink>>(sourceLinks);
+ // m_idLinks = std::make_shared<IdentifierLink>(identifierLinkMap);
+ m_measurements = std::make_shared<std::vector<Measurement>>(measurements);
+ m_initialSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ Acts::Vector3 {0, 0, origin}, Acts::Vector3{0, 0, -1});
+ m_clusters = std::make_shared<std::vector<const Tracker::FaserSCT_Cluster*>>(clusters);
+ m_spacePoints = std::make_shared<std::vector<const Tracker::FaserSCT_SpacePoint*>>(spacePoints);
+
+ /*
+ std::cout.precision(17);
+ for (auto &seed : selectedSeeds) {
+ std::cout << "np.array([";
+ for (const Acts::Vector3 &pos : seed.fakePositions) {
+ std::cout << "[" << pos.x() << ", " << pos.y() << ", " << pos.z() << "], ";
+ }
+ std::cout << "])" << std::endl;
+ std::cout << "chi2: " << seed.chi2 << ", momentum: " << seed.momentum << ", charge: " << seed.charge << std::endl;
+ std::cout << "fit = np.array([" << seed.c1 << ", " << seed.c0 << ", " << seed.cy << ", " << seed.cx << ", " << seed.r << "])" << std::endl;
+
+ std::vector<ParticleHitCount> particleHitCounts;
+ identifyContributingParticles(*simData, seed.clusters, particleHitCounts);
+ auto ip = particles.find(particleHitCounts.front().particleId);
+ if (ip != particles.end()) {
+ const HepMC::GenParticle* particle = ip->second;
+ HepMC::FourVector momentum = particle->momentum();
+ std::cout << "true momentum: " << momentum.rho() * 0.001 << std::endl;
+ }
+ for (const ParticleHitCount &hitCount : particleHitCounts) {
+ std::cout << hitCount.particleId << " : " << hitCount.hitCount << std::endl;
+ }
+ }
+ */
+
+ s_indexMap.clear();
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode CircleFitTrackSeedTool::finalize() {
+ return StatusCode::SUCCESS;
+}
+
+
+void CircleFitTrackSeedTool::go(const std::array<std::vector<Segment>, 4> &v,
+ std::vector<Segment> &combination,
+ std::vector<Seed> &seeds,
+ int offset, int k) {
+ if (k == 0) {
+ seeds.push_back(Seed(combination));
+ return;
+ }
+ for (std::size_t i = offset; i < v.size() + 1 - k; ++i) {
+ for (const auto& ve : v[i]) {
+ combination.push_back(ve);
+ go(v, combination, seeds, i+1, k-1);
+ combination.pop_back();
+ }
+ }
+}
+
+CircleFitTrackSeedTool::Segment::Segment(const Trk::Track* track, const FaserSCT_ID *idHelper) :
+ clusterSet(CircleFitTrackSeedTool::s_indexMap.size()) {
+ for (const Trk::TrackStateOnSurface* trackState : *(track->trackStateOnSurfaces())) {
+ auto clusterOnTrack = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*> (trackState->measurementOnTrack());
+ if (clusterOnTrack) {
+ const Tracker::FaserSCT_Cluster* cluster = clusterOnTrack->prepRawData();
+ Identifier id = cluster->identify();
+ clusters.push_back(cluster);
+ if (CircleFitTrackSeedTool::s_spacePointMap.count(id) > 0) {
+ const Tracker::FaserSCT_SpacePoint *sp = CircleFitTrackSeedTool::s_spacePointMap.at(cluster->identify());
+ if (std::find(spacePoints.begin(), spacePoints.end(), sp) == spacePoints.end()) {
+ spacePoints.push_back(sp);
+ }
+ }
+ station = idHelper->station(id);
+ clusterSet.set(CircleFitTrackSeedTool::s_indexMap.at(id));
+ auto fitParameters = track->trackParameters()->front();
+ position = fitParameters->position();
+ momentum = fitParameters->momentum();
+ }
+ }
+ fakePositions.push_back(position);
+ fakePositions.push_back(position - 30 * momentum.normalized());
+ fakePositions.push_back(position + 30 * momentum.normalized());
+}
+
+CircleFitTrackSeedTool::Seed::Seed(const std::vector<Segment> &segments) :
+ clusterSet(CircleFitTrackSeedTool::s_indexMap.size()) {
+ for (const Segment &seg : segments) {
+ clusters.insert(clusters.end(), seg.clusters.begin(), seg.clusters.end());
+ spacePoints.insert(spacePoints.end(), seg.spacePoints.begin(), seg.spacePoints.end());
+ positions.push_back(seg.position);
+ // TODO use reconstruct space points instead of fake positions
+ fakePositions.insert(fakePositions.end(), seg.fakePositions.begin(), seg.fakePositions.end());
+ for (size_t i = 0; i < seg.clusterSet.size(); ++i) {
+ if (seg.clusterSet[i]) clusterSet.set(i);
+ }
+ }
+
+ auto minCluster = std::min_element(
+ clusters.begin(), clusters.end(), [](const Tracker::FaserSCT_Cluster *lhs, const Tracker::FaserSCT_Cluster *rhs) {
+ return lhs->globalPosition().z() < rhs->globalPosition().z();
+ } );
+ minZ = (*minCluster)->globalPosition().z();
+
+ if (segments.size() > 1) {
+ direction = positions[1] - positions[0];
+ } else {
+ direction = segments[0].momentum;
+ }
+
+ std::vector<Acts::Vector2> pointsZX {};
+ for (const Acts::Vector3 &pos : fakePositions) {
+ pointsZX.push_back({pos.z(), pos.x()});
+ }
+ linearFit(pointsZX);
+
+ if (segments.size() > 1) {
+ fakeFit();
+ } else {
+ momentum = 9999999.;
+ charge = 1;
+ }
+
+ getChi2();
+ size = clusters.size();
+}
+
+
+void CircleFitTrackSeedTool::Seed::fit() {
+ CircleFit::CircleData circleData(positions);
+ CircleFit::Circle circle = CircleFit::circleFit(circleData);
+ momentum = circle.r > 0 ? circle.r * 0.001 * 0.3 * 0.55 : 9999999.;
+ charge = circle.cy < 0 ? -1 : 1;
+}
+
+void CircleFitTrackSeedTool::Seed::fakeFit(double B) {
+ CircleFit::CircleData circleData(fakePositions);
+ CircleFit::Circle circle = CircleFit::circleFit(circleData);
+ cx = circle.cx;
+ cy = circle.cy;
+ r = circle.r;
+ momentum = r * 0.3 * B * m_MeV2GeV;
+ charge = circle.cy > 0 ? 1 : -1;
+}
+
+void CircleFitTrackSeedTool::Seed::linearFit(const std::vector<Acts::Vector2> &points) {
+ auto [origin, dir] = LinearFit::linearFit(points);
+ c1 = dir[1]/dir[0];
+ c0 = origin[1] - origin[0] * c1;
+}
+
+
+double CircleFitTrackSeedTool::Seed::getY(double z) {
+ double sqt = std::sqrt(-cx*cx + 2*cx*z + r*r - z*z);
+ return abs(cy - sqt) < abs(cy + sqt) ? cy - sqt : cy + sqt;
+}
+
+
+double CircleFitTrackSeedTool::Seed::getX(double z) {
+ return c0 + z * c1;
+}
+
+void CircleFitTrackSeedTool::Seed::getChi2() {
+ chi2 = 0;
+ for (const Acts::Vector3 &pos : fakePositions) {
+ m_dy = pos.y() - getY(pos.z());
+ chi2 += (m_dy * m_dy) / (m_sigma_y * m_sigma_y);
+ }
+
+ for (const Acts::Vector3 &pos : positions) {
+ m_dx = pos.x() - getX(pos.x());
+ chi2 += (m_dx * m_dx) / (m_sigma_x * m_sigma_x);
+ }
+}
+
+Acts::CurvilinearTrackParameters CircleFitTrackSeedTool::Seed::get_params(double origin, Acts::BoundSymMatrix cov) const {
+ Acts::Vector3 pos = positions[0] - (positions[0].z() - origin)/direction.z() * direction;
+ Acts::Vector4 pos4 {pos.x(), pos.y(), pos.z(), 0};
+ return Acts::CurvilinearTrackParameters(pos4, direction.normalized(), momentum, charge, cov);
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
index 2ee66e278cd3488a4c99611c11d73a21bac2b3d2..71afcc4431fabb6064af9601d31dd783de497a4f 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
@@ -171,10 +171,10 @@ StatusCode CombinatorialKalmanFilterAlg::execute() {
// run the performance writer
if (m_statesWriter && !m_noDiagnostics) {
- ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, trajectories));
+ ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, trajectories, m_isMC));
}
if (m_summaryWriter && !m_noDiagnostics) {
- ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, trajectories));
+ ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, trajectories, m_isMC));
}
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool->write(geoctx, trajectories));
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6cfb07d143b8b49652ad54ee768a1af2b0671fae
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
@@ -0,0 +1,114 @@
+#include "FaserActsKalmanFilter/GhostBusters.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+
+
+GhostBusters::GhostBusters(const std::string &name, ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator), AthHistogramming(name),
+ m_histSvc("THistSvc/THistSvc", name) {}
+
+
+StatusCode GhostBusters::initialize() {
+ ATH_CHECK(m_trackCollection.initialize());
+ // ATH_CHECK(m_simDataCollectionKey.initialize());
+ ATH_CHECK(m_outputTrackCollection.initialize());
+ ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+
+ m_tree = new TTree("tree", "tree");
+ m_tree->Branch("event_number", &m_event_number, "event_number/I");
+ m_tree->Branch("station", &m_station, "station/I");
+ m_tree->Branch("x", &m_x, "x/D");
+ m_tree->Branch("y", &m_y, "y/D");
+ m_tree->Branch("z", &m_z, "z/D");
+ m_tree->Branch("chi2", &m_chi2, "chi2/D");
+ // m_tree->Branch("majorityHits", &m_majorityHits, "majorityHits/I");
+ m_tree->Branch("size", &m_size, "size/I");
+ m_tree->Branch("isGhost", &m_isGhost, "isGhost/B");
+ ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode GhostBusters::execute(const EventContext &ctx) const {
+ m_event_number = ctx.eventID().event_number();
+
+ SG::WriteHandle outputTrackCollection {m_outputTrackCollection, ctx};
+ std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+
+ // SG::ReadHandle<TrackerSimDataCollection> simData {m_simDataCollectionKey, ctx};
+ // ATH_CHECK(simData.isValid());
+
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+ ATH_CHECK(trackCollection.isValid());
+
+ std::map<int, std::vector<Segment>> segments {};
+ for (const Trk::Track* track : *trackCollection) {
+ auto segment = Segment(track, m_idHelper);
+ segments[segment.station()].push_back(segment);
+ }
+
+ for (const auto &station : segments) {
+ double nGoodSegments = 0;
+ std::vector<Segment> stationSegments = station.second;
+ std::sort(stationSegments.begin(), stationSegments.end(), [](const Segment &lhs, const Segment &rhs) { return lhs.y() > rhs.y(); });
+ double maxX = stationSegments.front().x();
+ double maxY = stationSegments.front().y();
+ double minX = stationSegments.back().x();
+ double minY = stationSegments.back().y();
+ double deltaY = maxY - minY;
+ double meanX = 0.5 * (minX + maxX);
+ double meanY = 0.5 * (minY + maxY);
+ double deltaX = deltaY / (2 * std::tan(0.02));
+ for (Segment &segment : stationSegments) {
+ bool isGhost = (segment.y() > meanY - m_yTolerance * deltaY) && (segment.y() < meanY + m_yTolerance * deltaY) && (
+ ((segment.x() > meanX - (1 + m_xTolerance) * deltaX) &&
+ (segment.x() < meanX - (1 - m_xTolerance) * deltaX)) ||
+ ((segment.x() > meanX + (1 - m_xTolerance) * deltaX) && (segment.x() < meanX + (1 + m_xTolerance) * deltaX)));
+ if (isGhost) segment.setGhost();
+ if (not isGhost && segment.size() >= 5) nGoodSegments++;
+ }
+ for (const Segment &segment : stationSegments) {
+ m_x = segment.x();
+ m_y = segment.y();
+ m_z = segment.z();
+ m_chi2 = segment.chi2();
+ m_station = segment.station();
+ m_size = segment.size();
+ // m_majorityHits = segment.majorityHits();
+ m_isGhost = segment.isGhost();
+ if (nGoodSegments >= 2 && segment.size() == 4) m_isGhost = true;
+ m_tree->Fill();
+ if (segment.isGhost()) continue;
+ if (nGoodSegments >= 2 && segment.size() == 4) continue;
+ outputTracks->push_back(new Trk::Track(*segment.track()));
+ }
+ }
+
+ ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode GhostBusters::finalize() {
+ return StatusCode::SUCCESS;
+}
+
+
+GhostBusters::Segment::Segment(const Trk::Track *track, const FaserSCT_ID *idHelper) {
+ m_track = track;
+ m_position = track->trackParameters()->front()->position();
+ m_chi2 = track->fitQuality()->chiSquared();
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters {};
+ for (const Trk::TrackStateOnSurface* trackState : *(track->trackStateOnSurfaces())) {
+ auto clusterOnTrack = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*> (trackState->measurementOnTrack());
+ if (clusterOnTrack) {
+ clusters.emplace_back(clusterOnTrack->prepRawData());
+ m_station = idHelper->station(clusterOnTrack->identify());
+ }
+ }
+ m_size = clusters.size();
+ // identifyContributingParticles(simDataCollection, clusters, m_particleHitCounts);
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
index 7e97a82dfc1bc93e95714c9bf00fcdbc75d405c7..eddf4cabd33698632f487d1f825388e514416fc6 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
@@ -41,7 +41,7 @@ StatusCode KalmanFitterTool::finalize() {
std::unique_ptr<Trk::Track>
KalmanFitterTool::fit(const EventContext &ctx, const Trk::Track &inputTrack,
std::vector<FaserActsRecMultiTrajectory> & /*trajectories*/,
- const Acts::BoundVector& inputVector = Acts::BoundVector::Zero()) const {
+ const Acts::BoundVector& inputVector = Acts::BoundVector::Zero(), bool isMC=false, double origin=0) const {
std::unique_ptr<Trk::Track> newTrack = nullptr;
std::vector<FaserActsRecMultiTrajectory> myTrajectories;
@@ -57,16 +57,16 @@ KalmanFitterTool::fit(const EventContext &ctx, const Trk::Track &inputTrack,
auto pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
- Acts::Vector3 {0, 0, 0}, Acts::Vector3{0, 0, -1});
+ Acts::Vector3 {0, 0, origin}, Acts::Vector3{0, 0, -1});
Acts::GeometryContext gctx = m_trackingGeometryTool->getNominalGeometryContext().context();
Acts::MagneticFieldContext mfContext = getMagneticFieldContext(ctx);
Acts::CalibrationContext calibContext = Acts::CalibrationContext();
auto [sourceLinks, measurements] = getMeasurementsFromTrack(inputTrack);
- auto trackParameters = getParametersFromTrack(inputTrack.trackParameters()->front(), inputVector);
+ auto trackParameters = getParametersFromTrack(inputTrack.trackParameters()->front(), inputVector, origin);
ATH_MSG_DEBUG("trackParameters: " << trackParameters.parameters().transpose());
- ATH_MSG_DEBUG("position: " << trackParameters.parameters().transpose());
+ ATH_MSG_DEBUG("position: " << trackParameters.position(gctx).transpose());
ATH_MSG_DEBUG("momentum: " << trackParameters.momentum().transpose());
ATH_MSG_DEBUG("charge: " << trackParameters.charge());
@@ -108,10 +108,10 @@ KalmanFitterTool::fit(const EventContext &ctx, const Trk::Track &inputTrack,
}
if (m_statesWriter && !m_noDiagnostics) {
- StatusCode statusStatesWriterTool = m_trajectoryStatesWriterTool->write(gctx, myTrajectories);
+ StatusCode statusStatesWriterTool = m_trajectoryStatesWriterTool->write(gctx, myTrajectories, isMC);
}
if (m_summaryWriter && !m_noDiagnostics) {
- StatusCode statusSummaryWriterTool = m_trajectorySummaryWriterTool->write(gctx, myTrajectories);
+ StatusCode statusSummaryWriterTool = m_trajectorySummaryWriterTool->write(gctx, myTrajectories, isMC);
}
return newTrack;
@@ -205,10 +205,10 @@ KalmanFitterTool::getMeasurementsFromTrack(const Trk::Track &track) const {
Acts::BoundTrackParameters
KalmanFitterTool::getParametersFromTrack(const Trk::TrackParameters *inputParameters,
- const Acts::BoundVector& inputVector) const {
+ const Acts::BoundVector& inputVector, double origin) const {
using namespace Acts::UnitLiterals;
std::shared_ptr<const Acts::Surface> pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
- Acts::Vector3 {0, 0, 0}, Acts::Vector3{0, 0, -1});
+ Acts::Vector3 {0, 0, origin}, Acts::Vector3{0, 0, -1});
auto atlasParam = inputParameters->parameters();
auto center = inputParameters->associatedSurface().center();
@@ -286,7 +286,7 @@ KalmanFitterTool::makeTrack(Acts::GeometryContext& geoCtx, TrackFitterResult& fi
state.smoothed(), state.smoothedCovariance());
actsSmoothedParam.push_back(std::make_unique<const Acts::BoundTrackParameters>(Acts::BoundTrackParameters(actsParam)));
// const auto& psurface=actsParam.referenceSurface();
- Acts::Vector2 local(actsParam.parameters()[Acts::eBoundLoc0], actsParam.parameters()[Acts::eBoundLoc1]);
+ // Acts::Vector2 local(actsParam.parameters()[Acts::eBoundLoc0], actsParam.parameters()[Acts::eBoundLoc1]);
// const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(actsParam.parameters()[Acts::eBoundPhi], actsParam.parameters()[Acts::eBoundTheta]);
// auto pos=actsParam.position(tgContext);
parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx);
@@ -310,7 +310,7 @@ KalmanFitterTool::makeTrack(Acts::GeometryContext& geoCtx, TrackFitterResult& fi
double nDoF = state.calibratedSize();
const Trk::FitQualityOnSurface *quality = new Trk::FitQualityOnSurface(state.chi2(), nDoF);
const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(measState, parm, quality, nullptr, typePattern);
- // If a state was succesfully created add it to the trajectory
+ // If a state was successfully created add it to the trajectory
if (perState) {
finalTrajectory->insert(finalTrajectory->begin(), perState);
}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.cxx
index 11c1483a7762679e42ce41542d1a5b165aaf4f0c..47922a13e99abca49446668f9f8c180a86ebecdf 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.cxx
@@ -119,7 +119,7 @@ void ResPlotTool::clear(ResPlotCache &resPlotCache) const {
}
void ResPlotTool::fill(
- ResPlotCache& resPlotCache, const Acts::GeometryContext& gctx,
+ ResPlotCache& resPlotCache, const Acts::GeometryContext& /*gctx*/,
std::unique_ptr<const Acts::BoundTrackParameters> truthParameters,
const Acts::BoundTrackParameters& fittedParameters) const {
using ParametersVector = Acts::BoundTrackParameters::ParametersVector;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
index 634e356a887d45afa6ddf7c06a1d8e4566d0ef8b..5378fc080cb49628248d4a7eda2e9206d6c016ba 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
@@ -13,6 +13,9 @@
#include <TFile.h>
#include <TTree.h>
+constexpr float NaNfloat = std::numeric_limits<float>::quiet_NaN();
+constexpr float NaNint = std::numeric_limits<int>::quiet_NaN();
+
using Acts::VectorHelpers::eta;
using Acts::VectorHelpers::perp;
using Acts::VectorHelpers::phi;
@@ -194,6 +197,7 @@ StatusCode RootTrajectoryStatesWriterTool::initialize() {
m_outputTree->Branch("chi2", &m_chi2);
}
+
return StatusCode::SUCCESS;
}
@@ -207,51 +211,46 @@ StatusCode RootTrajectoryStatesWriterTool::finalize() {
return StatusCode::SUCCESS;
}
-StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories) const {
+StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories, bool isMC) const {
if (m_outputFile == nullptr)
return StatusCode::SUCCESS;
+ // Get the event number
const EventContext& ctx = Gaudi::Hive::currentContext();
+ m_eventNr = ctx.eventID().event_number();
- 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;
- }
- ATH_MSG_VERBOSE("Found " << mcEvents->front()->particles_size() << " particles.");
+ std::vector<ParticleHitCount> particleHitCounts;
std::map<int, const HepMC::GenParticle*> particles {};
- for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
- particles[particle->barcode()] = particle;
- }
-
- SG::ReadHandle<TrackerSimDataCollection> simData {m_simDataCollectionKey, ctx};
- ATH_CHECK(simData.isValid());
-
- SG::ReadHandle<FaserSiHitCollection> siHitCollection {m_faserSiHitKey, ctx};
- ATH_CHECK(siHitCollection.isValid());
std::map<std::pair<int, Identifier>, const FaserSiHit*> siHitMap;
- for (const FaserSiHit& hit : *siHitCollection) {
- int barcode = hit.trackNumber();
- Identifier id = m_idHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(), hit.getModule(), hit.getSensor());
- siHitMap[std::make_pair(barcode, id)] = &hit;
- }
- // using HitParticlesMap = IndexMultimap<ActsFatras::Barcode>;
- // using HitSimHitsMap = IndexMultimap<Index>;
+ std::shared_ptr<TrackerSimDataCollection> simData {nullptr};
+
+ if (isMC) {
+ 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;
+ }
+ ATH_MSG_VERBOSE("Found " << mcEvents->front()->particles_size() << " particles.");
+ for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
+ particles[particle->barcode()] = particle;
+ }
- // Read additional input collections
- // const auto& particles = ctx.eventStore.get<SimParticleContainer>(m_cfg.inputParticles);
- // const auto& simHits = ctx.eventStore.get<SimHitContainer>(m_cfg.inputSimHits);
- // const auto& hitParticlesMap = ctx.eventStore.get<HitParticlesMap>(m_cfg.inputMeasurementParticlesMap);
- // const auto& hitSimHitsMap = ctx.eventStore.get<HitSimHitsMap>(m_cfg.inputMeasurementSimHitsMap);
+ SG::ReadHandle<TrackerSimDataCollection> simDataHandle {m_simDataCollectionKey, ctx};
+ ATH_CHECK(simDataHandle.isValid());
+ simData = std::make_shared<TrackerSimDataCollection>(*simDataHandle);
- // For each particle within a track, how many hits did it contribute
- std::vector<ParticleHitCount> particleHitCounts;
+ SG::ReadHandle<FaserSiHitCollection> siHitCollection {m_faserSiHitKey, ctx};
+ ATH_CHECK(siHitCollection.isValid());
+ for (const FaserSiHit& hit : *siHitCollection) {
+ int barcode = hit.trackNumber();
+ Identifier id = m_idHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(), hit.getModule(), hit.getSensor());
+ siHitMap[std::make_pair(barcode, id)] = &hit;
+ }
+ }
- // Get the event number
- m_eventNr = ctx.eventID().event_number();
// Loop over the trajectories
for (size_t itraj = 0; itraj < trajectories.size(); ++itraj) {
@@ -281,28 +280,38 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
m_nStates = trajState.nStates;
// Get the majority truth particle to this track
- int barcode;
- int truthQ = 1.;
- float truthMomentum = 1;
- identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
- if (not particleHitCounts.empty()) {
- // Get the barcode of the majority truth particle
- barcode = particleHitCounts.front().particleId;
- // Find the truth particle via the barcode
- auto ip = particles.find(barcode);
- if (ip != particles.end()) {
- const auto& particle = ip->second;
- ATH_MSG_DEBUG("Find the truth particle with barcode = " << barcode);
- // Get the truth particle charge
- // FIXME find better way to access charge of simulated particle, this does not work for
- // pions which have a positive pdg code (211) and positive charge
- truthQ = particle->pdg_id() > 0 ? -1 : 1;
- truthMomentum = particle->momentum().rho() * m_MeV2GeV;
- } else {
- ATH_MSG_WARNING("Truth particle with barcode = " << barcode << " not found!");
+ int barcode = NaNint;
+ int truthQ = NaNint;
+ float truthMomentum = NaNfloat;
+ float truthLOC0 = NaNfloat;
+ float truthLOC1 = NaNfloat;
+ float truthPHI = NaNfloat;
+ float truthTHETA = NaNfloat;
+ float truthQOP = NaNfloat;
+ float truthTIME = NaNfloat;
+
+ if (isMC) {
+ truthQ = 1;
+ truthMomentum = 1;
+ identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
+ if (not particleHitCounts.empty()) {
+ // Get the barcode of the majority truth particle
+ barcode = particleHitCounts.front().particleId;
+ // Find the truth particle via the barcode
+ auto ip = particles.find(barcode);
+ if (ip != particles.end()) {
+ const auto& particle = ip->second;
+ ATH_MSG_DEBUG("Find the truth particle with barcode = " << barcode);
+ // Get the truth particle charge
+ // FIXME find better way to access charge of simulated particle, this does not work for
+ // pions which have a positive pdg code (211) and positive charge
+ truthQ = particle->pdg_id() > 0 ? -1 : 1;
+ truthMomentum = particle->momentum().rho() * m_MeV2GeV;
+ } else {
+ ATH_MSG_WARNING("Truth particle with barcode = " << barcode << " not found!");
+ }
}
}
- float truthQOP = truthQ / truthMomentum;
// Get the trackStates on the trajectory
m_nParams = {0, 0, 0};
@@ -320,48 +329,63 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
// get the truth hits corresponding to this trackState
Identifier id = state.uncalibrated().hit()->identify();
Identifier waferId = m_idHelper->wafer_id(id);
- if (siHitMap.count(std::make_pair(barcode, waferId)) == 0) {
- ATH_MSG_WARNING("no FaserSiHit for hit with id " << id << " from particle " << barcode);
- return true;
+ // if (siHitMap.count(std::make_pair(barcode, waferId)) == 0) {
+ // ATH_MSG_WARNING("no FaserSiHit for hit with id " << id << " from particle " << barcode);
+ // return true;
+ // }
+
+ if (isMC && siHitMap.count(std::make_pair(barcode, waferId)) != 0) {
+ const FaserSiHit* siHit = siHitMap.find(std::make_pair(barcode, waferId))->second;
+ HepGeom::Point3D localStartPos = siHit->localStartPosition();
+ HepGeom::Point3D localEndPos = siHit->localEndPosition();
+ HepGeom::Point3D<double> localPos = 0.5 * (localEndPos + localStartPos);
+ auto truthLocal = Acts::Vector2(localPos.y(), localPos.z());
+ const TrackerDD::SiDetectorElement* element = m_detMgr->getDetectorElement(id);
+ const HepGeom::Point3D<double> globalStartPosition =
+ Amg::EigenTransformToCLHEP(element->transformHit()) * localStartPos;
+ const HepGeom::Point3D<double> globalEndPosition =
+ Amg::EigenTransformToCLHEP(element->transformHit()) * localEndPos;
+ auto globalPosition = 0.5 * (globalStartPosition + globalEndPosition);
+ auto globalDirection = globalEndPosition - globalStartPosition;
+ auto truthUnitDir = Acts::Vector3(globalDirection.x(), globalDirection.y(), globalDirection.z()).normalized();
+ auto truthPos = Acts::Vector3(globalPosition.x() , globalPosition.y(), globalPosition.z());
+ // FIXME get truthQOP for each state
+
+ // fill the truth hit info
+ m_t_x.push_back(truthPos[Acts::ePos0]);
+ m_t_y.push_back(truthPos[Acts::ePos1]);
+ m_t_z.push_back(truthPos[Acts::ePos2]);
+ m_t_dx.push_back(truthUnitDir[Acts::eMom0]);
+ m_t_dy.push_back(truthUnitDir[Acts::eMom1]);
+ m_t_dz.push_back(truthUnitDir[Acts::eMom2]);
+
+ // get the truth track parameter at this track State
+ float truthLOC0 = truthLocal[Acts::ePos0];
+ float truthLOC1 = truthLocal[Acts::ePos1];
+ float truthTIME = siHit->meanTime();
+ float truthPHI = phi(truthUnitDir);
+ float truthTHETA = theta(truthUnitDir);
+
+ // fill the truth track parameter at this track State
+ m_t_eLOC0.push_back(truthLOC0);
+ m_t_eLOC1.push_back(truthLOC1);
+ m_t_ePHI.push_back(truthPHI);
+ m_t_eTHETA.push_back(truthTHETA);
+ m_t_eQOP.push_back(truthQOP);
+ m_t_eT.push_back(truthTIME);
+ } else {
+ m_t_x.push_back(NaNfloat);
+ m_t_y.push_back(NaNfloat);
+ m_t_z.push_back(NaNfloat);
+ m_t_dx.push_back(NaNfloat);
+ m_t_dy.push_back(NaNfloat);
+ m_t_dz.push_back(NaNfloat);
+ m_t_eLOC0.push_back(NaNfloat);
+ m_t_eLOC1.push_back(NaNfloat);
+ m_t_ePHI.push_back(NaNfloat);
+ m_t_eTHETA.push_back(NaNfloat);
+ m_t_eT.push_back(NaNfloat);
}
- const FaserSiHit* siHit = siHitMap.find(std::make_pair(barcode, waferId))->second;
- HepGeom::Point3D localStartPos = siHit->localStartPosition();
- HepGeom::Point3D localEndPos = siHit->localEndPosition();
- HepGeom::Point3D<double> localPos = 0.5 * (localEndPos + localStartPos);
- auto truthLocal = Acts::Vector2(localPos.y(), localPos.z());
- const TrackerDD::SiDetectorElement* element = m_detMgr->getDetectorElement(id);
- const HepGeom::Point3D<double> globalStartPosition =
- Amg::EigenTransformToCLHEP(element->transformHit()) * localStartPos;
- const HepGeom::Point3D<double> globalEndPosition =
- Amg::EigenTransformToCLHEP(element->transformHit()) * localEndPos;
- auto globalPosition = 0.5 * (globalStartPosition + globalEndPosition);
- auto globalDirection = globalEndPosition - globalStartPosition;
- auto truthUnitDir = Acts::Vector3(globalDirection.x(), globalDirection.y(), globalDirection.z()).normalized();
- auto truthPos = Acts::Vector3(globalPosition.x() , globalPosition.y(), globalPosition.z());
- // FIXME get truthQOP for each state
-
- // fill the truth hit info
- m_t_x.push_back(truthPos[Acts::ePos0]);
- m_t_y.push_back(truthPos[Acts::ePos1]);
- m_t_z.push_back(truthPos[Acts::ePos2]);
- m_t_dx.push_back(truthUnitDir[Acts::eMom0]);
- m_t_dy.push_back(truthUnitDir[Acts::eMom1]);
- m_t_dz.push_back(truthUnitDir[Acts::eMom2]);
-
- // get the truth track parameter at this track State
- float truthLOC0 = truthLocal[Acts::ePos0];
- float truthLOC1 = truthLocal[Acts::ePos1];
- float truthTIME = siHit->meanTime();
- float truthPHI = phi(truthUnitDir);
- float truthTHETA = theta(truthUnitDir);
-
- // fill the truth track parameter at this track State
- m_t_eLOC0.push_back(truthLOC0);
- m_t_eLOC1.push_back(truthLOC1);
- m_t_ePHI.push_back(truthPHI);
- m_t_eTHETA.push_back(truthTHETA);
- m_t_eQOP.push_back(truthQOP);
- m_t_eT.push_back(truthTIME);
// get the geometry ID
auto geoID = surface.geometryId();
@@ -433,20 +457,20 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
auto res = state.effectiveCalibrated() - H * parameters;
m_res_x_hit.push_back(res[Acts::eBoundLoc0]);
// m_res_y_hit.push_back(res[Acts::eBoundLoc1]);
- m_res_y_hit.push_back(-99.);
+ m_res_y_hit.push_back(NaNfloat);
m_err_x_hit.push_back(
sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
// m_err_y_hit.push_back(
// sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
- m_err_x_hit.push_back(-99.);
- m_res_y_hit.push_back(-99.);
+ m_err_x_hit.push_back(NaNfloat);
+ m_res_y_hit.push_back(NaNfloat);
m_pull_x_hit.push_back(
res[Acts::eBoundLoc0] /
sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
// m_pull_y_hit.push_back(
// res[Acts::eBoundLoc1] /
// sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
- m_pull_y_hit.push_back(-99);
+ m_pull_y_hit.push_back(NaNfloat);
m_dim_hit.push_back(state.calibratedSize());
}
@@ -459,52 +483,94 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
m_eT[ipar].push_back(parameters[Acts::eBoundTime]);
// track parameters residual
- m_res_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0] -
- truthLOC0);
- m_res_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1] -
- truthLOC1);
- float resPhi = Acts::detail::difference_periodic<float>(
- parameters[Acts::eBoundPhi], truthPHI,
- static_cast<float>(2 * M_PI));
- m_res_ePHI[ipar].push_back(resPhi);
- m_res_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta] -
- truthTHETA);
- m_res_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP] -
- truthQOP);
- m_res_eT[ipar].push_back(parameters[Acts::eBoundTime] - truthTIME);
-
- // track parameters error
- m_err_eLOC0[ipar].push_back(
- sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
- m_err_eLOC1[ipar].push_back(
- sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
- m_err_ePHI[ipar].push_back(
- sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
- m_err_eTHETA[ipar].push_back(
- sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
- m_err_eQOP[ipar].push_back(
- sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
- m_err_eT[ipar].push_back(
- sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
-
- // track parameters pull
- m_pull_eLOC0[ipar].push_back(
- (parameters[Acts::eBoundLoc0] - truthLOC0) /
- sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
- m_pull_eLOC1[ipar].push_back(
- (parameters[Acts::eBoundLoc1] - truthLOC1) /
- sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
- m_pull_ePHI[ipar].push_back(
- resPhi / sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
- m_pull_eTHETA[ipar].push_back(
- (parameters[Acts::eBoundTheta] - truthTHETA) /
- sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
- m_pull_eQOP[ipar].push_back(
- (parameters[Acts::eBoundQOverP] - truthQOP) /
- sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
- m_pull_eT[ipar].push_back(
- (parameters[Acts::eBoundTime] - truthTIME) /
- sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
+ float resPhi;
+ if (isMC) {
+ m_res_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0] - truthLOC0);
+ m_res_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1] - truthLOC1);
+ resPhi = Acts::detail::difference_periodic<float>( parameters[Acts::eBoundPhi], truthPHI,
+ static_cast<float>(2 * M_PI));
+ m_res_ePHI[ipar].push_back(resPhi);
+ m_res_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta] - truthTHETA);
+ m_res_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP] - truthQOP);
+ m_res_eT[ipar].push_back(parameters[Acts::eBoundTime] - truthTIME);
+
+ // track parameters error
+ m_err_eLOC0[ipar].push_back(
+ sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+ m_err_eLOC1[ipar].push_back(
+ sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
+ m_err_ePHI[ipar].push_back(
+ sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
+ m_err_eTHETA[ipar].push_back(
+ sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
+ m_err_eQOP[ipar].push_back(
+ sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
+ m_err_eT[ipar].push_back(
+ sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
+
+ // track parameters pull
+ m_pull_eLOC0[ipar].push_back(
+ (parameters[Acts::eBoundLoc0] - truthLOC0) /
+ sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+ m_pull_eLOC1[ipar].push_back(
+ (parameters[Acts::eBoundLoc1] - truthLOC1) /
+ sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
+ m_pull_ePHI[ipar].push_back(
+ resPhi / sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
+ m_pull_eTHETA[ipar].push_back(
+ (parameters[Acts::eBoundTheta] - truthTHETA) /
+ sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
+ m_pull_eQOP[ipar].push_back(
+ (parameters[Acts::eBoundQOverP] - truthQOP) /
+ sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
+ m_pull_eT[ipar].push_back(
+ (parameters[Acts::eBoundTime] - truthTIME) /
+ sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
+ } else {
+ if (ipar == 0) {
+ // push default values if no track parameters
+ m_res_x_hit.push_back(NaNfloat);
+ m_res_y_hit.push_back(NaNfloat);
+ m_err_x_hit.push_back(NaNfloat);
+ m_err_y_hit.push_back(NaNfloat);
+ m_pull_x_hit.push_back(NaNfloat);
+ m_pull_y_hit.push_back(NaNfloat);
+ m_dim_hit.push_back(NaNint);
+ }
+ // push default values if no track parameters
+ // m_eLOC0[ipar].push_back(NaNfloat);
+ // m_eLOC1[ipar].push_back(NaNfloat);
+ // m_ePHI[ipar].push_back(NaNfloat);
+ // m_eTHETA[ipar].push_back(NaNfloat);
+ // m_eQOP[ipar].push_back(NaNfloat);
+ // m_eT[ipar].push_back(NaNfloat);
+ m_res_eLOC0[ipar].push_back(NaNfloat);
+ m_res_eLOC1[ipar].push_back(NaNfloat);
+ m_res_ePHI[ipar].push_back(NaNfloat);
+ m_res_eTHETA[ipar].push_back(NaNfloat);
+ m_res_eQOP[ipar].push_back(NaNfloat);
+ m_res_eT[ipar].push_back(NaNfloat);
+ m_err_eLOC0[ipar].push_back(NaNfloat);
+ m_err_eLOC1[ipar].push_back(NaNfloat);
+ m_err_ePHI[ipar].push_back(NaNfloat);
+ m_err_eTHETA[ipar].push_back(NaNfloat);
+ m_err_eQOP[ipar].push_back(NaNfloat);
+ m_err_eT[ipar].push_back(NaNfloat);
+ m_pull_eLOC0[ipar].push_back(NaNfloat);
+ m_pull_eLOC1[ipar].push_back(NaNfloat);
+ m_pull_ePHI[ipar].push_back(NaNfloat);
+ m_pull_eTHETA[ipar].push_back(NaNfloat);
+ m_pull_eQOP[ipar].push_back(NaNfloat);
+ m_pull_eT[ipar].push_back(NaNfloat);
+ // m_x[ipar].push_back(NaNfloat);
+ // m_y[ipar].push_back(NaNfloat);
+ // m_z[ipar].push_back(NaNfloat);
+ // m_px[ipar].push_back(NaNfloat);
+ // m_py[ipar].push_back(NaNfloat);
+ // m_pz[ipar].push_back(NaNfloat);
+ // m_pT[ipar].push_back(NaNfloat);
+ // m_eta[ipar].push_back(NaNfloat);
+ }
// further track parameter info
Acts::FreeVector freeParams =
@@ -524,47 +590,47 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
} else {
if (ipar == 0) {
// push default values if no track parameters
- m_res_x_hit.push_back(-99.);
- m_res_y_hit.push_back(-99.);
- m_err_x_hit.push_back(-99.);
- m_err_y_hit.push_back(-99.);
- m_pull_x_hit.push_back(-99.);
- m_pull_y_hit.push_back(-99.);
- m_dim_hit.push_back(-99.);
+ m_res_x_hit.push_back(NaNfloat);
+ m_res_y_hit.push_back(NaNfloat);
+ m_err_x_hit.push_back(NaNfloat);
+ m_err_y_hit.push_back(NaNfloat);
+ m_pull_x_hit.push_back(NaNfloat);
+ m_pull_y_hit.push_back(NaNfloat);
+ m_dim_hit.push_back(NaNint);
}
// push default values if no track parameters
- m_eLOC0[ipar].push_back(-99.);
- m_eLOC1[ipar].push_back(-99.);
- m_ePHI[ipar].push_back(-99.);
- m_eTHETA[ipar].push_back(-99.);
- m_eQOP[ipar].push_back(-99.);
- m_eT[ipar].push_back(-99.);
- m_res_eLOC0[ipar].push_back(-99.);
- m_res_eLOC1[ipar].push_back(-99.);
- m_res_ePHI[ipar].push_back(-99.);
- m_res_eTHETA[ipar].push_back(-99.);
- m_res_eQOP[ipar].push_back(-99.);
- m_res_eT[ipar].push_back(-99.);
- m_err_eLOC0[ipar].push_back(-99);
- m_err_eLOC1[ipar].push_back(-99);
- m_err_ePHI[ipar].push_back(-99);
- m_err_eTHETA[ipar].push_back(-99);
- m_err_eQOP[ipar].push_back(-99);
- m_err_eT[ipar].push_back(-99);
- m_pull_eLOC0[ipar].push_back(-99.);
- m_pull_eLOC1[ipar].push_back(-99.);
- m_pull_ePHI[ipar].push_back(-99.);
- m_pull_eTHETA[ipar].push_back(-99.);
- m_pull_eQOP[ipar].push_back(-99.);
- m_pull_eT[ipar].push_back(-99.);
- m_x[ipar].push_back(-99.);
- m_y[ipar].push_back(-99.);
- m_z[ipar].push_back(-99.);
- m_px[ipar].push_back(-99.);
- m_py[ipar].push_back(-99.);
- m_pz[ipar].push_back(-99.);
- m_pT[ipar].push_back(-99.);
- m_eta[ipar].push_back(-99.);
+ m_eLOC0[ipar].push_back(NaNfloat);
+ m_eLOC1[ipar].push_back(NaNfloat);
+ m_ePHI[ipar].push_back(NaNfloat);
+ m_eTHETA[ipar].push_back(NaNfloat);
+ m_eQOP[ipar].push_back(NaNfloat);
+ m_eT[ipar].push_back(NaNfloat);
+ m_res_eLOC0[ipar].push_back(NaNfloat);
+ m_res_eLOC1[ipar].push_back(NaNfloat);
+ m_res_ePHI[ipar].push_back(NaNfloat);
+ m_res_eTHETA[ipar].push_back(NaNfloat);
+ m_res_eQOP[ipar].push_back(NaNfloat);
+ m_res_eT[ipar].push_back(NaNfloat);
+ m_err_eLOC0[ipar].push_back(NaNfloat);
+ m_err_eLOC1[ipar].push_back(NaNfloat);
+ m_err_ePHI[ipar].push_back(NaNfloat);
+ m_err_eTHETA[ipar].push_back(NaNfloat);
+ m_err_eQOP[ipar].push_back(NaNfloat);
+ m_err_eT[ipar].push_back(NaNfloat);
+ m_pull_eLOC0[ipar].push_back(NaNfloat);
+ m_pull_eLOC1[ipar].push_back(NaNfloat);
+ m_pull_ePHI[ipar].push_back(NaNfloat);
+ m_pull_eTHETA[ipar].push_back(NaNfloat);
+ m_pull_eQOP[ipar].push_back(NaNfloat);
+ m_pull_eT[ipar].push_back(NaNfloat);
+ m_x[ipar].push_back(NaNfloat);
+ m_y[ipar].push_back(NaNfloat);
+ m_z[ipar].push_back(NaNfloat);
+ m_px[ipar].push_back(NaNfloat);
+ m_py[ipar].push_back(NaNfloat);
+ m_pz[ipar].push_back(NaNfloat);
+ m_pT[ipar].push_back(NaNfloat);
+ m_eta[ipar].push_back(NaNfloat);
}
// fill the track parameters status
m_hasParams[ipar].push_back(hasParams[ipar]);
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
index cb19d4dd9dbdce58901a7f1a002a7f3e1a986c9e..76d76029674dcccbe48c41f5cf31b82b7292ba1b 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
@@ -15,7 +15,6 @@
#include <TTree.h>
/// NaN values for TTree variables
-constexpr double NaNdouble = std::numeric_limits<double>::quiet_NaN();
constexpr float NaNfloat = std::numeric_limits<float>::quiet_NaN();
constexpr float NaNint = std::numeric_limits<int>::quiet_NaN();
@@ -129,22 +128,27 @@ StatusCode RootTrajectorySummaryWriterTool::finalize() {
StatusCode RootTrajectorySummaryWriterTool::write(
- const Acts::GeometryContext& geoContext, const TrajectoriesContainer& trajectories) const {
+ const Acts::GeometryContext& geoContext, const TrajectoriesContainer& trajectories, bool isMC) const {
EventContext ctx = Gaudi::Hive::currentContext();
- SG::ReadHandle<TrackerSimDataCollection> simData {m_simDataCollectionKey, ctx};
- ATH_CHECK(simData.isValid());
+ std::shared_ptr<TrackerSimDataCollection> simData {nullptr};
+ std::map<int, const HepMC::GenParticle*> particles {};
- 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;
- }
+ if (isMC) {
+ SG::ReadHandle<TrackerSimDataCollection> simDataHandle {m_simDataCollectionKey, ctx};
+ ATH_CHECK(simDataHandle.isValid());
+ simData = std::make_shared<TrackerSimDataCollection>(*simDataHandle);
- std::map<int, const HepMC::GenParticle*> particles {};
- for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
- particles[particle->barcode()] = particle;
+ 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;
+ }
+
+ for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
+ particles[particle->barcode()] = particle;
+ }
}
// For each particle within a track, how many hits did it contribute
@@ -219,69 +223,71 @@ StatusCode RootTrajectorySummaryWriterTool::write(
pSurface = const_cast<Acts::Surface*>(&boundParam.referenceSurface());
}
- // Get the majority truth particle to this track
- ATH_MSG_VERBOSE("get majority truth particle");
- identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
- for (const auto& particle : particleHitCounts) {
- ATH_MSG_VERBOSE(particle.particleId << ": " << particle.hitCount << " hits");
- }
+ if (isMC) {
+ // Get the majority truth particle to this track
+ ATH_MSG_VERBOSE("get majority truth particle");
+ identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
+ for (const auto& particle : particleHitCounts) {
+ ATH_MSG_VERBOSE(particle.particleId << ": " << particle.hitCount << " hits");
+ }
- bool foundMajorityParticle = false;
- // Get the truth particle info
- if (not particleHitCounts.empty()) {
- // Get the barcode of the majority truth particle
- majorityParticleId = particleHitCounts.front().particleId;
- nMajorityHits = particleHitCounts.front().hitCount;
-
- // Find the truth particle via the barcode
- auto ip = particles.find(majorityParticleId);
- if (ip != particles.end()) {
- foundMajorityParticle = true;
-
- const HepMC::GenParticle* particle = ip->second;
- ATH_MSG_DEBUG("Find the truth particle with barcode = " << majorityParticleId);
-
- // extrapolate parameters from vertex to reference surface at origin.
- std::unique_ptr<const Acts::BoundTrackParameters> truthParameters
- = extrapolateToReferenceSurface(ctx, particle);
- if (!truthParameters) {
- continue;
- }
- // Get the truth particle info at vertex
- // const HepMC::GenVertex* vertex = particle->production_vertex();
- t_p = truthParameters->momentum().mag();
- t_charge = truthParameters->charge();
- t_time = truthParameters->time();
- t_vx = truthParameters->position(geoContext).x();
- t_vy = truthParameters->position(geoContext).y();
- t_vz = truthParameters->position(geoContext).z();
- t_px = truthParameters->momentum().x();
- t_py = truthParameters->momentum().y();
- t_pz = truthParameters->momentum().z();
- t_theta = theta(truthParameters->momentum().normalized());
- t_phi = phi(truthParameters->momentum().normalized());
- t_eta = eta(truthParameters->momentum().normalized());
- t_pT = t_p * perp(truthParameters->momentum().normalized());
-
- auto unitDirection = Acts::Vector3(t_px, t_py, t_pz).normalized();
- auto positon = Acts::Vector3 (t_vx, t_vy, t_vz);
-
- if (pSurface) {
- // get the truth perigee parameter
- auto lpResult = pSurface->globalToLocal(geoContext, positon, unitDirection);
- if (lpResult.ok()) {
- t_vlx = lpResult.value()[Acts::BoundIndices::eBoundLoc0];
- t_vly = lpResult.value()[Acts::BoundIndices::eBoundLoc1];
- } else {
- ATH_MSG_ERROR("Global to local transformation did not succeed.");
+ bool foundMajorityParticle = false;
+ // Get the truth particle info
+ if (not particleHitCounts.empty()) {
+ // Get the barcode of the majority truth particle
+ majorityParticleId = particleHitCounts.front().particleId;
+ nMajorityHits = particleHitCounts.front().hitCount;
+
+ // Find the truth particle via the barcode
+ auto ip = particles.find(majorityParticleId);
+ if (ip != particles.end()) {
+ foundMajorityParticle = true;
+
+ const HepMC::GenParticle* particle = ip->second;
+ ATH_MSG_DEBUG("Find the truth particle with barcode = " << majorityParticleId);
+
+ // extrapolate parameters from vertex to reference surface at origin.
+ std::unique_ptr<const Acts::BoundTrackParameters> truthParameters
+ = extrapolateToReferenceSurface(ctx, particle);
+ if (!truthParameters) {
+ continue;
}
+ // Get the truth particle info at vertex
+ // const HepMC::GenVertex* vertex = particle->production_vertex();
+ t_p = truthParameters->momentum().mag();
+ t_charge = truthParameters->charge();
+ t_time = truthParameters->time();
+ t_vx = truthParameters->position(geoContext).x();
+ t_vy = truthParameters->position(geoContext).y();
+ t_vz = truthParameters->position(geoContext).z();
+ t_px = truthParameters->momentum().x();
+ t_py = truthParameters->momentum().y();
+ t_pz = truthParameters->momentum().z();
+ t_theta = theta(truthParameters->momentum().normalized());
+ t_phi = phi(truthParameters->momentum().normalized());
+ t_eta = eta(truthParameters->momentum().normalized());
+ t_pT = t_p * perp(truthParameters->momentum().normalized());
+
+ auto unitDirection = Acts::Vector3(t_px, t_py, t_pz).normalized();
+ auto positon = Acts::Vector3 (t_vx, t_vy, t_vz);
+
+ if (pSurface) {
+ // get the truth perigee parameter
+ auto lpResult = pSurface->globalToLocal(geoContext, positon, unitDirection);
+ if (lpResult.ok()) {
+ t_vlx = lpResult.value()[Acts::BoundIndices::eBoundLoc0];
+ t_vly = lpResult.value()[Acts::BoundIndices::eBoundLoc1];
+ } else {
+ ATH_MSG_ERROR("Global to local transformation did not succeed.");
+ }
+ }
+ } else {
+ ATH_MSG_WARNING("Truth particle with barcode = " << majorityParticleId << " not found in the input collection!");
}
- } else {
- ATH_MSG_WARNING("Truth particle with barcode = " << majorityParticleId << " not found in the input collection!");
}
- }
- if (not foundMajorityParticle) {
- ATH_MSG_WARNING("Truth particle for mj " << itraj << " subtraj " << isubtraj << " not found!");
+ if (not foundMajorityParticle) {
+ ATH_MSG_WARNING("Truth particle for mj " << itraj << " subtraj " << isubtraj << " not found!");
+ }
}
// Push the corresponding truth particle info for the track.
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/SeedingAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/SeedingAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5325703b0ba3599a1fdcdafb350850ddc621c555
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/SeedingAlg.cxx
@@ -0,0 +1,22 @@
+#include "FaserActsKalmanFilter/SeedingAlg.h"
+
+SeedingAlg::SeedingAlg(const std::string& name, ISvcLocator* pSvcLocator) :
+ AthAlgorithm(name, pSvcLocator) {}
+
+
+StatusCode SeedingAlg::initialize() {
+ ATH_CHECK(m_trackSeedTool.retrieve());
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode SeedingAlg::execute() {
+ // const EventContext& ctx = Gaudi::Hive::currentContext();
+ ATH_CHECK(m_trackSeedTool->run());
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode SeedingAlg::finalize() {
+ return StatusCode::SUCCESS;
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
index 785bf2be17d9fe0c612b0250c6eae67ae801626e..6b780fc368e874f4520af430e63d1cc7c4e02575 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
@@ -45,18 +45,49 @@ void identifyContributingParticles(
if (not state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
return true;
}
+ std::vector<int> barcodes {};
// register all particles that generated this hit
- std::vector<Identifier> ids = state.uncalibrated().hit()->rdoList();
- for (const Identifier &id : ids) {
+ for (const Identifier &id : state.uncalibrated().hit()->rdoList()) {
if (simDataCollection.count(id) == 0) {
return true;
}
const auto &deposits = simDataCollection.find(id)->second.getdeposits();
for (const TrackerSimData::Deposit &deposit : deposits) {
- increaseHitCount(particleHitCounts, deposit.first->barcode());
+ int barcode = deposit.first->barcode();
+ if (std::find(barcodes.begin(), barcodes.end(), barcode) == barcodes.end()) {
+ barcodes.push_back(barcode);
+ increaseHitCount(particleHitCounts, deposit.first->barcode());
+ }
}
}
return true;
});
sortHitCount(particleHitCounts);
}
+
+/* Identify all particles that contribute to a trajectory.
+ * If a cluster consists of multiple RDOs we check for each from which particle(s) it has been created.
+ * And if multiple particles created a RDO we increase the hit count for each of them.
+ */
+void identifyContributingParticles(
+ const TrackerSimDataCollection& simDataCollection,
+ const std::vector<const Tracker::FaserSCT_Cluster*> clusters,
+ std::vector<ParticleHitCount>& particleHitCounts) {
+ particleHitCounts.clear();
+ for (const Tracker::FaserSCT_Cluster *cluster : clusters) {
+ std::vector<int> barcodes {};
+ for (const Identifier &id : cluster->rdoList()) {
+ if (simDataCollection.count(id) == 0) continue;
+ const auto &deposits = simDataCollection.find(id)->second.getdeposits();
+ for (const TrackerSimData::Deposit &deposit : deposits) {
+ int barcode = deposit.first->barcode();
+ if (std::find(barcodes.begin(), barcodes.end(), barcode) == barcodes.end()) {
+ barcodes.push_back(barcode);
+ increaseHitCount(particleHitCounts, deposit.first->barcode());
+ }
+ }
+ }
+ }
+ sortHitCount(particleHitCounts);
+}
+
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
index 8d56370abefe0f7399095ac3b62d5bc387ad2f35..4c23cb7febdd977df01a7a4a4de997267dbfd4c6 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
@@ -25,7 +25,9 @@
#include "FaserActsKalmanFilter/CKF2.h"
#include "FaserActsKalmanFilter/KalmanFitterTool.h"
#include "FaserActsKalmanFilter/MyTrackSeedTool.h"
-
+#include "FaserActsKalmanFilter/SeedingAlg.h"
+#include "FaserActsKalmanFilter/CircleFitTrackSeedTool.h"
+#include "FaserActsKalmanFilter/GhostBusters.h"
DECLARE_COMPONENT(FaserActsKalmanFilterAlg)
DECLARE_COMPONENT(CombinatorialKalmanFilterAlg)
@@ -50,3 +52,6 @@ DECLARE_COMPONENT(ActsTrackSeedTool)
DECLARE_COMPONENT(CKF2)
DECLARE_COMPONENT(KalmanFitterTool)
DECLARE_COMPONENT(MyTrackSeedTool)
+DECLARE_COMPONENT(SeedingAlg)
+DECLARE_COMPONENT(CircleFitTrackSeedTool)
+DECLARE_COMPONENT(GhostBusters)
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/CKF2.py b/Tracking/Acts/FaserActsKalmanFilter/test/CKF2.py
index 5d41a03a33273d3003677932920b059df871b3d5..98059c9e65857ac9f702615daa263b9c307f7877 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/test/CKF2.py
+++ b/Tracking/Acts/FaserActsKalmanFilter/test/CKF2.py
@@ -7,32 +7,37 @@ from AthenaCommon.Configurable import Configurable
from CalypsoConfiguration.AllConfigFlags import ConfigFlags
from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
-# from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.GhostBustersConfig import GhostBustersCfg
from FaserActsKalmanFilter.CKF2Config import CKF2Cfg
log.setLevel(DEBUG)
Configurable.configurableRun3Behavior = True
ConfigFlags.Input.Files = ['my.RDO.pool.root']
+ConfigFlags.addFlag("Output.xAODFileName", f"CKF.xAOD.root")
ConfigFlags.Output.ESDFileName = "CKF.ESD.pool.root"
-ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
-ConfigFlags.GeoModel.FaserVersion = "FASER-01"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02"
+ConfigFlags.GeoModel.FaserVersion = "FASERNU-03"
ConfigFlags.GeoModel.Align.Dynamic = False
ConfigFlags.Beam.NumberOfCollisions = 0.
ConfigFlags.TrackingGeometry.MaterialSource = "Input"
+ConfigFlags.Input.isMC = True
ConfigFlags.lock()
acc = MainServicesCfg(ConfigFlags)
acc.merge(PoolReadCfg(ConfigFlags))
-# acc.merge(PoolWriteCfg(ConfigFlags))
+acc.merge(PoolWriteCfg(ConfigFlags))
acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
-acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.51, MinClustersPerFit=5, TanThetaCut=0.25))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.61, MinClustersPerFit=5, TanThetaXZCut=0.083))
+acc.merge(GhostBustersCfg(ConfigFlags))
acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=True))
+# acc.getEventAlgo("CKF2").OutputLevel = DEBUG
# logging.getLogger('forcomps').setLevel(VERBOSE)
# acc.foreach_component("*").OutputLevel = VERBOSE
@@ -42,5 +47,17 @@ acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=True))
# acc.printConfig(withDetails=True)
# ConfigFlags.dump()
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+itemList = [
+ "xAOD::EventInfo#*",
+ "xAOD::EventAuxInfo#*",
+ "xAOD::FaserTriggerData#*",
+ "xAOD::FaserTriggerDataAux#*",
+ "FaserSCT_RDO_Container#*",
+ "Tracker::FaserSCT_ClusterContainer#*",
+ "TrackCollection#*",
+]
+acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList))
+
sc = acc.run(maxEvents=-1)
sys.exit(not sc.isSuccess())
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/GhostBusters.py b/Tracking/Acts/FaserActsKalmanFilter/test/GhostBusters.py
new file mode 100644
index 0000000000000000000000000000000000000000..297f76e61fd138c325b67e98d0ea195e062eb5f3
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/test/GhostBusters.py
@@ -0,0 +1,48 @@
+#!/usr/bin/env python
+
+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 CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.SeedingConfig import SeedingCfg
+from FaserActsKalmanFilter.GhostBustersConfig import GhostBustersCfg
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+ConfigFlags.Input.Files = ['my.RDO.pool.root']
+ConfigFlags.Output.ESDFileName = "ghosts.ESD.pool.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+ConfigFlags.GeoModel.FaserVersion = "FASER-01"
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+# ConfigFlags.TrackingGeometry.MaterialSource = "Input"
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+acc.merge(PoolWriteCfg(ConfigFlags))
+acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.61, MinClustersPerFit=5, TanThetaXZCut=0.083))
+# acc.getEventAlgo("Tracker::SegmentFitAlg").OutputLevel = VERBOSE
+acc.merge(GhostBustersCfg(ConfigFlags, xTolerance=0.5, yTolerance=0.5))
+acc.getEventAlgo("GhostBusters").OutputLevel = DEBUG
+
+# logging.getLogger('forcomps').setLevel(VERBOSE)
+# acc.foreach_component("*").OutputLevel = VERBOSE
+# acc.foreach_component("*ClassID*").OutputLevel = VERBOSE
+# acc.getService("StoreGateSvc").Dump = True
+# acc.getService("ConditionStore").Dump = True
+# acc.printConfig(withDetails=True)
+# ConfigFlags.dump()
+
+sc = acc.run(maxEvents=-1)
+sys.exit(not sc.isSuccess())
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/SeedingDbg.py b/Tracking/Acts/FaserActsKalmanFilter/test/SeedingDbg.py
new file mode 100644
index 0000000000000000000000000000000000000000..01e04b803a6f1ed60ddd54db02fcea952ad8849c
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/test/SeedingDbg.py
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+
+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 CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.SeedingConfig import SeedingCfg
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+ConfigFlags.Input.Files = ['my.RDO.pool.root']
+ConfigFlags.Output.ESDFileName = "circleFitSeeding.ESD.pool.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+ConfigFlags.GeoModel.FaserVersion = "FASER-01"
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+# ConfigFlags.TrackingGeometry.MaterialSource = "Input"
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+acc.merge(PoolWriteCfg(ConfigFlags))
+acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.51, MinClustersPerFit=5, TanThetaCut=0.25))
+acc.merge(SeedingCfg(ConfigFlags))
+acc.getEventAlgo("SeedingAlg").OutputLevel = VERBOSE
+
+# logging.getLogger('forcomps').setLevel(VERBOSE)
+# acc.foreach_component("*").OutputLevel = VERBOSE
+# acc.foreach_component("*ClassID*").OutputLevel = VERBOSE
+# acc.getService("StoreGateSvc").Dump = True
+# acc.getService("ConditionStore").Dump = True
+# acc.printConfig(withDetails=True)
+# ConfigFlags.dump()
+
+sc = acc.run(maxEvents=10)
+sys.exit(not sc.isSuccess())
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/TI12CKF2.py b/Tracking/Acts/FaserActsKalmanFilter/test/TI12CKF2.py
new file mode 100644
index 0000000000000000000000000000000000000000..f6c5a5af53581d986d57ee32ad2e4d5706c6f6e1
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/test/TI12CKF2.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+
+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 CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+# from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.GhostBustersConfig import GhostBustersCfg
+from FaserActsKalmanFilter.TI12CKF2Config import TI12CKF2Cfg
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+ConfigFlags.Input.Files = ['threeStationRun6833.raw']
+ConfigFlags.Output.ESDFileName = "CKF.ESD.pool.root"
+ConfigFlags.addFlag("Output.xAODFileName", f"CKF.xAOD.root")
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02"
+ConfigFlags.IOVDb.DatabaseInstance = "OFLP200"
+ConfigFlags.Input.ProjectName = "data22"
+ConfigFlags.Input.isMC = False
+ConfigFlags.GeoModel.FaserVersion = "FASER-02"
+ConfigFlags.Common.isOnline = False
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+ConfigFlags.Detector.GeometryFaserSCT = True
+ConfigFlags.TrackingGeometry.MaterialSource = "Input"
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolWriteCfg(ConfigFlags))
+acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags))
+acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="SCT_LEVELMODE_RDOs", ClusterToolTimingPattern="X1X"))
+acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.61, MinClustersPerFit=5, TanThetaXZCut=0.083))
+acc.merge(GhostBustersCfg(ConfigFlags))
+acc.merge(TI12CKF2Cfg(ConfigFlags, noDiagnostics=False))
+acc.getEventAlgo("CKF2").OutputLevel = DEBUG
+
+# logging.getLogger('forcomps').setLevel(VERBOSE)
+# acc.foreach_component("*").OutputLevel = VERBOSE
+# acc.foreach_component("*ClassID*").OutputLevel = VERBOSE
+# acc.getService("StoreGateSvc").Dump = True
+# acc.getService("ConditionStore").Dump = True
+# acc.printConfig(withDetails=True)
+# ConfigFlags.dump()
+
+# 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
+
+sc = acc.run(maxEvents=-1)
+sys.exit(not sc.isSuccess())
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/TI12_CKF.py b/Tracking/Acts/FaserActsKalmanFilter/test/TI12_CKF.py
deleted file mode 100644
index 98acb56f00ace95d54277726768250e9543b06ab..0000000000000000000000000000000000000000
--- a/Tracking/Acts/FaserActsKalmanFilter/test/TI12_CKF.py
+++ /dev/null
@@ -1,56 +0,0 @@
-#!/usr/bin/env python
-
-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 CalypsoConfiguration.MainServicesConfig import MainServicesCfg
-from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg
-from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
-from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
-from FaserActsKalmanFilter.CombinatorialKalmanFilterConfig import CombinatorialKalmanFilterCfg
-
-log.setLevel(DEBUG)
-Configurable.configurableRun3Behavior = True
-
-ConfigFlags.Input.Files = ['/home/tboeckh/tmp/Faser-Physics-006470-00093.raw_middleStation.SPs']
-ConfigFlags.Output.ESDFileName = "CKF.ESD.pool.root"
-ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02"
-ConfigFlags.GeoModel.FaserVersion = "FASER-02"
-ConfigFlags.IOVDb.DatabaseInstance = "OFLP200"
-ConfigFlags.Input.ProjectName = "data21"
-ConfigFlags.Input.isMC = False
-ConfigFlags.Common.isOnline = False
-ConfigFlags.GeoModel.Align.Dynamic = False
-ConfigFlags.Beam.NumberOfCollisions = 0.
-ConfigFlags.Detector.GeometryFaserSCT = True
-ConfigFlags.lock()
-
-acc = MainServicesCfg(ConfigFlags)
-acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags))
-acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, name="LevelClustering", DataObjectName="SCT_LEVELMODE_RDOs", ClusterToolTimingPattern="X1X"))
-acc.merge(SegmentFitAlgCfg(ConfigFlags, name=f"LevelFit", MaxClusters=44))
-# acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="SCT_RDOs", ClusterToolTimingPattern="XXX"))
-# acc.merge(SegmentFitAlgCfg(ConfigFlags))
-acc.merge(CombinatorialKalmanFilterCfg(ConfigFlags))
-acc.getEventAlgo("CombinatorialKalmanFilterAlg").OutputLevel = VERBOSE
-
-replicaSvc = acc.getService("DBReplicaSvc")
-replicaSvc.COOLSQLiteVetoPattern = ""
-replicaSvc.UseCOOLSQLite = True
-replicaSvc.UseCOOLFrontier = False
-replicaSvc.UseGeomSQLite = True
-
-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()
-
-sc = acc.run(maxEvents=-1)
-sys.exit(not sc.isSuccess())