Merge branch 'missing_hits_fix' into 'master'

Fix missing hits

See merge request !362

.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

PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx

@@ -225,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);
@@ -912,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();
@@ -933,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));
@@ -1418,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;

PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h

@@ -24,10 +24,13 @@
 #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;
@@ -225,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;

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 )

Tracking/Acts/FaserActsGeometry/src/FaserActsGeometryBoundaryTestAlg.cxx

+/*
+  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;
+}

Tracking/Acts/FaserActsGeometry/src/FaserActsGeometryBoundaryTestAlg.h

+/*
+  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

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.;
 

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 )

Tracking/Acts/FaserActsGeometry/test/FaserActsGeometryBoundary_test.py

+#!/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()))