Merge branch 'new_geometry' of github.com:pkurash/AliceO2 into new_geometry

merge

Author: pkurash

Common/DCAFitter/DCAFitterN_derivation.md

@@ -479,9 +479,11 @@ I_XX -> I_XX + 1/sigma_X,prior^2 .
 In code this is implemented as
 
 ```cpp
-constexpr float XRegErrFactor = 10.f;
-const float sigmaX2 = C_YY * XRegErrFactor;
-sxx += 1.f / sigmaX2;
+static constexpr float XRegErrFactor = 10.f;
+...
+if (xRegErrFactor > 0.f) {
+  sxx += 1.f / (cyy * xRegErrFactor);
+}
 ```
 
 This is different from multiplying `I_XX` by a number below one.  A reduction of
@@ -491,7 +493,15 @@ not alter `I_YY`, `I_YZ`, or `I_ZZ`.
 
 The regularization should remain weak.  It is a numerical stabilizer for badly
 conditioned geometries, not an additional detector measurement of the local
-track `X` coordinate.
+track `X` coordinate.  For this reason it is applied **only** where an invertible
+single-track contribution is actually required, i.e. for the `I_i` entering the
+chi2 minimization (`mTrcEInv`, hence `calcInverseWeight()`, `calcPCACoefs()`,
+`calcChi2()` and the Newton Hessian).  `calcPCACovMatrix()` rebuilds the `I_i`
+with `TrackCovI::XRegNone`: there the prior is not needed (the sum over prongs is
+inverted, not the individual terms, and a genuinely ill-conditioned sum is
+detected and replaced by a loose dummy covariance), and including it would make
+the reported longitudinal vertex error follow the dummy `XRegErrFactor * C_YY`
+instead of the track slopes.
 
 ### `calcPCACovMatrix()`
 

Common/DCAFitter/include/DCAFitter/DCAFitterN.h

