diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 0aef82debb458e44c3b9bfa1dceb4d2f79548bee..521e9b519646214a6154fc49f8195c60ed6806b2 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -16,10 +16,11 @@ build_image:
# description: triggers a build of the project as a Docker image,
# each branch will have an individual Docker image that will be used
# in the following stages of the pipeline for testing the code
+ image:
+ name: gitlab-registry.cern.ch/linuxsupport/cc7-base:latest
stage: build
tags:
- - cvmfs
- - docker
+ - k8s-cvmfs
script:
- yum -y install redhat-lsb redhat-lsb-core man uuid-devel libuuid libuuid-devel mesa-libGL-devel libXpm-devel
- mkdir build
@@ -33,12 +34,13 @@ build_image:
- build/
test_unittest:
+ image:
+ name: gitlab-registry.cern.ch/linuxsupport/cc7-base:latest
stage: test
tags:
- - cvmfs
- - docker
+ - k8s-cvmfs
script:
- - yum -y install man
+ - yum -y install man git make cmake3 gcc-c++ gcc binutils libX11-devel libXpm-devel libXft-devel libXext-devel python openssl-devel
- cd build
- set +e && source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh; set -e
- set +e && asetup --input=../../calypso/asetup.faser Athena,22.0.49; set -e
@@ -46,3 +48,6 @@ test_unittest:
- ctest -j3
dependencies:
- build_image
+ artifacts:
+ paths:
+ - LastTest.log
diff --git a/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py b/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py
index 934433e2f375d6f0883a5a66e66dba70f1515427..8986c1165ece677bb37a42998bc2e0d05642cd8e 100644
--- a/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py
+++ b/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py
@@ -36,4 +36,4 @@ class TruthEventDumperAlg(PyAthena.Alg):
# print("Mid part: ", self.maxMid, " out of ", 1000000 - 200000, " (",100*self.maxMid/(1000000-200000),"% of overflow")
# print("Hi part: ", self.maxHi, " out of ", (1<<31)//1000000, " (", 100*self.maxHi/((1<<31)//1000000),"% of overflow")
- return StatusCode.Success
\ No newline at end of file
+ return StatusCode.Success
diff --git a/DetectorDescription/GeoModel/FaserGeoModel/data/geomDB.sql b/DetectorDescription/GeoModel/FaserGeoModel/data/geomDB.sql
index 046c5462f7dde50bdaa43a4851b9b303c82a4770..d4f4890493841ab030f6570cecb52742fba53f69 100644
--- a/DetectorDescription/GeoModel/FaserGeoModel/data/geomDB.sql
+++ b/DetectorDescription/GeoModel/FaserGeoModel/data/geomDB.sql
@@ -1038,6 +1038,8 @@ INSERT INTO "HVS_TAG2NODE" VALUES (41, "Ecal-TB01", 107857, NULL, 0, 0, 16858
INSERT INTO "HVS_TAG2NODE" VALUES (410, "EcalTopLevel-00", 100058, NULL, 0, 0, 1599350400000000000, NULL, 22);
INSERT INTO "HVS_TAG2NODE" VALUES (410, "EcalTopLevel-02", 107797, NULL, 0, 0, 1619222400000000000, NULL, 22);
INSERT INTO "HVS_TAG2NODE" VALUES (410, "EcalTopLevel-TB00", 107822, NULL, 0, 0, 1627862400000000000, NULL, 22);
+INSERT INTO "HVS_TAG2NODE" VALUES (410, "EcalTopLevel-04", 107886, NULL, 0, 0, 1695600000000000000, NULL, 22);
+INSERT INTO "HVS_TAG2NODE" VALUES (410, "EcalTopLevel-TB01", 107887, NULL, 0, 0, 1695600000000000000, NULL, 22);
INSERT INTO "HVS_TAG2NODE" VALUES (411, "EcalRowGeneral-00", 100059, NULL, 0, 0, 1599350400000000000, NULL, 22);
INSERT INTO "HVS_TAG2NODE" VALUES (411, "EcalRowGeneral-TB00", 107830, NULL, 0, 0, 1627862400000000000, NULL, 22);
INSERT INTO "HVS_TAG2NODE" VALUES (414, "EcalSwitches-00", 100057, NULL, 0, 0, 1593907200000000000, NULL, 22);
@@ -2009,7 +2011,7 @@ INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "DipoleSwitches", "DipoleSw
INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "DipoleWrappingGeneral", "DipoleWrappingGeneral-01", 107885);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "Calorimeter", "Calorimeter-04", 107858);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "Ecal", "Ecal-04", 107856);
-INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "EcalTopLevel", "EcalTopLevel-02", 107797);
+INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "EcalTopLevel", "EcalTopLevel-04", 107886);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "EcalRowGeneral", "EcalRowGeneral-00", 100059);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "EcalSwitches", "EcalSwitches-01", 107855);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASERNU-04", "VetoTopLevel", "VetoTopLevel-02", 107798);
@@ -2087,7 +2089,7 @@ INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "DipoleSwitches", "DipoleSw
INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "DipoleWrappingGeneral", "DipoleWrappingGeneral-01", 107885);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "Calorimeter", "Calorimeter-TB01", 107859);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "Ecal", "Ecal-TB01", 107857);
-INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "EcalTopLevel", "EcalTopLevel-TB00", 107822);
+INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "EcalTopLevel", "EcalTopLevel-TB00", 107887);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "EcalRowGeneral", "EcalRowGeneral-TB00", 107830);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "EcalSwitches", "EcalSwitches-01", 107855);
INSERT INTO "HVS_TAGCACHE" VALUES ("FASER-TB01", "VetoTopLevel", "VetoTopLevel-TB00", 107823);
@@ -2856,13 +2858,17 @@ INSERT INTO "PRESHOWERWRAPPINGGENERAL_DATA2TAG" VALUES (107875, 0);
INSERT INTO "PRESHOWERWRAPPINGGENERAL_DATA2TAG" VALUES (107876, 1);
--
--
-INSERT INTO "ECALTOPLEVEL_DATA" VALUES (0, 0.0, 0.0, 3099.2, 0.0, 0.0, 0.0, 321, "Ecal");
-INSERT INTO "ECALTOPLEVEL_DATA" VALUES (1, 0.0, -71.6, 0.0, 0.0, 2.8, 0.0, 321, "BottomRow");
-INSERT INTO "ECALTOPLEVEL_DATA" VALUES (2, 0.0, 49.6, 0.0, 0.0, 2.8, 0.0, 321, "TopRow");
-INSERT INTO "ECALTOPLEVEL_DATA" VALUES (3, 0.0, 0.0, 3098.9, 0.0, 0.0, 0.0, 321, "Ecal");
-INSERT INTO "ECALTOPLEVEL_DATA" VALUES (4, 0.0, -71.6, 0.0, 0.0, 0.0, 0.0, 321, "BottomRow");
-INSERT INTO "ECALTOPLEVEL_DATA" VALUES (5, 0.0, 49.6, 0.0, 0.0, 0.0, 0.0, 321, "TopRow");
-INSERT INTO "ECALTOPLEVEL_DATA" VALUES (6, -17.22, -22.80, 1088.66, -2.8, -2.8, 0.0, 321, "Ecal");
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (0, 0.0, 0.0, 3099.2, 0.0, 0.0, 0.0, 321, "Ecal"); -- Old TI12
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (1, 0.0, -71.6, 0.0, 0.0, 2.8, 0.0, 321, "BottomRow"); -- TI12 (with overlap)
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (2, 0.0, 49.6, 0.0, 0.0, 2.8, 0.0, 321, "TopRow"); -- TI12 (with overlap)
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (3, 0.0, 0.0, 3098.9, 0.0, 0.0, 0.0, 321, "Ecal"); -- TI12
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (4, 0.0, -71.6, 0.0, 0.0, 0.0, 0.0, 321, "BottomRow"); -- TB00 (with overlap)
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (5, 0.0, 49.6, 0.0, 0.0, 0.0, 0.0, 321, "TopRow"); -- TB00 (with overlap)
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (6, -17.22, -22.80, 1088.66, -2.8, -2.8, 0.0, 321, "Ecal"); -- test beam
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (7, 0.0, -72.0, 0.0, 0.0, 2.8, 0.0, 321, "BottomRow"); -- TI12 (overlap fixed)
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (8, 0.0, 50.0, 0.0, 0.0, 2.8, 0.0, 321, "TopRow"); -- TI12 (overlap fixed)
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (9, 0.0, -72.0, 0.0, 0.0, 0.0, 0.0, 321, "BottomRow"); -- TB01 (overlap fixed)
+INSERT INTO "ECALTOPLEVEL_DATA" VALUES (10, 0.0, 50.0, 0.0, 0.0, 0.0, 0.0, 321, "TopRow"); -- TB01 (overlap fixed)
INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (100058, 0);
INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (100058, 1);
INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (100058, 2);
@@ -2872,6 +2878,13 @@ INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107797, 3);
INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107822, 4);
INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107822, 5);
INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107822, 6);
+INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107886, 3);
+INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107886, 7);
+INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107886, 8);
+INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107887, 6);
+INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107887, 9);
+INSERT INTO "ECALTOPLEVEL_DATA2TAG" VALUES (107887, 10);
+
--
--
INSERT INTO "ECALROWGENERAL_DATA" VALUES (0, 2, 0.0, 2.8);
diff --git a/Generators/TestBeamReader/CMakeLists.txt b/Generators/TestBeamReader/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..64e5ec13d29414fcbf2468e64170029a8cf6c44b
--- /dev/null
+++ b/Generators/TestBeamReader/CMakeLists.txt
@@ -0,0 +1,11 @@
+################################################################################
+# Package: ReadTestBeam
+################################################################################
+
+# Declare the package name:
+atlas_subdir( TestBeamReader )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
+
+atlas_install_joboptions( share/*.py )
\ No newline at end of file
diff --git a/Generators/TestBeamReader/python/TestBeamReaderAlgH2_CMS.py b/Generators/TestBeamReader/python/TestBeamReaderAlgH2_CMS.py
new file mode 100644
index 0000000000000000000000000000000000000000..959242ab0aaa0639c672ba669b9c16400f5f726d
--- /dev/null
+++ b/Generators/TestBeamReader/python/TestBeamReaderAlgH2_CMS.py
@@ -0,0 +1,95 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+from GeneratorModules.EvgenAlg import EvgenAlg
+from AthenaPython.PyAthena import StatusCode, EventInfo, EventID, EventType
+from AthenaCommon.SystemOfUnits import GeV, m
+import ROOT
+import math
+
+
+class TestBeamReader(EvgenAlg):
+ def __init__(self, name="TestBeamReader", MCEventKey="BeamTruthEvent"):
+ super(TestBeamReader,self).__init__(name=name)
+ self.McEventKey = MCEventKey
+ return
+
+ def fillEvent(self, evt):
+ try:
+ from AthenaPython.PyAthena import HepMC3 as HepMC
+ except ImportError:
+ from AthenaPython.PyAthena import HepMC as HepMC
+ evt.weights().push_back(1.0)
+ eventNumber = 0
+ treeEventId = 0
+ runNumber = 100
+ mcEventType = EventType()
+ mcEventType.add_type(EventType.IS_SIMULATION)
+
+ mcEventId = EventID(run_number = runNumber, event_number = eventNumber)
+ mcEventInfo = EventInfo(id = mcEventId, type = mcEventType)
+ self.evtStore.record(mcEventInfo, "McEventInfo", True, False)
+ ROOT.SetOwnership(mcEventType, False)
+ ROOT.SetOwnership(mcEventId, False)
+ ROOT.SetOwnership(mcEventInfo, False)
+ pos = HepMC.FourVector(13.8*0.001*m, 22.7*0.001*m, 0*m, 0)
+ gv = HepMC.GenVertex(pos)
+ ROOT.SetOwnership(gv, False)
+ evt.add_vertex(gv)
+
+ nParticles = list(self.evtStore["N"])
+
+ pdgclol = list(self.evtStore["vecPDG"])
+ px =list(self.evtStore["vecPx"])
+ py = list(self.evtStore["vecPy"])
+ pz = list(self.evtStore["vecPz"])
+
+ for i in range(nParticles[0]):
+ gp = HepMC.GenParticle()
+ pdgc = int()
+
+ pdgc = math.floor(pdgclol[i])
+ if (pdgc ==-11 or pdgc ==11):
+ M = (0.511*0.001)
+ elif (pdgc == 22):
+ M = 0.0
+ else:
+ M = 0.0
+ E = math.sqrt(((px[i])*(px[i]))+((py[i])*(py[i]))+((pz[i])*(pz[i]))+((M)*(M)))
+ mom = HepMC.FourVector(px[i]*0.001*GeV, py[i]*0.001*GeV, pz[i]*0.001*GeV, E*0.001*GeV)
+ gp.set_momentum(mom)
+ gp.set_generated_mass(M*GeV)
+ gp.set_pdg_id(pdgc)
+
+ hepmc_status = 4
+ gp.set_status(hepmc_status)
+ ROOT.SetOwnership(gp, False)
+ gv.add_particle_in(gp)
+
+
+
+ for i in range(nParticles[0]):
+ gp = HepMC.GenParticle()
+
+ pdgc = int()
+ pdgc = math.floor(pdgclol[i])
+ if (pdgc ==-11 or pdgc ==11):
+ M = (0.511*0.001)
+ elif (pdgc == 22):
+ M = 0.0
+ else:
+ M = 0.0
+ E = math.sqrt(((px[i])*(px[i]))+((py[i])*(py[i]))+((pz[i])*(pz[i]))+((M)*(M)))
+ mom = HepMC.FourVector(px[i]*0.001*GeV, py[i]*0.001*GeV, pz[i]*0.001*GeV, E*0.001*GeV)
+ gp.set_momentum(mom)
+ gp.set_generated_mass(M*GeV)
+ gp.set_pdg_id(pdgc)
+
+ hepmc_status = 1
+ gp.set_status(hepmc_status)
+ ROOT.SetOwnership(gp, False)
+
+ gv.add_particle_out(gp)
+
+
+ return StatusCode.Success
diff --git a/Generators/TestBeamReader/python/__init__.py b/Generators/TestBeamReader/python/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..d13ae164caa4e93bf2899cea4e67d0ec515784a2
--- /dev/null
+++ b/Generators/TestBeamReader/python/__init__.py
@@ -0,0 +1 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
diff --git a/Generators/TestBeamReader/share/Numpy_H2.py b/Generators/TestBeamReader/share/Numpy_H2.py
new file mode 100644
index 0000000000000000000000000000000000000000..9644c3d584d8157b8aa6b1acdcba8d1b7fad5b32
--- /dev/null
+++ b/Generators/TestBeamReader/share/Numpy_H2.py
@@ -0,0 +1,90 @@
+#!/usr/bin/env python
+
+# Set up (Py)ROOT.
+import ROOT
+import datetime
+from array import array
+import math
+import datetime
+import pickle
+import pandas
+import uproot3
+import uproot
+import numpy as np
+import awkward as ak
+t1 = uproot3.open("/eos/project/f/faser-preshower/simulations/H2_beamline/1M_events/FASER_e-_v46_80GeV.root")["NTuples/GoodParticle"]
+
+
+FASERdf= t1.pandas.df(["FASER_Px", "FASER_Py","FASER_Pz" ,"FASER_PDGid", "FASER_EventID", "FASER_x", "FASER_y"] )
+pandas.set_option("display.max_rows", None, "display.max_columns", None)
+
+
+FASERdf_new = FASERdf.groupby('FASER_EventID').agg(list)
+
+
+FASERdf_new =FASERdf_new.drop(FASERdf_new.index[0])
+
+
+data = {key: FASERdf_new[key].values for key in [ "FASER_Px", "FASER_Py","FASER_Pz" ,"FASER_PDGid", "FASER_x", "FASER_y"]}
+
+a = ak.Array(data["FASER_Px"])
+print(a)
+#ta =a
+b = ak.Array(data["FASER_Py"])
+print(b)
+c = ak.Array(data["FASER_Pz"])
+print(c)
+d = ak.Array(data["FASER_PDGid"])
+print(d)
+e = ak.Array(data["FASER_x"])
+print(e)
+f = ak.Array(data["FASER_y"])
+print(f)
+
+with uproot.recreate("80GeV_first.root") as rootfile:
+ rootfile["tree1"] = ak.zip({"FASER_Px": a, "FASER_Py": b, "FASER_Pz": c, "FASER_PDGid":d, "FASER_x":e, "FASER_y": f})
+
+f = ROOT.TFile.Open("80GeV_e-_calypso_H2input.root", "recreate")
+t = ROOT.TTree("t","t")
+
+f1 = ROOT.TFile("80GeV_first.root")
+tree1 = f1.Get("tree1")
+
+
+
+
+vecPx = ROOT.std.vector(float)()
+vecPy = ROOT.std.vector(float)()
+vecPz = ROOT.std.vector(float)()
+vecPDG = ROOT.std.vector(float)()
+Num = ROOT.std.vector(int)()
+t.Branch("vecPx", vecPx)
+t.Branch("vecPy", vecPy)
+t.Branch("vecPz", vecPz)
+t.Branch("vecPDG", vecPDG)
+t.Branch("N", Num)
+print(tree1.GetEntries())
+for x in range(tree1.GetEntries()):
+ tree1.GetEntry(x)
+ px= list(tree1.FASER_Px)
+ py= list(tree1.FASER_Py)
+ pz=list(tree1.FASER_Pz)
+ pdg=list(tree1.FASER_PDGid)
+ n = (tree1.n)
+ j = 0
+ vecPx.clear()
+ vecPy.clear()
+ vecPz.clear()
+ vecPDG.clear()
+ Num.clear()
+ for j in range(n):
+ vecPx.push_back(px[j])
+ vecPy.push_back(py[j])
+ vecPz.push_back(pz[j])
+ vecPDG.push_back(pdg[j])
+ Num.push_back(n)
+ t.Fill()
+
+
+f.Write()
+f.Close()
diff --git a/Generators/TestBeamReader/share/TestBeamReader_jobOptionsH2_CMS.py b/Generators/TestBeamReader/share/TestBeamReader_jobOptionsH2_CMS.py
new file mode 100644
index 0000000000000000000000000000000000000000..3656ceefec3cfa79d329341d492e2f9c7e0e0e29
--- /dev/null
+++ b/Generators/TestBeamReader/share/TestBeamReader_jobOptionsH2_CMS.py
@@ -0,0 +1,39 @@
+#
+# Usage: athena.py -c'INPUT=["/eos/project/f/faser-preshower/simulations/H2_beamline/1M_events/FASER_e-_v46_100GeV.root"]; OUTPUT="100GeV.H2.EVNT.pool.root"' TestBeamReader/TestBeamReader_jobOptionsH2_CMS.py
+# INPUT and OUTPUT are mandatory (INPUT can be a list, as shown above)
+# EVTMAX can be omitted; then all input files are read to the end
+#
+
+if not 'INPUT' in dir():
+ print("Missing INPUT parameter")
+ exit()
+
+if not isinstance(INPUT, (list,tuple)):
+ INPUT = [INPUT]
+ pass
+
+if not 'OUTPUT' in dir():
+ print("Missing OUTPUT parameter")
+ exit()
+
+import AthenaRootComps.ReadAthenaRoot
+
+svcMgr.EventSelector.InputCollections = INPUT
+svcMgr.EventSelector.TupleName = "t"
+
+from TestBeamReader.TestBeamReaderAlgH2_CMS import TestBeamReader
+
+from AthenaCommon.AlgSequence import AlgSequence
+job = AlgSequence()
+job += TestBeamReader()
+
+from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream
+ostream = AthenaPoolOutputStream( "StreamEVGEN" , OUTPUT, noTag=True )
+ostream.ItemList.remove("EventInfo#*")
+ostream.ItemList += [ "EventInfo#McEventInfo",
+ "McEventCollection#*" ]
+print(ostream.ItemList)
+if not 'EVTMAX' in dir():
+ EVTMAX = -1
+
+theApp.EvtMax = EVTMAX
diff --git a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
index 607c58d892d201731098f64829a82d556f1a302a..dd8faf032f48b2c3d78453b84494e1903d64367d 100644
--- a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
+++ b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
@@ -7,7 +7,7 @@ atlas_add_component(
src/NtupleDumperAlg.h
src/NtupleDumperAlg.cxx
src/component/NtupleDumper_entries.cxx
- LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserCalorimeter xAODFaserTrigger xAODFaserLHC ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserCalorimeter xAODFaserTrigger xAODFaserLHC ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib FaserActsmanVertexingLib
PRIVATE_LINK_LIBRARIES nlohmann_json::nlohmann_json
)
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 6e4c2aae531c6feee89df6102a90d118f6276dff..0e88d1e46f6d84ed87bda1b669c9f4ca35aa3ea6 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -17,7 +17,9 @@
#include <cmath>
#include <TH1F.h>
#include <numeric>
+#include <limits>
+constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
using namespace Acts::UnitLiterals;
NtupleDumperAlg::NtupleDumperAlg(const std::string &name,
@@ -119,6 +121,7 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(m_trackingGeometryTool.retrieve());
ATH_CHECK(m_trackTruthMatchingTool.retrieve());
ATH_CHECK(m_fiducialParticleTool.retrieve());
+ ATH_CHECK(m_vertexingTool.retrieve());
ATH_CHECK(m_spacePointContainerKey.initialize());
@@ -222,6 +225,9 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("longTracks", &m_longTracks, "longTracks/I");
m_tree->Branch("Track_Chi2", &m_Chi2);
m_tree->Branch("Track_nDoF", &m_DoF);
+ m_tree->Branch("Track_module_eta0", &m_module_eta0);
+ m_tree->Branch("Track_module_phi0", &m_module_phi0);
+ m_tree->Branch("Track_hitSet", &m_hitSet);
m_tree->Branch("Track_x0", &m_xup);
m_tree->Branch("Track_y0", &m_yup);
m_tree->Branch("Track_z0", &m_zup);
@@ -289,6 +295,8 @@ StatusCode NtupleDumperAlg::initialize()
//TRUTH
m_tree->Branch("t_pdg", &m_t_pdg);
m_tree->Branch("t_barcode", &m_t_barcode);
+ m_tree->Branch("t_pdg_parent", &m_t_pdg_parent);
+ m_tree->Branch("t_barcode_parent", &m_t_barcode_parent);
m_tree->Branch("t_truthHitRatio", &m_t_truthHitRatio);
m_tree->Branch("t_prodVtx_x", &m_t_prodVtx_x);
@@ -375,6 +383,19 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("truthd1_y", &m_truthd1_y);
m_tree->Branch("truthd1_z", &m_truthd1_z);
+ m_tree->Branch("vertex_x", &m_vertex_x, "vertex_x/D");
+ m_tree->Branch("vertex_y", &m_vertex_y, "vertex_y/D");
+ m_tree->Branch("vertex_z", &m_vertex_z, "vertex_z/D");
+ m_tree->Branch("vertex_chi2", &m_vertex_chi2, "vertex_chi2/D");
+ m_tree->Branch("vertex_px0", &m_vertex_px0, "vertex_px0/D");
+ m_tree->Branch("vertex_py0", &m_vertex_py0, "vertex_py0/D");
+ m_tree->Branch("vertex_pz0", &m_vertex_pz0, "vertex_pz0/D");
+ m_tree->Branch("vertex_p0", &m_vertex_p0, "vertex_p0/D");
+ m_tree->Branch("vertex_px1", &m_vertex_px1, "vertex_px1/D");
+ m_tree->Branch("vertex_py1", &m_vertex_py1, "vertex_py1/D");
+ m_tree->Branch("vertex_pz1", &m_vertex_pz1, "vertex_pz1/D");
+ m_tree->Branch("vertex_p1", &m_vertex_p1, "vertex_p1/D");
+
ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
// Register histograms
@@ -885,6 +906,8 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
ATH_CHECK(trackCollection.isValid());
+ // create list of good tracks in fiducial region (r < 95mm) and chi2/nDoF < 25, #Layers >= 7, #Hits >= 12
+ std::vector<const Trk::Track*> goodTracks {};
for (const Trk::Track* track : *trackCollection)
{
if (track == nullptr) continue;
@@ -892,16 +915,20 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
std::set<std::pair<int, int>> layerMap;
std::set<int> stationMap;
// Check for hit in the three downstream stations
+ HitSet hitSet {24};
for (auto measurement : *(track->measurementsOnTrack())) {
const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
if (cluster != nullptr) {
Identifier id = cluster->identify();
int station = m_sctHelper->station(id);
int layer = m_sctHelper->layer(id);
+ int side = m_sctHelper->side(id);
stationMap.emplace(station);
layerMap.emplace(station, layer);
+ hitSet.set(23 - (station * 6 + layer * 2 + side));
}
}
+ m_hitSet.push_back(hitSet.to_ulong());
if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue;
const Trk::TrackParameters* upstreamParameters = track->trackParameters()->front();
@@ -913,6 +940,12 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_Chi2.push_back(track->fitQuality()->chiSquared());
m_DoF.push_back(track->fitQuality()->numberDoF());
+ const Trk::MeasurementBase *measurement = track->measurementsOnTrack()->front();
+ const Tracker::FaserSCT_ClusterOnTrack* cluster =
+ dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+ Identifier id = cluster->identify();
+ m_module_eta0.push_back(m_sctHelper->eta_module(id));
+ m_module_phi0.push_back(m_sctHelper->phi_module(id));
m_nHit0.push_back(stationMap.count(0));
m_nHit1.push_back(stationMap.count(1));
@@ -951,6 +984,13 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_y_atMaxRadius.push_back(positionAtMaxRadius.y());
m_z_atMaxRadius.push_back(positionAtMaxRadius.z());
+ if ((m_Chi2.back() / m_DoF.back() < m_maxChi2NDOF) &&
+ (m_nLayers.back() >= m_minLayers) &&
+ (m_DoF.back() + 5 >= m_minHits) &&
+ (m_r_atMaxRadius.back() < 95)) {
+ goodTracks.push_back(track);
+ }
+
if (isMC) { // if simulation, store track truth info as well
auto [truthParticle, hitCount] = m_trackTruthMatchingTool->getTruthParticle(track);
if (truthParticle != nullptr) {
@@ -958,6 +998,13 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
truthParticleCount[truthParticle->barcode()] += 1;
m_t_pdg.push_back(truthParticle->pdgId());
m_t_barcode.push_back(truthParticle->barcode());
+ if (truthParticle->nParents() > 0){
+ m_t_pdg_parent.push_back(truthParticle->parent(0)->pdgId());
+ m_t_barcode_parent.push_back(truthParticle->parent(0)->barcode());
+ } else {
+ m_t_pdg_parent.push_back(0);
+ m_t_barcode_parent.push_back(-1);
+ }
// the track fit eats up 5 degrees of freedom, thus the number of hits on track is m_DoF + 5
m_t_truthHitRatio.push_back(hitCount / (m_DoF.back() + 5));
m_isFiducial.push_back(m_fiducialParticleTool->isFiducial(truthParticle->barcode()));
@@ -1152,6 +1199,41 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_longTracks++;
}
+ if (m_runVertexing && goodTracks.size() >= 2) {
+ // sort tracks my momentum and reconstruct vertex
+ std::sort(goodTracks.begin(), goodTracks.end(), [](const Trk::Track *lhs, const Trk::Track *rhs){
+ return lhs->trackParameters()->front()->momentum().z() < rhs->trackParameters()->front()->momentum().z();
+ });
+ std::vector<const Trk::Track*> lowMomentumTracks {goodTracks[0], goodTracks[1]};
+
+ std::optional<FaserTracking::POCA> vertex = m_vertexingTool->getVertex(lowMomentumTracks);
+ if (vertex) {
+ Eigen::Vector3d position = vertex->position;
+ m_vertex_x = position.x();
+ m_vertex_y = position.y();
+ m_vertex_z = position.z();
+ m_vertex_chi2 = vertex->chi2;
+
+ // get track parameters of two tracks with lowest momentum at vertex
+ auto vertexTrackParameters0 = m_vertexingTool->extrapolateTrack(goodTracks[0], m_vertex_z);
+ if (vertexTrackParameters0 != nullptr) {
+ // convert momentum from GeV to MeV
+ m_vertex_px0 = vertexTrackParameters0->momentum().x() * 1e3;
+ m_vertex_py0 = vertexTrackParameters0->momentum().y() * 1e3;
+ m_vertex_pz0 = vertexTrackParameters0->momentum().z() * 1e3;
+ m_vertex_p0 = std::sqrt(m_vertex_px0*m_vertex_px0 + m_vertex_py0*m_vertex_py0 + m_vertex_pz0*m_vertex_pz0);
+ }
+ auto vertexTrackParameters1 = m_vertexingTool->extrapolateTrack(goodTracks[1], m_vertex_z);
+ if (vertexTrackParameters1 != nullptr) {
+ // convert momentum from GeV to MeV
+ m_vertex_px1 = vertexTrackParameters1->momentum().x() * 1e3;
+ m_vertex_py1 = vertexTrackParameters1->momentum().y() * 1e3;
+ m_vertex_pz1 = vertexTrackParameters1->momentum().z() * 1e3;
+ m_vertex_p1 = std::sqrt(m_vertex_px1*m_vertex_px1 + m_vertex_py1*m_vertex_py1 + m_vertex_pz1*m_vertex_pz1);
+ }
+ }
+ }
+
if (!isMC) {
if (m_doTrackFilter) { // filter events: colliding bunches have at least one long track, non-colliding bunches have at least one long track or one calo module with raw_peak > 3 mV
if (m_longTracks == 0 && abs(m_distanceToCollidingBCID) <= 1) {
@@ -1349,6 +1431,9 @@ NtupleDumperAlg::clearTree() const
m_Chi2.clear();
m_DoF.clear();
+ m_module_eta0.clear();
+ m_module_phi0.clear();
+ m_hitSet.clear();
m_charge.clear();
m_nLayers.clear();
m_longTracks = 0;
@@ -1400,6 +1485,8 @@ NtupleDumperAlg::clearTree() const
m_t_pdg.clear();
m_t_barcode.clear();
+ m_t_pdg_parent.clear();
+ m_t_barcode_parent.clear();
m_t_truthHitRatio.clear();
m_t_prodVtx_x.clear();
m_t_prodVtx_y.clear();
@@ -1466,11 +1553,25 @@ NtupleDumperAlg::clearTree() const
m_truthd1_y.clear();
m_truthd1_z.clear();
+ m_vertex_x = NaN;
+ m_vertex_y = NaN;
+ m_vertex_z = NaN;
+ m_vertex_chi2 = NaN;
+ m_vertex_px0 = NaN;
+ m_vertex_py0 = NaN;
+ m_vertex_pz0 = NaN;
+ m_vertex_p0 = NaN;
+ m_vertex_px1 = NaN;
+ m_vertex_py1 = NaN;
+ m_vertex_pz1 = NaN;
+ m_vertex_p1 = NaN;
}
void NtupleDumperAlg::clearTrackTruth() const {
m_t_pdg.push_back(0);
m_t_barcode.push_back(-1);
+ m_t_pdg_parent.push_back(0);
+ m_t_barcode_parent.push_back(-1);
m_t_truthHitRatio.push_back(0);
m_t_prodVtx_x.push_back(999999);
m_t_prodVtx_y.push_back(999999);
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 67f3bacf8c8f4a6ef134741637d093fdea60c24b..a741bcbe55097ba754a256fe3804f0e6cc170db9 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -23,10 +23,14 @@
#include "TrackerSimEvent/FaserSiHitCollection.h"
#include "xAODEventInfo/EventInfo.h"
#include "StoreGate/ReadDecorHandle.h"
+#include "FaserActsVertexing/IVertexingTool.h"
+#include <boost/dynamic_bitset.hpp>
#include <vector>
#include <nlohmann/json.hpp>
+using HitSet = boost::dynamic_bitset<>;
+
class TTree;
class TH1;
class FaserSCT_ID;
@@ -93,6 +97,7 @@ private:
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
ToolHandle<ITrackTruthMatchingTool> m_trackTruthMatchingTool {this, "TrackTruthMatchingTool", "TrackTruthMatchingTool"};
ToolHandle<IFiducialParticleTool> m_fiducialParticleTool {this, "FiducialParticleTool", "FiducialParticleTool"};
+ ToolHandle<FaserTracking::IVertexingTool> m_vertexingTool { this, "VertexingTool", "FaserTracking::PointOfClosestApproachSearchTool"};
const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
@@ -124,6 +129,12 @@ private:
BooleanProperty m_applyGoodRunsList {this, "ApplyGoodRunsList", false, "Only write out events passing GRL (data only)"};
StringProperty m_goodRunsList {this, "GoodRunsList", "", "GoodRunsList in json format"};
+ // track quality cuts
+ UnsignedIntegerProperty m_minLayers{this, "minLayers", 7, "Miminimum number of layers of a track."};
+ UnsignedIntegerProperty m_minHits{this, "minHits", 12, "Miminimum number of hits of a track."};
+ DoubleProperty m_maxChi2NDOF{this, "maxChi2NDOF", 25, "Maximum chi2 per degree of freedom."};
+ BooleanProperty m_runVertexing{ this, "RunVertexing", true, "Run the vertexing, defaults to true." };
+
// json object to hold data read from GRL file (or empty if not)
nlohmann::json m_grl;
@@ -217,6 +228,9 @@ private:
mutable int m_longTracks;
mutable int m_propagationError;
+ mutable std::vector<int> m_module_eta0;
+ mutable std::vector<int> m_module_phi0;
+ mutable std::vector<unsigned long> m_hitSet;
mutable std::vector<double> m_Chi2;
mutable std::vector<double> m_DoF;
mutable std::vector<double> m_xup;
@@ -274,6 +288,8 @@ private:
mutable std::vector<int> m_t_pdg; // pdg code of the truth matched particle
mutable std::vector<int> m_t_barcode; // barcode of the truth matched particle
+ mutable std::vector<int> m_t_pdg_parent; // pdg code of the parent of the truth matched particle
+ mutable std::vector<int> m_t_barcode_parent; // barcode of the parent of the truth matched particle
mutable std::vector<double> m_t_truthHitRatio; // ratio of hits on track matched to the truth particle over all hits on track
mutable std::vector<double> m_t_prodVtx_x; // x component of the production vertex in mm
mutable std::vector<double> m_t_prodVtx_y; // y component of the production vertex in mm
@@ -355,6 +371,19 @@ private:
mutable int m_eventsPassed = 0;
mutable int m_eventsFailedGRL = 0;
+
+ mutable double m_vertex_x; // components of reconstructed vertex in mm
+ mutable double m_vertex_y;
+ mutable double m_vertex_z;
+ mutable double m_vertex_chi2; // chi2 value of vertex fit
+ mutable double m_vertex_px0; // components of lowest momentum track extrapolated to vertex in MeV
+ mutable double m_vertex_py0;
+ mutable double m_vertex_pz0;
+ mutable double m_vertex_p0;
+ mutable double m_vertex_px1; // components of second lowest momentum track extrapolated to vertex in MeV
+ mutable double m_vertex_py1;
+ mutable double m_vertex_pz1;
+ mutable double m_vertex_p1;
};
inline const ServiceHandle <ITHistSvc> &NtupleDumperAlg::histSvc() const {
diff --git a/Scintillator/ScintDetDescr/VetoGeoModel/src/VetoStation.cxx b/Scintillator/ScintDetDescr/VetoGeoModel/src/VetoStation.cxx
index 7c0723d54966f8c65a8a26f322fa2c460af89c9b..05e4a8b0bcb628a699b917cfecf9612d7910aeae 100644
--- a/Scintillator/ScintDetDescr/VetoGeoModel/src/VetoStation.cxx
+++ b/Scintillator/ScintDetDescr/VetoGeoModel/src/VetoStation.cxx
@@ -102,7 +102,7 @@ VetoStation::build(VetoIdentifier id)
id.setPlate(iPlate);
GeoAlignableTransform* transform = new GeoAlignableTransform(GeoTrf::Translate3D(0.0,
0.0,
- (plateThickness - activeDepth)/2 + iPlate * m_platePitch));
+ plateCenterZ));
station->add(transform);
GeoVPhysVol* platePV = m_plate->build(id);
station->add(platePV);
diff --git a/Simulation/G4Faser/G4FaserAlg/CMakeLists.txt b/Simulation/G4Faser/G4FaserAlg/CMakeLists.txt
index 979897a94c4b476433d0ef3fe2ceeb31105f9ef5..ec2babac5f7c92340ffff8b39f318e7e7013fde2 100644
--- a/Simulation/G4Faser/G4FaserAlg/CMakeLists.txt
+++ b/Simulation/G4Faser/G4FaserAlg/CMakeLists.txt
@@ -31,21 +31,32 @@ atlas_add_test( G4FaserAlgConfig_TestFaser
PROPERTIES TIMEOUT 300
PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
-atlas_add_test( G4FaserAlgConfig_TestFaserNu
- SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/G4FaserAlgConfigNew_Test.py GeoModel.FaserVersion="'FASERNU-03'" IOVDb.GlobalTag="'OFLCOND-FASER-02'" Output.HITSFileName='faserNu.HITS.pool.root'
+atlas_add_test( G4FaserAlgConfig_TestFaserNu03
+ SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/G4FaserAlgConfigNew_Test.py GeoModel.FaserVersion="'FASERNU-03'" IOVDb.GlobalTag="'OFLCOND-FASER-02'" Output.HITSFileName='faserNu03.HITS.pool.root'
PROPERTIES TIMEOUT 300
PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
-atlas_add_test( G4FaserAlgConfig_TestFaserNu_NewField
- SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/G4FaserAlgConfigNew_Test.py GeoModel.FaserVersion="'FASERNU-03'" IOVDb.GlobalTag="'OFLCOND-FASER-03'" Output.HITSFileName='faserNu.HITS.pool.root'
+atlas_add_test( G4FaserAlgConfig_TestFaserNu03_NewField
+ SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/G4FaserAlgConfigNew_Test.py GeoModel.FaserVersion="'FASERNU-03'" IOVDb.GlobalTag="'OFLCOND-FASER-03'" Output.HITSFileName='faserNu03_NewField.HITS.pool.root'
PROPERTIES TIMEOUT 300
PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
-atlas_add_test( G4FaserAlgConfig_TestTestbeam
- SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/G4FaserAlgConfigNew_Test.py GeoModel.FaserVersion="'FASER-TB00'" IOVDb.GlobalTag="'OFLCOND-FASER-TB00'" Output.HITSFileName='tb.HITS.pool.root'
+atlas_add_test( G4FaserAlgConfig_TestFaserNu04
+ SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/G4FaserAlgConfigNew_Test.py GeoModel.FaserVersion="'FASERNU-04'" IOVDb.GlobalTag="'OFLCOND-FASER-03'" Output.HITSFileName='faserNu04.HITS.pool.root'
PROPERTIES TIMEOUT 300
PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
+atlas_add_test( G4FaserAlgConfig_TestTestbeam00
+ SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/G4FaserAlgConfigNew_Test.py GeoModel.FaserVersion="'FASER-TB00'" IOVDb.GlobalTag="'OFLCOND-FASER-TB00'" Output.HITSFileName='tb00.HITS.pool.root'
+ PROPERTIES TIMEOUT 300
+ PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
+
+atlas_add_test( G4FaserAlgConfig_TestTestbeam01
+ SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/G4FaserAlgConfigNew_Test.py GeoModel.FaserVersion="'FASER-TB01'" IOVDb.GlobalTag="'OFLCOND-FASER-TB00'" Output.HITSFileName='tb01.HITS.pool.root'
+ PROPERTIES TIMEOUT 300
+ PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
+
+
# Install files from the package:
atlas_install_python_modules( python/*.py )
atlas_install_scripts( test/*.py )
diff --git a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
index df6abb7ab1a0fec74ad0da4b2a4fc0cad8ca84a8..fb07e943ab87584ec8990cf7e1dbdb2d75359795 100755
--- a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
+++ b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
@@ -48,7 +48,7 @@ if __name__ == '__main__':
ConfigFlags.addFlag("Sim.Beam.xshift", 0) # Potential beam shift
ConfigFlags.addFlag("Sim.Beam.yshift", 0)
- ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Geometry set-up
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-04" # Geometry set-up
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Conditions set-up
ConfigFlags.addFlag("Input.InitialTimeStamp", 0) # To avoid autoconfig
ConfigFlags.GeoModel.Align.Dynamic = False
@@ -172,6 +172,11 @@ if __name__ == '__main__':
#cfg.getEventAlgo("OutputStreamHITS").ItemList += ["McEventCollection#BeamTruthEvent_ATLASCoord"]
#
+# Uncomment to check volumes for overlap - will cause CTest to fail due to overwriting file
+#
+# from G4DebuggingTools.G4DebuggingToolsConfigNew import VolumeDebugger
+# cfg.merge(VolumeDebugger(ConfigFlags, name="G4UA::UserActionSvc", TargetVolume="", Verbose=True))
+#
# Dump config
#
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/Tracker/TrackerDetDescr/DipoleGeoModel/src/DipoleFactory.cxx b/Tracker/TrackerDetDescr/DipoleGeoModel/src/DipoleFactory.cxx
index 6129c49ccd609137bd92e9b22f0eb2a161a7cf67..ccfda716028494d56d7c2cac629c42113252b1cd 100644
--- a/Tracker/TrackerDetDescr/DipoleGeoModel/src/DipoleFactory.cxx
+++ b/Tracker/TrackerDetDescr/DipoleGeoModel/src/DipoleFactory.cxx
@@ -93,13 +93,16 @@ void DipoleFactory::create(GeoPhysVol *world )
double upstreamThickness = wrappingParameters->thickness1();
double downstreamThickness = wrappingParameters->thickness2();
double wrappingRadius = wrappingParameters->radius();
+ msg(MSG::ALWAYS) << "Upstream wrapping thickness: " << upstreamThickness << "; Downstream wrapping thickness: " << downstreamThickness << endmsg;
+ msg(MSG::ALWAYS) << "Upstream wrapping position: " << -(parameters->longLength() + upstreamThickness)/2 << "; Downstream wrapping position: " << (parameters->longLength() + downstreamThickness)/2 << endmsg;
GeoTrf::Transform3D dipoleTransform = parameters->partTransform("Dipole");
const GeoTube* shortShape = new GeoTube(parameters->innerRadius(), parameters->outerRadius(), parameters->shortLength()/2);
GeoLogVol* shortLog = new GeoLogVol("ShortDipole", shortShape, NdFeB);
const GeoTube* longShape = new GeoTube(parameters->innerRadius(), parameters->outerRadius(), parameters->longLength()/2);
GeoLogVol* longLog = new GeoLogVol("LongDipole", longShape, NdFeB);
- const GeoTube* wrappedShape = new GeoTube(parameters->innerRadius(), std::max(parameters->outerRadius(), wrappingRadius), parameters->longLength()/2 + std::max(upstreamThickness, downstreamThickness));
+ // const GeoTube* wrappedShape = new GeoTube(parameters->innerRadius(), std::max(parameters->outerRadius(), wrappingRadius), parameters->longLength()/2 + std::max(upstreamThickness, downstreamThickness));
+ const GeoTube* wrappedShape = new GeoTube(0.0, std::max(parameters->outerRadius(), wrappingRadius), parameters->longLength()/2 + std::max(upstreamThickness, downstreamThickness));
GeoLogVol* wrappedLog = new GeoLogVol("WrappedLongDipole", wrappedShape, air);
GeoPhysVol* wrappedPV = new GeoPhysVol(wrappedLog);
GeoPhysVol* longPV = new GeoPhysVol(longLog);
@@ -113,7 +116,7 @@ void DipoleFactory::create(GeoPhysVol *world )
const GeoTube* downstreamShape = new GeoTube(0.0, wrappingRadius, downstreamThickness/2);
GeoLogVol* downstreamLog = new GeoLogVol("DownstreamWrapping", downstreamShape, downstreamMaterial);
GeoPhysVol* downstreamPV = new GeoPhysVol(downstreamLog);
- wrappedPV->add(new GeoTransform(GeoTrf::Translate3D(0.0, 0.0, (parameters->longLength() + upstreamThickness)/2)));
+ wrappedPV->add(new GeoTransform(GeoTrf::Translate3D(0.0, 0.0, (parameters->longLength() + downstreamThickness)/2)));
wrappedPV->add(downstreamPV);
std::vector<std::string> partNames {"UpstreamDipole", "CentralDipole", "DownstreamDipole"};
diff --git a/Tracking/Acts/FaserActsGeometry/CMakeLists.txt b/Tracking/Acts/FaserActsGeometry/CMakeLists.txt
index e15adcec74c6efaee85c11e388517683262c6522..58c40ebec78ecb3ab4585b9dc744392c4aa19042 100755
--- a/Tracking/Acts/FaserActsGeometry/CMakeLists.txt
+++ b/Tracking/Acts/FaserActsGeometry/CMakeLists.txt
@@ -56,6 +56,8 @@ atlas_add_component( FaserActsGeometry
src/FaserActsAlignmentCondAlg.cxx
src/NominalAlignmentCondAlg.cxx
src/FaserActsVolumeMappingTool.cxx
+ src/FaserActsGeometryBoundaryTestAlg.h
+ src/FaserActsGeometryBoundaryTestAlg.cxx
src/components/*.cxx
PUBLIC_HEADERS FaserActsGeometry
INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} ${BOOST_INCLUDE_DIRS}
@@ -84,3 +86,7 @@ atlas_install_headers( FaserActsGeometry )
atlas_install_python_modules( python/*.py )
atlas_install_scripts( test/*.py )
+atlas_add_test( FaserActsGeometryBoundary_test
+ SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/test/FaserActsGeometryBoundary_test.py
+ PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
+ PROPERTIES TIMEOUT 300 )
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
index d31fe8d3c11e804e115aef2cc6dfc08a2354b9e7..14d2699822166993a61b0de47a74c1fc5b139f65 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
@@ -71,7 +71,7 @@ FaserActsExtrapolationTool::initialize()
ATH_MSG_INFO("Initializing ACTS extrapolation");
- m_logger = Acts::getDefaultLogger("ExtrapolationTool", Acts::Logging::INFO);
+ m_logger = Acts::getDefaultLogger("ExtrapolationTool", Acts::Logging::ERROR);
// m_logger = makeActsAthenaLogger(this, "Prop", "FaserActsExtrapTool");
ATH_CHECK( m_trackingGeometryTool.retrieve() );
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsGeometryBoundaryTestAlg.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsGeometryBoundaryTestAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..57c1e3a2fe68bea133d167ea00dd539edfc0f3fe
--- /dev/null
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsGeometryBoundaryTestAlg.cxx
@@ -0,0 +1,97 @@
+/*
+ Copyright (C) 2002-2023 CERN for the benefit of the ATLAS and FASER collaborations
+*/
+
+
+#include "FaserActsGeometryBoundaryTestAlg.h"
+#include "FaserActsGeometry/FaserActsGeometryContext.h"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/Geometry/DetectorElementBase.hpp"
+#include "Acts/Geometry/TrackingVolume.hpp"
+
+
+FaserActsGeometryBoundaryTestAlg::FaserActsGeometryBoundaryTestAlg(const std::string &name, ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator) {}
+
+
+StatusCode FaserActsGeometryBoundaryTestAlg::initialize() {
+ ATH_CHECK(m_trackingGeometryTool.retrieve());
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode FaserActsGeometryBoundaryTestAlg::execute(const EventContext &ctx) const {
+
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
+ = m_trackingGeometryTool->trackingGeometry();
+
+ const Acts::GeometryContext gctx =
+ m_trackingGeometryTool->getGeometryContext(ctx).context();
+
+ // loop over all tracking layer and check if contained detector elements are insie boundaries
+ // return StatusCode::Failure so that test fails, if the edge of any detector element is outside boundaries
+ const Acts::TrackingVolume *trackingVolume = trackingGeometry->highestTrackingVolume();
+ for (auto volume : trackingVolume->confinedVolumes()->arrayObjects()) {
+ for (const auto& layer : volume->confinedLayers()->arrayObjects()) {
+ if (layer->layerType() == Acts::LayerType::active) {
+ // get inner and outer boundaries from approach descriptor:
+ // - innerBoundary is approach suface with minimum z position
+ // - outerBoundary is approach suface with maximum z position
+ Acts::SurfaceVector approachSurfaces = layer->approachDescriptor()->containedSurfaces();
+ const Acts::Surface *innerApproachSurface = *std::min_element(
+ approachSurfaces.begin(), approachSurfaces.end(),
+ [&gctx](const auto &lhs, const auto &rhs) {
+ return lhs->center(gctx).z() < rhs->center(gctx).z();
+ }
+ );
+ const Acts::Surface *outerApproachSurface = *std::max_element(
+ approachSurfaces.begin(), approachSurfaces.end(),
+ [&gctx](const auto &lhs, const auto &rhs) {
+ return lhs->center(gctx).z() < rhs->center(gctx).z();
+ }
+ );
+ double zInnerBoundary = innerApproachSurface->center(gctx).z();
+ double zOuterBoundary = outerApproachSurface->center(gctx).z();
+
+ // loop over surface array and check if all edges are between inner and outer boundary
+ for (const Acts::Surface *surface : layer->surfaceArray()->surfaces()) {
+ auto planeSurface = dynamic_cast<const Acts::PlaneSurface *>(surface);
+ // make sure surface has a associated detector element (there are other active detector elements
+ // e.g. containing material or at the tracking volume boundaries)
+ if (surface->associatedDetectorElement() != nullptr) {
+ auto bounds = dynamic_cast<const Acts::RectangleBounds*>(&surface->bounds());
+ const Acts::Vector2 min = bounds->min();
+ const Acts::Vector2 max = bounds->max();
+ // create dummpy momentum vector: local to global transformation requires momentum vector,
+ // which is ignored for PlaneSurface
+ Acts::Vector3 dummyMomentum {1., 1., 1.};
+ // get global position at all edges of the surface
+ std::vector<Acts::Vector3> edges {};
+ edges.push_back(planeSurface->localToGlobal(gctx, {min.x(), min.y()}, dummyMomentum));
+ edges.push_back(planeSurface->localToGlobal(gctx, {min.x(), max.y()}, dummyMomentum));
+ edges.push_back(planeSurface->localToGlobal(gctx, {max.x(), min.y()}, dummyMomentum));
+ edges.push_back(planeSurface->localToGlobal(gctx, {max.x(), max.y()}, dummyMomentum));
+ for (const Acts::Vector3 &edgePosition : edges) {
+ if ((edgePosition.z() < zInnerBoundary) || (edgePosition.z() > zOuterBoundary)) {
+ std::cout << "?? surface outside boundaries\n";
+ std::cout << "inner Boundary: " << zInnerBoundary << std::endl;
+ std::cout << "outer Boundary: " << zOuterBoundary << std::endl;
+ std::cout << "edge: " << edgePosition.x() << ", " << edgePosition.y() << ", " << edgePosition.z() << std::endl;
+ return StatusCode::FAILURE;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode FaserActsGeometryBoundaryTestAlg::finalize() {
+ return StatusCode::SUCCESS;
+}
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsGeometryBoundaryTestAlg.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsGeometryBoundaryTestAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2c0bf3c73b35ac1d0b9857a8c0aae72cfb3c56c
--- /dev/null
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsGeometryBoundaryTestAlg.h
@@ -0,0 +1,25 @@
+/*
+ Copyright (C) 2002-2023 CERN for the benefit of the ATLAS and FASER collaborations
+*/
+
+#ifndef FASERACTSGEOMETRY_FASERACTSGEOMETRYBOUNDARYTESTALG_H
+#define FASERACTSGEOMETRY_FASERACTSGEOMETRYBOUNDARYTESTALG_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+
+
+class FaserActsGeometryBoundaryTestAlg : public AthReentrantAlgorithm {
+public:
+ FaserActsGeometryBoundaryTestAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual StatusCode initialize() override final;
+ virtual StatusCode execute(const EventContext &ctx) const override final;
+ virtual StatusCode finalize() override final;
+
+private:
+ ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {
+ this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
+};
+
+
+#endif // FASERACTSGEOMETRY_FASERACTSGEOMETRYBOUNDARYTESTALG_H
\ No newline at end of file
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
index c5f6a68ab01b8da95f3d62149d23bdf41684263d..649f2503cb69c93b96b8e454727998656898eccf 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx
@@ -206,7 +206,8 @@ FaserActsLayerBuilder::buildLayers(const Acts::GeometryContext& gctx,
std::shared_ptr<const Acts::ProtoSurfaceMaterial> materialProxy = nullptr;
double layerZ = 0.5 * (pl.min(Acts::binZ) + pl.max(Acts::binZ));
- double layerThickness = (pl.max(Acts::binZ) - pl.min(Acts::binZ));
+ // increase layerThickness by 4 mm so that after alignment shifts all modules are still within the boundaries
+ double layerThickness = (pl.max(Acts::binZ) - pl.min(Acts::binZ)) + 4_mm;
double layerZInner = layerZ - layerThickness/2.;
double layerZOuter = layerZ + layerThickness/2.;
diff --git a/Tracking/Acts/FaserActsGeometry/src/components/FaserActsGeometry_entries.cxx b/Tracking/Acts/FaserActsGeometry/src/components/FaserActsGeometry_entries.cxx
index e7b7fb2fe49208911b0d98c2ba4a6d6b960df89f..a103c7a9466eb66446036cc2978e82dba8889fee 100755
--- a/Tracking/Acts/FaserActsGeometry/src/components/FaserActsGeometry_entries.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/components/FaserActsGeometry_entries.cxx
@@ -16,6 +16,7 @@
#include "../FaserActsMaterialMapping.h"
#include "../FaserActsSurfaceMappingTool.h"
#include "../FaserActsMaterialTrackWriterSvc.h"
+#include "../FaserActsGeometryBoundaryTestAlg.h"
DECLARE_COMPONENT( FaserActsTrackingGeometrySvc )
DECLARE_COMPONENT( FaserActsTrackingGeometryTool )
@@ -30,3 +31,4 @@ DECLARE_COMPONENT( FaserActsMaterialJsonWriterTool )
DECLARE_COMPONENT( FaserActsMaterialTrackWriterSvc )
DECLARE_COMPONENT( FaserActsMaterialMapping )
DECLARE_COMPONENT( FaserActsSurfaceMappingTool )
+DECLARE_COMPONENT( FaserActsGeometryBoundaryTestAlg )
diff --git a/Tracking/Acts/FaserActsGeometry/test/FaserActsGeometryBoundary_test.py b/Tracking/Acts/FaserActsGeometry/test/FaserActsGeometryBoundary_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..b6a841642b17b6ecd412fc55ced498524ece9ca5
--- /dev/null
+++ b/Tracking/Acts/FaserActsGeometry/test/FaserActsGeometryBoundary_test.py
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+"""
+Copyright (C) 2002-2023 CERN for the benefit of the ATLAS and FASER collaboration
+"""
+
+
+def FaserActsGeometryBoundaryTestCfg(flags, name="FaserActsGeometryBoundaryTestAlg", **kwargs):
+
+ from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+ from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+ from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometryToolCfg
+ from AthenaConfiguration.ComponentFactory import CompFactory
+
+ acc = FaserSCT_GeometryCfg(flags)
+ acc.merge(MagneticFieldSvcCfg(flags))
+ result, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(flags)
+ test_alg = CompFactory.FaserActsGeometryBoundaryTestAlg
+ test_alg.TrackingGeometryTool = actsTrackingGeometryTool
+ acc.merge(result)
+ acc.addEventAlgo(test_alg(name, **kwargs))
+ return acc
+
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+
+ Configurable.configurableRun3Behavior = True
+ ConfigFlags.Input.isMC = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-04"
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03"
+ ConfigFlags.Detector.GeometryFaserSCT = True
+ ConfigFlags.lock()
+
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(FaserGeometryCfg(ConfigFlags))
+ acc.merge(FaserActsGeometryBoundaryTestCfg(ConfigFlags))
+
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
+
+ sys.exit(int(acc.run(maxEvents=1).isFailure()))
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..61908a643b1d9f03eaa3b4ebf9c6eb2d34786a6a
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/CMakeLists.txt
@@ -0,0 +1,44 @@
+atlas_subdir(FaserActsVertexing)
+
+find_package( Eigen )
+find_package( Acts COMPONENTS Core )
+
+atlas_add_library( FaserActsmanVertexingLib
+ PUBLIC_HEADERS FaserActsVertexing
+ INCLUDE_DIRS ${EIGEN_INCLUDE_DIRS}
+ LINK_LIBRARIES ${EIGEN_LIBRARIES}
+ ActsCore
+ AthenaBaseComps
+ FaserActsGeometryLib
+ FaserActsKalmanFilterLib
+ MagFieldConditions
+ StoreGateLib
+ TrackerPrepRawData
+ TrackerRIO_OnTrack
+ TrackerSimEvent
+ TrackerSpacePoint
+ TrkTrack
+ xAODTruth
+)
+
+atlas_add_component(
+ FaserActsVertexing
+ FaserActsVertexing/IVertexingTool.h
+ src/PointOfClosestApproachSearchTool.h
+ src/PointOfClosestApproachSearchTool.cxx
+ src/components/FaserActsVertexing_entries.cxx
+ INCLUDE_DIRS ${EIGEN_INCLUDE_DIRS}
+ LINK_LIBRARIES ${EIGEN_LIBRARIES}
+ ActsCore
+ AthenaBaseComps
+ FaserActsGeometryLib
+ FaserActsKalmanFilterLib
+ MagFieldConditions
+ StoreGateLib
+ TrackerPrepRawData
+ TrackerRIO_OnTrack
+ TrackerSimEvent
+ TrackerSpacePoint
+ TrkTrack
+ xAODTruth
+)
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/FaserActsVertexing/IVertexingTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/FaserActsVertexing/IVertexingTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..3c6f93714b6cf07f2428a80c157318399d7b632e
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/FaserActsVertexing/IVertexingTool.h
@@ -0,0 +1,37 @@
+#ifndef FASERACTSVERTEXING_IVERTEXINGTOOL_H
+#define FASERACTSVERTEXING_IVERTEXINGTOOL_H
+
+#include "Acts/EventData/TrackParameters.hpp"
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/IInterface.h"
+#include "TrkTrack/Track.h"
+#include <vector>
+
+namespace FaserTracking {
+
+// Point of closest approach between n tracks. Contains the three-dimensional
+// reconstructed position, the chi2 value of the vertex fit and the tracks that
+// originate from this vertex.
+struct POCA {
+ POCA(Eigen::Vector3d vertex, double chi2, std::vector<const Trk::Track *> tracks)
+ : position(vertex), chi2(chi2), tracks(tracks) {}
+ Eigen::Vector3d position;
+ double chi2;
+ std::vector<const Trk::Track *> tracks;
+};
+
+class IVertexingTool : virtual public IAlgTool {
+public:
+ DeclareInterfaceID(IVertexingTool, 1, 0);
+
+ // Run vertex fit and return the fitted POCA.
+ virtual std::optional<POCA> getVertex(const std::vector<const Trk::Track *> &tracks) const = 0;
+
+ // Get bound track parameters at given z position.
+ virtual std::unique_ptr<const Acts::BoundTrackParameters>
+ extrapolateTrack(const Trk::Track *track, double targetPosition) const = 0;
+};
+
+} // namespace FaserTracking
+
+#endif // FASERACTSVERTEXING_IVERTEXINGTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ea423d703372dbfcc954dc2089b8169282cdc634
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.cxx
@@ -0,0 +1,212 @@
+#include "PointOfClosestApproachSearchTool.h"
+#include "TrkParameters/TrackParameters.h"
+#include <climits>
+
+namespace FaserTracking {
+
+PointOfClosestApproachSearchTool::PointOfClosestApproachSearchTool(const std::string &type,
+ const std::string &name,
+ const IInterface *parent)
+ : base_class(type, name, parent) {}
+
+StatusCode PointOfClosestApproachSearchTool::initialize() {
+ ATH_CHECK(m_trackingGeometryTool.retrieve());
+ ATH_CHECK(m_extrapolationTool.retrieve());
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PointOfClosestApproachSearchTool::finalize() { return StatusCode::SUCCESS; }
+
+std::optional<POCA>
+PointOfClosestApproachSearchTool::getVertex(const std::vector<const Trk::Track *> &tracks) const {
+
+ if (tracks.size() < 2) {
+ ATH_MSG_DEBUG("Found " << tracks.size() << " tracks, but require at least 2 tracks.");
+ return {};
+ } else if (tracks.size() > 3) {
+ ATH_MSG_DEBUG("Found " << tracks.size() << " tracks, but require at most 3 tracks.");
+ return {};
+ }
+
+ clear();
+
+ m_ctx = Gaudi::Hive::currentContext();
+ m_gctx = m_trackingGeometryTool->getNominalGeometryContext().context();
+
+ // convert to Acts::BoundTrackParameters
+ std::vector<Acts::BoundTrackParameters> trackParameters{};
+ for (const Trk::Track *track : tracks) {
+ trackParameters.push_back(getParametersFromTrack(track->trackParameters()->front()));
+ }
+
+ // extrapolate to nSteps z positions betweeen zMin and zMax and calculate center of tracks and
+ // chi2 for each position.
+ size_t nSteps = (m_zMax - m_zMin) / m_stepSize + 1;
+ m_allVertexCandidates.reserve(nSteps);
+ for (int z = m_zMin; z <= m_zMax; z += m_stepSize) {
+ auto vertexCandidate = getVertexCandidate(trackParameters, z);
+ if (vertexCandidate) {
+ m_allVertexCandidates.push_back(*vertexCandidate);
+ }
+ }
+
+ // find vertex candidates with minium (local) chi2 value
+ for (size_t i = 1; i < nSteps - 1; ++i) {
+ if ((m_allVertexCandidates[i - 1].chi2 > m_allVertexCandidates[i].chi2) and
+ (m_allVertexCandidates[i + 1].chi2 > m_allVertexCandidates[i].chi2)) {
+ m_goodVertexCandidates.push_back(m_allVertexCandidates[i]);
+ }
+ }
+
+ // Check if there are one or two local minima and calculate mean in the case of two minima.
+ if (m_goodVertexCandidates.size() == 0) {
+ ATH_MSG_DEBUG("Cannot find minimum. Extrapolation failed. This is likely an error.");
+ return {};
+ } else if (m_goodVertexCandidates.size() == 1) {
+ const VertexCandidate &vx = m_goodVertexCandidates.front();
+ return POCA(vx.position, vx.chi2, tracks);
+ } else if (m_goodVertexCandidates.size() == 2) {
+ const VertexCandidate &vx0 = m_goodVertexCandidates.at(0);
+ const VertexCandidate &vx1 = m_goodVertexCandidates.at(1);
+ return POCA(0.5 * (vx0.position + vx1.position), 0.5 * (vx0.chi2 + vx1.chi2), tracks);
+ } else {
+ ATH_MSG_DEBUG("Cannot find global minimum. Expect a maximum of two local chi2 minima but found "
+ << m_goodVertexCandidates.size() << " minima.");
+ return {};
+ }
+}
+
+// Calculate center of tracks
+Eigen::Vector2d PointOfClosestApproachSearchTool::getCenter(
+ const std::vector<PointOfClosestApproachSearchTool::TrackPosition> &trackPositions) const {
+ Eigen::Vector2d center{0, 0};
+ Eigen::Vector2d sumInverse{0, 0};
+ for (const TrackPosition &value : trackPositions) {
+ center(0) += value.inv(0, 0) * value.pos(0);
+ sumInverse(0) += value.inv(0, 0);
+ center(1) += value.inv(1, 1) * value.pos(1);
+ sumInverse(1) += value.inv(1, 1);
+ }
+ center(0) /= sumInverse(0);
+ center(1) /= sumInverse(1);
+ return center;
+}
+
+// Extrapolate track to given z positon.
+std::unique_ptr<const Acts::BoundTrackParameters>
+PointOfClosestApproachSearchTool::extrapolateTrack(const Trk::Track *track,
+ double targetPosition) const {
+ Acts::BoundTrackParameters startParameters =
+ getParametersFromTrack(track->trackParameters()->front());
+ auto targetSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ Acts::Vector3(0, 0, targetPosition), Acts::Vector3(0, 0, 1));
+ Acts::NavigationDirection navigationDirection =
+ targetPosition > startParameters.referenceSurface().center(m_gctx).z() ? Acts::forward
+ : Acts::backward;
+ return m_extrapolationTool->propagate(m_ctx, startParameters, *targetSurface,
+ navigationDirection);
+}
+
+std::optional<PointOfClosestApproachSearchTool::TrackPosition>
+PointOfClosestApproachSearchTool::extrapolateTrackParameters(
+ const Acts::BoundTrackParameters &startParameters, double z) const {
+ auto targetSurface =
+ Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, z), Acts::Vector3(0, 0, 1));
+ Acts::NavigationDirection navigationDirection =
+ z > startParameters.referenceSurface().center(m_gctx).z() ? Acts::forward : Acts::backward;
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters =
+ m_extrapolationTool->propagate(m_ctx, startParameters, *targetSurface, navigationDirection);
+ if (targetParameters != nullptr && targetParameters->covariance().has_value()) {
+ // extrapolate covariance matrix?
+ // return TrackPosition{targetParameters->position(m_gctx),
+ // startParameters.covariance().value().topLeftCorner(2, 2)};
+ return TrackPosition(targetParameters->position(m_gctx),
+ targetParameters->covariance().value().topLeftCorner(2, 2));
+ } else {
+ ATH_MSG_DEBUG("Extrapolation failed.");
+ return std::nullopt;
+ }
+}
+
+std::optional<PointOfClosestApproachSearchTool::VertexCandidate>
+PointOfClosestApproachSearchTool::getVertexCandidate(
+ std::vector<Acts::BoundTrackParameters> &trackParameters, double z) const {
+ // extrapolate tracks to given z position
+ std::vector<TrackPosition> trackPositions{};
+ for (const Acts::BoundTrackParameters &startParameters : trackParameters) {
+ auto tp = extrapolateTrackParameters(startParameters, z);
+ if (tp) {
+ trackPositions.push_back(*tp);
+ }
+ }
+ if (trackPositions.size() < 2) {
+ ATH_MSG_DEBUG(
+ "Extrapolation failed for too many tracks. Require atleast 2 extrapolated positions.");
+ return std::nullopt;
+ }
+
+ // calculate center weighted by inverse covariance matrix
+ Eigen::Vector2d trackCenter = getCenter(trackPositions);
+
+ // either use weighted x and y position or only y position to calculate chi2 value
+ double vertexChi2 = chi2Y(trackPositions, trackCenter);
+
+ return PointOfClosestApproachSearchTool::VertexCandidate(z, vertexChi2, trackCenter);
+}
+
+double PointOfClosestApproachSearchTool::chi2(const std::vector<TrackPosition> &trackPositions,
+ const Eigen::Vector2d &trackCenter) const {
+ double chi2 = 0;
+ for (const TrackPosition &value : trackPositions) {
+ Eigen::Vector2d delta = trackCenter - value.pos;
+ chi2 += delta.dot(value.inv.diagonal().cwiseProduct(delta));
+ }
+ return chi2;
+}
+
+double PointOfClosestApproachSearchTool::chi2Y(const std::vector<TrackPosition> &trackPositions,
+ const Eigen::Vector2d &trackCenter) const {
+ double chi2 = 0;
+ for (const TrackPosition &value : trackPositions) {
+ Eigen::Vector2d delta = trackCenter - value.pos;
+ chi2 += delta.y() * value.inv(1, 1) * delta.y();
+ }
+ return chi2;
+}
+
+Acts::BoundTrackParameters PointOfClosestApproachSearchTool::getParametersFromTrack(
+ const Trk::TrackParameters *trackParameters) const {
+ using namespace Acts::UnitLiterals;
+
+ Acts::GeometryContext geometryContext =
+ m_trackingGeometryTool->getNominalGeometryContext().context();
+ Acts::AngleAxis3 rotZ(M_PI / 2, Acts::Vector3(0., 0., 1.));
+ Acts::Translation3 translation{0., 0., trackParameters->associatedSurface().center().z()};
+ auto transform = Acts::Transform3(translation * rotZ);
+ auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(transform);
+
+ auto bound = Acts::detail::transformFreeToBoundParameters(
+ trackParameters->position(), 0, trackParameters->momentum(),
+ trackParameters->charge() / trackParameters->momentum().mag() / 1_MeV, *surface,
+ geometryContext);
+ if (!bound.ok()) {
+ ATH_MSG_WARNING("FreeToBound parameter transformation failed");
+ }
+ // convert covariance matrix to GeV
+ Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Identity();
+ cov.topLeftCorner(5, 5) = *trackParameters->covariance();
+ for (int i = 0; i < cov.rows(); i++) {
+ cov(i, 4) = cov(i, 4) / 1_MeV;
+ }
+ for (int i = 0; i < cov.cols(); i++) {
+ cov(4, i) = cov(4, i) / 1_MeV;
+ }
+ return Acts::BoundTrackParameters{surface, bound.value(), trackParameters->charge(), cov};
+}
+
+void PointOfClosestApproachSearchTool::clear() const {
+ m_allVertexCandidates.clear();
+ m_goodVertexCandidates.clear();
+}
+
+} // namespace FaserTracking
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..9ed794948483ee01fcfb13203971ba7a897657da
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.h
@@ -0,0 +1,84 @@
+#ifndef FASERACTSVERTEXING_POINTOFCLOSESTAPPROACHSEARCHTOOL_H
+#define FASERACTSVERTEXING_POINTOFCLOSESTAPPROACHSEARCHTOOL_H
+
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+#include "FaserActsVertexing/IVertexingTool.h"
+#include "TrkTrack/Track.h"
+
+namespace FaserTracking {
+
+class PointOfClosestApproachSearchTool : public extends<AthAlgTool, IVertexingTool> {
+public:
+ PointOfClosestApproachSearchTool(const std::string &type, const std::string &name,
+ const IInterface *parent);
+ virtual ~PointOfClosestApproachSearchTool() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+
+ std::optional<POCA> getVertex(const std::vector<const Trk::Track *> &tracks) const override;
+
+ std::unique_ptr<const Acts::BoundTrackParameters>
+ extrapolateTrack(const Trk::Track *track, double targetPosition) const override;
+
+private:
+ struct VertexCandidate {
+ VertexCandidate(double z, double chi2, const Eigen::Vector2d ¢er) : chi2(chi2) {
+ position = Eigen::Vector3d(center.x(), center.y(), z);
+ }
+ double chi2;
+ Eigen::Vector3d position;
+ };
+
+ struct TrackPosition {
+ TrackPosition(const Acts::Vector3 &position, const Acts::BoundSymMatrix &covariance) {
+ pos << position.x(), position.y();
+ cov << covariance(1, 1), 0, 0, covariance(0, 0);
+ inv = cov.inverse();
+ }
+ Eigen::Vector2d pos;
+ Eigen::Matrix2d cov;
+ Eigen::Matrix2d inv;
+ };
+
+ ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool{
+ this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool{
+ this, "ExtrapolationTool", "FaserActsExtrapolationTool"};
+ DoubleProperty m_stepSize{this, "StepSize", 50};
+ DoubleProperty m_zMin{this, "ZMin", -2000};
+ DoubleProperty m_zMax{this, "ZMax", 500};
+ StringProperty m_debugLevel{this, "DebugLevel", "INFO"};
+
+ void clear() const;
+
+ std::optional<TrackPosition>
+ extrapolateTrackParameters(const Acts::BoundTrackParameters &startParameters, double z) const;
+
+ Eigen::Vector2d getCenter(const std::vector<TrackPosition> &trackPositions) const;
+
+ Acts::BoundTrackParameters
+ getParametersFromTrack(const Trk::TrackParameters *trackParameters) const;
+
+ double chi2(const std::vector<TrackPosition> &trackPositions,
+ const Eigen::Vector2d &trackCenter) const;
+
+ double chi2Y(const std::vector<TrackPosition> &trackPositions,
+ const Eigen::Vector2d &trackCenter) const;
+
+ std::optional<VertexCandidate>
+ getVertexCandidate(std::vector<Acts::BoundTrackParameters> &trackParameters, double z) const;
+
+ mutable EventContext m_ctx;
+ mutable Acts::GeometryContext m_gctx;
+
+ mutable std::vector<VertexCandidate> m_allVertexCandidates;
+ mutable std::vector<VertexCandidate> m_goodVertexCandidates;
+};
+
+} // namespace FaserTracking
+
+#endif /* FASERACTSVERTEXING_POINTOFCLOSESTAPPROACHSEARCHTOOL_H */
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/components/FaserActsVertexing_entries.cxx b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/components/FaserActsVertexing_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4f27a5a3d95ba97b48d36ac0b43c634c1d99d9ae
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/components/FaserActsVertexing_entries.cxx
@@ -0,0 +1,3 @@
+#include "../PointOfClosestApproachSearchTool.h"
+
+DECLARE_COMPONENT(FaserTracking::PointOfClosestApproachSearchTool)
\ No newline at end of file