@@ -171,6 +171,7 @@ inline void printStat(const FitStatusArray& a)
 
 BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
 {
+  constexpr bool oldMode = false; // if true, use the old mode of DCAFitterN, which is less correct but faster
   constexpr int NTest = 10000;
   o2::utils::TreeStreamRedirector outStream("dcafitterNTest.root");
 
@@ -196,6 +197,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
     std::memset(fitstat.data(), 0, sizeof(fitstat));
 
     o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
+    ft.setOldMode(oldMode);          // use the old mode of DCAFitterN
     ft.setBz(bz);
     ft.setPropagateToPCA(true);  // After finding the vertex, propagate tracks to the DCA. This is default anyway
     ft.setMaxR(200);             // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
@@ -280,6 +282,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
     std::memset(fitstat.data(), 0, sizeof(fitstat));
 
     o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
+    ft.setOldMode(oldMode);          // use the old mode of DCAFitterN
     ft.setBz(bz);
     ft.setPropagateToPCA(true);  // After finding the vertex, propagate tracks to the DCA. This is default anyway
     ft.setMaxR(200);             // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
@@ -366,6 +369,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
     std::memset(fitstat.data(), 0, sizeof(fitstat));
 
     o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
+    ft.setOldMode(oldMode);          // use the old mode of DCAFitterN
     ft.setBz(bz);
     ft.setPropagateToPCA(true);  // After finding the vertex, propagate tracks to the DCA. This is default anyway
     ft.setMaxR(200);             // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
@@ -451,6 +455,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
     std::memset(fitstat.data(), 0, sizeof(fitstat));
 
     o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
+    ft.setOldMode(oldMode);          // use the old mode of DCAFitterN
     ft.setBz(bz);
     ft.setPropagateToPCA(true);  // After finding the vertex, propagate tracks to the DCA. This is default anyway
     ft.setMaxR(200);             // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
@@ -535,6 +540,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
     std::memset(fitstat.data(), 0, sizeof(fitstat));
 
     o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter
+    ft.setOldMode(oldMode);          // use the old mode of DCAFitterN
     ft.setBz(bz);
     ft.setPropagateToPCA(true);  // After finding the vertex, propagate tracks to the DCA. This is default anyway
     ft.setMaxR(200);             // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway

Common/DCAFitter/test/testDCAFitterN.cxx

@@ -34,6 +34,19 @@ class GlobalFwdTrack : public o2::track::TrackParCovFwd, public o2::dataformats:
   GlobalFwdTrack(o2::track::TrackParCovFwd const& t) { *this = t; }
   ~GlobalFwdTrack() = default;
 
+  GlobalFwdTrack& operator=(const TrackParCovFwd& rhs)
+  {
+    o2::track::TrackParCovFwd::operator=(rhs);
+    return *this;
+  }
+
+  GlobalFwdTrack& operator=(const GlobalFwdTrack& rhs)
+  {
+    o2::track::TrackParCovFwd::operator=(rhs);
+    o2::dataformats::MatchInfoFwd::operator=(rhs);
+    return *this;
+  }
+
   SMatrix5 computeResiduals2Cov(const o2::track::TrackParCovFwd& t) const
   {
     SMatrix5 Residuals2Cov;

DataFormats/Reconstruction/include/ReconstructionDataFormats/GlobalFwdTrack.h

@@ -37,7 +37,6 @@ o2_add_library(DetectorsBase
                        src/GlobalParams.cxx
                        src/O2Tessellated.cxx
                        src/TGeoGeometryUtils.cxx
-                       src/CADGeometryUtils.cxx
                PUBLIC_LINK_LIBRARIES FairRoot::Base
                                      O2::CommonUtils
                                      O2::DetectorsCommonDataFormats
@@ -90,6 +89,13 @@ o2_add_test(
   PUBLIC_LINK_LIBRARIES O2::DetectorsBase
   LABELS detectorsbase)
 
+o2_add_test(
+  O2Tessellated
+  SOURCES test/testO2Tessellated.cxx
+  COMPONENT_NAME DetectorsBase
+  PUBLIC_LINK_LIBRARIES O2::DetectorsBase
+  LABELS detectorsbase)
+
 if(BUILD_SIMULATION)
   if (NOT APPLE)
     o2_add_test(

Detectors/Base/CMakeLists.txt

@@ -88,6 +88,10 @@ class O2Tessellated : public TGeoBBox
   const TBuffer3D& GetBuffer3D(int reqSections, Bool_t localFrame) const override;
   void GetMeshNumbers(int& nvert, int& nsegs, int& npols) const override;
   int GetNmeshVertices() const override { return fNvert; }
+
+  /// Fill \a array with \a npoints points on this solid's boundary: every vertex, then deterministic R2 samples on facet interiors.
+  Bool_t GetPointsOnSegments(Int_t npoints, Double_t* array) const override;
+
   void InspectShape() const override {}
   TBuffer3D* MakeBuffer3D() const override;
   void Print(Option_t* option = "") const override;

Detectors/Base/include/DetectorsBase/O2Tessellated.h

@@ -32,6 +32,7 @@
 
 #ifdef O2_WITH_VECGEOM
 #include "TGeo2VecGeom/RootGeoManager.h"
+#include <VecGeom/base/Version.h>
 #include <VecGeom/management/GeoManager.h>
 #include <VecGeom/management/ABBoxManager.h>
 #include <VecGeom/management/BVHManager.h>
@@ -40,7 +41,11 @@
 #include <VecGeom/navigation/NewSimpleNavigator.h>
 #include <VecGeom/navigation/BVHNavigator.h>
 #include <VecGeom/navigation/SimpleLevelLocator.h>
+#if VECGEOM_VERSION >= 0x020000
+#include <VecGeom/navigation/SimpleABBoxLevelLocator.h>
+#else
 #include <VecGeom/navigation/BVHLevelLocator.h>
+#endif
 #include <VecGeom/navigation/VNavigator.h>
 #include <VecGeom/volumes/LogicalVolume.h>
 #include <mutex>
@@ -557,6 +562,14 @@ void GeometryManager::loadGeometry(std::string_view simPrefix, bool applyMisalig
 
 namespace
 {
+/// Volumes with very few daughters are cheaper to brute-force than to accelerate. Defined once
+/// because two places must agree on it: where the navigators and locators are attached below,
+/// and where vecGeomMaterialBudget() decides how to take a step.
+bool usesBvhAcceleration(vecgeom::LogicalVolume const* vol)
+{
+  return vol->GetDaughtersp()->size() > 2;
+}
+
 /// Converts the currently loaded TGeo geometry to VecGeom and sets up navigators, once per
 /// process, the first time the VecGeom backend is requested. Not part of loadGeometry(),
 /// which every job calls regardless of whether it ever uses the VecGeom backend.
@@ -578,16 +591,35 @@ void ensureVecGeomWorldBuilt()
     vecgeom::BVHManager::Init();
 
     // For each logical volume, set both a navigator (used for ComputeStep) and a matched
-    // level locator (used for point relocation after a boundary crossing via GlobalLocator);
-    // volumes with very few daughters are cheaper to brute-force than to accelerate.
+    // level locator (used for point relocation after a boundary crossing via GlobalLocator).
     for (auto& lvol : vecgeom::GeoManager::Instance().GetLogicalVolumesMap()) {
       auto* vol = lvol.second;
-      if (vol->GetDaughtersp()->size() <= 2) {
+      if (!usesBvhAcceleration(vol)) {
         vol->SetNavigator(vecgeom::NewSimpleNavigator<>::Instance());
         vol->SetLevelLocator(vecgeom::SimpleLevelLocator::GetInstance());
       } else {
+#if VECGEOM_VERSION >= 0x020000
+        // VecGeom 2 turned BVHNavigator into a plain class with static entry points instead of a
+        // VNavigator singleton, so there is nothing to attach: vecGeomMaterialBudget() calls it
+        // directly.
+        //
+        // The locator changes too, and not by choice. BVHLevelLocator does not compile in 2.1.0 or
+        // 2.1.1 -- the header is byte-identical in both -- because its four LevelLocate() calls
+        // have no match among the single templated BVH::LevelLocate(int exclude_item_id, ...) that
+        // v2 ships. It survived two releases because nothing in VecGeom includes that header
+        // except itself, so upstream CI never compiles it; O2 appears to be its only consumer.
+        //
+        // SimpleABBoxLevelLocator is the accelerated stand-in, using the ABBoxes built just above.
+        // Three of the four methods could be rebuilt on the templated API (the idiom is in
+        // BVHNavigator itself: bvh->LevelInside<BVHNavigator>(exclude_id, point, id, dlp)), but
+        // the direction-aware LevelLocateExclVol has no v2 counterpart at all, so this stays a
+        // fallback rather than a reimplementation. Revert to BVHLevelLocator once upstream fixes
+        // or removes it, and measure: whether ABBox location costs anything real here is unknown.
+        vol->SetLevelLocator(vecgeom::SimpleABBoxLevelLocator::GetInstance());
+#else
         vol->SetNavigator(vecgeom::BVHNavigator<>::Instance());
         vol->SetLevelLocator(vecgeom::BVHLevelLocator::GetInstance());
+#endif
       }
     }
   });
@@ -614,9 +646,24 @@ o2::base::MatBudget GeometryManager::vecGeomMaterialBudget(float x0, float y0, f
     dir[i] *= invlen;
   }
 
+  // Only the allocation differs between VecGeom versions; everything below works on pointers in
+  // both, which also keeps the std::swap() at the end of the loop a pointer swap rather than a
+  // copy of the state itself.
+  //
+  // VecGeom 1 builds NavigationState as NavStatePath, a variable-size object that must be told the
+  // maximum depth at construction and can only be made through MakeInstance(). VecGeom 2 dropped
+  // NavStatePath and MakeInstance with it: NavigationState is NavStateIndex or NavStateTuple, both
+  // fixed-size value types, so a thread_local object is the direct equivalent.
+#if VECGEOM_VERSION >= 0x020000
+  thread_local static vecgeom::NavigationState newnavstateStorage, currnavstateStorage, startCacheStorage;
+  thread_local static vecgeom::NavigationState* newnavstate = &newnavstateStorage;
+  thread_local static vecgeom::NavigationState* currnavstate = &currnavstateStorage;
+  thread_local static vecgeom::NavigationState* startCache = &startCacheStorage;
+#else
   thread_local static vecgeom::NavigationState* newnavstate = vecgeom::NavigationState::MakeInstance(vecgeom::GeoManager::Instance().getMaxDepth());
   thread_local static vecgeom::NavigationState* currnavstate = vecgeom::NavigationState::MakeInstance(vecgeom::GeoManager::Instance().getMaxDepth());
   thread_local static vecgeom::NavigationState* startCache = vecgeom::NavigationState::MakeInstance(vecgeom::GeoManager::Instance().getMaxDepth());
+#endif
   thread_local static bool startCacheValid = false;
 
   Vector3D currPoint(x0, y0, z0);
@@ -649,9 +696,18 @@ o2::base::MatBudget GeometryManager::vecGeomMaterialBudget(float x0, float y0, f
   Int_t nzero = 0;
   while (remainingDist > 1.E-10) {
     auto* lvol = currnavstate->Top()->GetLogicalVolume();
-    accountMaterial(static_cast<TGeoMaterial*>(lvol->GetMaterialPtr()), budStep);
-    vecgeom::VNavigator const* navigator = lvol->GetNavigator();
-    double step = static_cast<double>(navigator->ComputeStepAndPropagatedState(currPoint, dirr, remainingDist, *currnavstate, *newnavstate));
+    // Not LogicalVolume::GetMaterialPtr(): VecGeom 2 dropped the material slot from the logical
+    // volume. TGeo2VecGeom keeps what its conversion hook returned, indexed by logical volume id,
+    // and serves it for both VecGeom versions.
+    accountMaterial(static_cast<TGeoMaterial*>(tgeo2vecgeom::RootGeoManager::Instance().GetMaterialPtr(lvol)), budStep);
+#if VECGEOM_VERSION >= 0x020000
+    const double step =
+      usesBvhAcceleration(lvol)
+        ? static_cast<double>(vecgeom::BVHNavigator::ComputeStepAndPropagatedState(currPoint, dirr, remainingDist, *currnavstate, *newnavstate))
+        : static_cast<double>(lvol->GetNavigator()->ComputeStepAndPropagatedState(currPoint, dirr, remainingDist, *currnavstate, *newnavstate));
+#else
+    const double step = static_cast<double>(lvol->GetNavigator()->ComputeStepAndPropagatedState(currPoint, dirr, remainingDist, *currnavstate, *newnavstate));
+#endif
     if (step < 2.E-10) {
       nzero++;
     } else {